Thermal loading on crystals in an X-ray free-electron laser oscillator

X-ray free electron laser oscillators (XFELO) is future light source to produce fully coherent hard X-ray pulses. The X-rays circulate in an optical cavity built from multiple Bragg reflecting mirrors, which has a high reflectance in a bandwidth of ten meV level. The X-ray crystal mirrors exposed to intense X-ray beams in the cavity are subjects to thermal deformations that shift and distort the Bragg reflection. Therefore, the stability of the XFELO operation relies on the abilities of mirrors to preserve the Bragg reflection under such heat load. A new approach was used to analyze the heat load of mirrors and the XFELO operation. The essential light-matter interaction is simulated by the GEANT4 with a dedicated Bragg-reflection physical process to obtain the precise absorption information of the XFELO pulse in the crystals. The transient thermal conduction is analyzed by the finite-element analysis software upon the energy absorption information extract from GEANT4 simulation. A simplified heat-load model is then developed to integrate the heat load in the XFELO. With the help of the heat-load model, the analysis of XFELO operating with several cryogenically cooled diamond mirrors is conducted. The results indicate that the heat load would induce an oscillation when XFELO operates without enough cooling.


INTRODUCTION
In the hard X-ray regime, the operating free electron lasers (FELs) are based on the SASE (self-amplified spontaneous emission) mode [1][2][3][4][5], which generates X-ray pulses with unique characteristics, such as high spatial coherence, ultra-high peak power, and ultrashort pulse duration. However, the stochastic nature of SASE leads to low temporal coherence and poor pulse-to-pulse stability in the intensity of produced X-ray pulses. In order to overcome these drawbacks, an X-ray free electron laser oscillator (XFELO) has been proposed [6][7][8], where the X-ray pulses circulate in a low-loss optical cavity formed by multiple Bragg-reflecting crystals. The crystal is used because crystals could reflect hard X-rays with close to 100% reflectivity and a10 meV bandwidth [9,10], whereas conventional mirrors are unavailable in the short wavelength region.
For the cavity having a reasonable length, the repetition rate of the electron bunch has to be about 1 MHz or higher. With intra-cavity pulse energy of about 800 µJ and a beam radius of about 50 µm at a repetition rate of MHz level, the intense x-ray pulses would impose a high heat load to the mirrors, which results in the lattice distortion. The lattice distortion would decrease the peak of Bragg reflectivity, shift the wavelength range, and increase sidebands [11,12], as the Bragg reflection originates from the X-rays scattering of a periodic atomic lattice. Therefore, the realization of the optical cavity with Bragg reflecting mirrors relies highly on the stability of the crystal lattice. Thus, the heat load of the crystal mir- * denghaixiao@zjlab.org.cn ror under such intense X-ray pulses is an essential issue for building an XFELO.
In order to investigate the influence of heat load in the crystal , two types of experiment methods are used to test performance limitations of the crystal thermal deformation [13][14][15][16][17][18][19]. The first one is to utilize a high power conventional laser in the long-wavelength regime for modelling laser-mirror interactions. The advantage of this method is the possibility of precisely controlling the pulse duration and pulse energy, which is crucial for lightmaterial interactions. However, this approach produces a penetration depth at the nanometer level [14], which is far shorter than those produced by X-rays. Besides, the repetition rate of the high power conventional laser is challenging to reach MHz or higher. Another method is to apply the synchrotron radiation light sources for the heat load experiments. The synchrotron radiation could reach the hard X-ray regime with enough timeaveraged power. However, the synchrotron radiation is polychromatic, and the monochromator is required to obtain high spectral purity. With the monochromator, the pulse energy of the synchrotron radiation light source decreases significantly to the level that is far lower than XFELO [15]. Besides, if the X-ray beam is micro-focused to the FEL beam level [18], its angular spreads would be much larger than Bragg reflection acceptance. Thus, in general, the experimental conditions do not meet the requirements of XFELO operation, and it is still necessary to numerically model the heat load of mirrors and its feedback on XFELO operation.
In previous studies [20][21][22], the numerical evaluation of the X-ray absorption in the crystal uses an exponential attenuation model that describes the intensity of the X-rays decreases to 1/e of initial value through one attenuation length [23]. The exponential attenuation model is practicable in most of the X-ray interactions but not satisfied in the X-ray Bragg reflection [24]. Instead, a simulation model describing the physical process of Bragg reflection in GEANT4, a powerful particle tracker, is developed in this paper. Combining the other implemented light-matter physical process in GEANT4, a comprehensive set of interactions is provided into the currently simulating processes. As the heat source in the crystal is obtained, a start-to-end like transient thermal analysis can be conducted with the finite-element analysis (FEA) software [25]. FEA simulations can determine the temperature map of the crystal, and therefore the strain field in a heat-distorted crystal. Furthermore, by coupling the thermal loading into the XFELO simulation, it can gain the ability to investigate the continuous influence of crystal thermal loading on XFELO operation. This paper is organized as follows. The Sec. II presents the characteristics of Bragg reflections relating to the thermal load. In Sec. III, the light-matter interaction is simulated within GEANT4, and the following thermal behavior is analyzed. In Sec. IV, a simplified model to couple the thermal loading in XFELO simulation is presented. Then the simulation of thermal loading coupled XFELO is shown. A summary is presented in Sec. V.

II. CHARACTERISTICS OF BRAGG REFLECTIONS
An X-ray free electron laser oscillator is a low-gain device, in which the X-ray pulse is amplified by a gain G for each round trip when it overlaps with an electron bunch in the undulator. In addition, the X-rays circulate in a low-loss optical cavity formed by multiple Bragg mirrors, corresponding to a total reflectance R. In general, when (1 + G)R > 1, the X-rays could evaluate from the initial noise to intense coherent radiation. And the oscillator reaches steady state when (1 + G)R = 1.
Overall, the Bragg reflectivity should be sufficiently high so that the initial exponential gain of the intracavity pulse energy can be sustained for a reasonable set of the electron beam and undulator parameters. As the single-pass gain is fixed by the electron beam qualities and undulator configuration, the ability to achieve the stable XFELO operation relies on the reflectivity R that is sensitive to the lattice structure of the crystals. Adverse effects, including the thermal expansion of the crystal and mounting vibration, may emerge to affect the lattice and to disrupt the Bragg reflection. Thus, the relationship between the distortion of the crystal lattice and the Bragg reflection is discussed because of the unique characteristics of the Bragg reflection.
The Bragg reflection or diffraction originates from the X-rays scattering of atoms in a periodic structure [26], which can be imagined as a reflection of imaginary "mirrors" formed by atomic planes in the crystal lattice, seen in Fig. 1. The Bragg condition can be expressed as where n is the diffraction order, θ B is the angle of incidence, λ B is the radiation wavelength at which the reflection occurs, d H is the distance between the atomic planes and T the temperature. The lattice planes is described by the Miller indices, h, k, l, three integers. The spacing between the lattice planes is given by where a is the lattice constant. For a diamond crystal, a cubic lattice, the lattice constant a is 3.567 [27]. Diamond is the preferred material employed by a high power XFELO to form the X-ray cavity. Because it is a unique combination of outstanding thermal and optical properties, including a high thermal diffusivity, a low thermal expansion, and a high Bragg reflectivity of x-rays. In the following discussions, the crystal mirror is assumed to be diamond for the design and the optimization of the XFELO configuration. As expected, the high reflectivity in Bragg reflection is only possible if the lattice structure is perfect for a large amount of diffracting plane. In general, the change of crystal lattice acts on both the position and the amplitude of the Bragg reflectivity [11,12]. Fig. 2 shows three kinds of lattice distortions (up panel) and corresponding effect (bottom panel). In each condition, the undistorted lattice (black dot) is compared with the strained lattice (red dot). The lower plots present the influence of these distortions on the reflectivity curve over the spectrum. a) Random disorder decreases the peak reflectivity of the Bragg reflection. b) Homogeneous linear strains shift the reflection. c) Longitudinal non-linear strains create sidebands and decrease the peak value of the Bragg reflectivity.
When a high-brilliance X-ray pulse is reflected by a Bragg reflecting crystal, a certain part of the energy is deposited. The energy from light is first absorbed by the electrons and then transferred to the crystal lattice. As the position and the amount of deposited energy are different, various thermal distortions of the lattice would be presented in the crystal. The type a) distortion occurs for a very high energy density that imposes radiation damage upon crystal to lose the long-range order, e.g., ionization. Far below the damage threshold, the absorbed energy due to the light-matter interactions translate into a thermal strain that leads to the other disorders, type b) and c), which are the major issues in the XFELO heat load. The deformed reflectivity curve of the large thermal deformation can be obtained by solving the Takagi-Taupin equations with the strain profile owing to the X-ray heating [28]. In principle, the relative change of the lattice constant can be described by the thermal expansion coefficient, β. For a cubic crystal like the diamond, the relative change owing to the thermal expansion can be expressed as: While the temperature change is small, β can be assumed constant and relative change of lattice constant can be approximated to a(T )−a(T0) with the first order Taylor series. The thermal expansion coefficient β(T ) can be obtained by an empirical formula from Ref. [27], which can be written as β = 4.25 × 10 −14 T 3 .
As equation (3) predicted, a small expansion coefficient is expected to reduce the change of the crystal and to increase the tolerance of the heat load on the crystal. It is therefore preferable to operate the mirrors at a low temperature, below 100 K. Because in the cryogenic temperature, the thermal expansion coefficient drops by a few orders of magnitude compared with the room temperature value of ≈ 1 × 10 −6 . Besides, the thermal conductivity of diamond would reach the peak value in the temperature range from 60 K to 100 K [29].

III. LIGHT-MATTER INTERACTIONS AND HEAT CONDUCTIONS
In order to investigate the lattice distortion of mirrors under heat load, it is fundamental to obtain the accurate and reliable information of the absorbed part of the FEL X-ray pulse. This means that the interaction between the XFELO radiation and the crystal mirror is an essential limitation in the realization of thermal loading calculation. In this section, with the help of a dedicated Bragg reflection model, the light-matter interactions have been numerically simulated by using the GEANT4 toolkit [30], which is developed and actively used for particle experiments and detector designs.
GEANT4 is a Monte Carlo based toolkit for simulating the particle interacting with matter. It includes many functionalities like tracking, geometry, physics models, and hits. For convenience, it already implements many widely used physics models to describe the fundamental light-matter interactions. Additionally, many other physics processes and low-energy models are included, covering a wide range of interactions, e.g., electromagnetic, hadronic, and a large set of long-lived particles.
Based on the Low Energy Livermore Model in the GEANT4 [31], an additional Bragg reflection process has been implemented to extend the available range of particle interactions. The implemented Bragg reflection model could change the direction of the photon in the crystal, which is satisfied with the Bragg condition. Thus, this model naturally enjoys the future of the multiple reflecting, which is the most important consideration for Bragg reflection and is absent from the previous exponential attenuation model. The crystal is built in geometry construction with its orientations, lattice information. As the Bragg reflection occurs in a range of ten meV level, the table that stores the cross section of each interaction in the GEANT4 kernel must be extended to reach the meV level. The cross section of the Bragg reflection, the most important factor in the GEANT4 simulation, is calculated with the help of the dynamical theory of X-ray diffraction. The cross section is given by where 1/µ is the excitation length calculated by the dynamical theory [9,32], m u = 1.66053886 × 10 −24 g is the atomic mass unit, and A is the relative atomic mass of the target element. Beyond the Darwin width, the Bragg reflection filter out the spectral components. Therefore, it is assumed for simplicity that the cross section beyond the Bragg width is zero in the GEANT4 simulations. As an example, the Bragg width is about 11 meV while the photon energy is 14.33 keV at the diamond (3 3 7) reflection. Besides, a parameter defined as the probability of the reflected X-ray to be reflected again can be used to control the multiple Bragg reflection and to correct the reflectivity in a very thick crystal so as to satisfy the predictions of the dynamical theory. According to the earlier simulation results of XFELO mode for SHINE (Shanghai HIgh repetition rate XFEL aNd Extreme light facility), the first hard X-ray free electron laser in China [33,34], a set of baseline parameters are utilized for the generation of primary photons in the GEANT4 simulation. The cavity is formed by four diamond mirrors with (3 3 7) diffracting planes while the cavity configuration is optimized to amplify the X-ray pulse with a photon energy of 14.33 keV. The intra-cavity pulse energy is assumed to be 600 µJ with a 1 MHz repetition rate [35]. Due to the excellent longitudinal coherence, constant phonon energy of 14.33 keV is adopted. The temporal distribution of X-ray photons is assumed to a Gaussian distribution with an RMS length of 12 µm. The total length of the X-ray pulses, defined as the distance from the first photon to the last one, is fixed to be 100 µm in the photon generator. Meanwhile, the pulse has a Gaussian transverse profile with an RMS size of 25 µm, while the divergence angle is assumed to be 0.7 µrad, which is much smaller than Bragg reflection widths. In addition, the crystal size is 1000µm×1000µm×70µm in the GEANT4 simulation.
Overall, in the view of implemented and defined components, all aspects of the Bragg reflection simulation have been involved: 1) the geometry of the diamond crystal, 2) the primary particle (photons), 3) the generation of particles, 4) the tracking of particles through materials, 5) the physics processes modeling particle interactions (photoelectric effect, Compton scattering, pair production, bremsstrahlung, and Bragg reflection.). Fig. 3 shows a snapshot of the particle tracking in the crystals. As it is shown, the Bragg reflections govern the interaction of the X-rays with the crystal. The absorptions of the X-rays mainly result from the photoelectric effect and Compton scattering. The Bragg reflectivity could be optimized by adjusting the thickness of the crystal. According to the statistical data from the GEANT4 simulation at 14.33 keV, a diamond mirror with a thickness of 70 µm allows reflecting 87% of the incident X-ray while the transported part remains 11.7%. In other words, only 1.3% of the incident X-ray has been absorbed. Since the dynamical theory of X-ray diffraction predicts the reflectance to be 88% and transmittance to be 11.7%, the GEANT4 simulation results are reasonable and acceptable. Consider an intra-cavity X-ray FEL pulse with a pulse energy of 600 µJ at 14.33 keV reflecting by a diamond, the total absorbed energy is about 7.3 µJ in the diamond crystal during the reflection. For the diamond (3 3 7) reflection at 14.33 keV, the Bragg angle is 83 • nearing the exact backscattering. In this condition, the transverse distribution of the absorbed energy is mainly fixed by the transverse profile of the incident X-ray pulse, seen in Fig. 4. With an RMS spot size of about 25 µm, the maximum energy density reaches 2400 µJ/mm 2 corresponding to time-averaged power of 2.4 kW/mm 2 at 1 MHz. Since the crystal rotates about 7 • along the x-direction to fulfill the Bragg condition, the irradiation spot size in Fig. 4 on y-direction is a little larger than that of x-direction. Because the reflection of a non-normal incidence would cause a long tail with relating to transverse shift distance as radiation penetrated in crystal. For a small Bragg angle, this effect would significantly affect the transverse distribution of the absorbed energy. Moreover, it is insufficient to evaluate the multiple reflections in a theoretic formula but is simple in this GEANT4 model.
One advantage of this model is the ability to obtain the three-dimensional absorption information of the Xray pulse while Bragg reflection is considered. The distributions of the absorbed energy on x-z plane (top) and y-z plane (bottom) are shown in Fig. 5. The incident direction of the X-rays points to the negative value of the z-axis. The tilted penetration satisfying the Bragg condition is demonstrated in the bottom plot of Fig. 5. Due to the short excitation length related to a high reflectance, the incident X-ray pulse is attenuated significantly. Since it can give more details about the absorbed energy, this model could be a powerful tool to investigate the thermal loading of XFELO. Another advantage is that this solution could study dynamical absorption, the absorbed energy at different times, seen in Fig. 6. In this case, the time of zero is the time when the first photon enters crystal. At the time of 500 fs, most of the interactions of X-rays with crystal are finished. The time as light-matter interaction finished is determined by the X-ray pulse duration and the crystal thickness. The peak value of the absorbed energies is inwards crystal, as the absorption always occurs after penetration. Unlike the exponential attenuation model, which predicts a near-uniform trend for the thickness of only 70 µm and a large attenuation length of the diamond, this model indicates that the absorptions concentrate on the part near the incident surface of the crystal owing the high reflectance. For a thinner crystal, the X-ray would penetrate the crystal with few multiple reflections, and the absorptions are quite similar to the prediction of the exponential attenuation model. Further, with the help of molecular dynamics simulation software, the dynamical absorption can be extended to investigate the possible dynamical change of the crystal over the time of the pulse duration. Since the XFELO pulses may be too short on the time scale to influence their own Bragg reflection [36,37], one can assume that the crystal temperature change acts only on the subsequent pulses at present research. Therefore, the crystal temperature seen by the next pulse is much more crucial in order to maintain the stability of the cavity. Additionally, the time between sequential pulses is far longer compared to the response time of crystal thermalization. Thus, it can be assumed that the X-ray pulses induce instantaneous heating in the crystal. In other words, the energy depositions resulting from lightmatter interactions translate into the transient temperature increase.
In order to study the thermal behavior of the diamond crystal after the intense X-rays irradiation, the commercial finite element analysis software ANSYS is used [25]. The initial temperature is obtained from the absorption information of the GEANT4 simulation. In ANSYS transient thermal analysis, the crystal size is 800µm×800µm×70µm, and the environmental temperature is 70 K. The work in Ref. [29,38,39] is adopted to obtain the thermal conductivity and specific heat capability, and peak value of thermal conductivity can be found between 60 K and 100 K. As it is various as a function of the temperature, a fitting linear polynomial formula is utilized for simplicity. The mesh size is optimized, and the time step is fixed to be 0.5 ns.
Calculated results of the heat conduction are presented in Fig. 7 and Fig. 8. Fig. 7 shows snapshots of the temperature filed at various time of 1, 3, 5, 10, 250 and 1000 ns after irradiation. Fig. 8  As shown in Fig. 7, the initial temperature is higher on the incident surface than on the back surface due to the Bragg reflection. The highest temperature is about 1300 K near the center of the incident surface. It drops rapidly in the first 50 ns, and then gradually slow in following 100 ns. Because the thermal conductivity increases rapidly when the temperature decreases in the first 50 ns, and then heat diffusion is slowing when the temperature gradient decreases gradually. The thermal relaxation time across the whole thickness is considered to be about 20 ns, as suggested by the nearly identical temperature on the incident and back surface in Fig. 7 and Fig. 8. Meanwhile, the thermal relaxation time in the radial direction is over 150 ns, as the temperature on the incident and back surface reaches the volume-averaged temperature. Indeed, the temperature averaged over the whole crystal is nearly homogeneous after 200 ns. The maximum temperature difference over the entire crystal is about 0.1 K at 250 ns, and it reduces to 0.05 K at 500 ns. In fact, the averaged temperature increase to about 74 K, with a total absorbed energy of 7.3 µJ. In the following 1000 ns, the averaged temperature decrease to about 71.8 K, corresponding to a relative change of lattice parameter of 1 × 10 −8 that meets the requirement of cavity stability. At 2000 ns, the averaged temperature decreases to 0.6 K. It still have a downward trend but would continue to decrease over a longer period.
The multiple-pulse result is present in Fig. 9. All the parameters except time step are the same as those in the above case, and the energy of individual laser pulses is added periodically. The time step is fixed to be 1 ns. The simulations are undertaken over 20 µs, 20 pulses at 1 MHz. The averaged temperature rises instantaneously upon each pulse. Its evolution over the interpulse period is similar to that in the above single-shot case. However, the averaged temperature keeps on rising, for some residual heat is left behind in the crystal and the finite thermal conductivity. Thus, the temperature increment consists of two parts: 1) the nearly instantaneous heating due to the X-rays irradiation; 2) long-term accumulation of the remaining heat in crystal. Since heat transmission depends on the temperature difference, the rising temperature would improve the capability of transferring heat to the cooling system. Thus, the remnant tempera-7LPH>µV@ 7HPSHUDWXUH>.@ FIG. 9. Time histories of the averaged temperature (minus 70 K). The repetition rate is 1 MHz, and the periodic input laser pulse is the same as the single-pulse case.
ture drops against temperature increase, and the system would eventually evolve to a steady-state with a stable accumulated temperature increment. In our case, the temperature increment averaged in the whole crystal is about 4 K, and the accumulated one at the steady state is about 3.5 K. The heat accumulation raises a practical requirement for the cooling system that ensures the accumulated temperature is trivial. The solution should be established in a specific case where an XFELO operates, and it is not the purpose here. The specific cooling system design can take part experiences from the monochromator design of the synchrotron radiation light sources.

IV. THERMAL LOADING COUPLED XFELO SIMULATION
Since it is complicated to run the time-dependent XFELO simulation considering the thermal loading, a simplified model to describe the thermal loading is needed to be developed. The thermal loading mainly affects the Bragg reflection through the distortion of the lattice. As mentioned in Sec. II, the pulse energy determines that the type 1) distortion is negligible in XFELO. Furthermore, based on the transient thermal simulation above, the type 3) distortion that needs a temperature gradient across the thickness could also be negligible. The thermal relaxation time across the thickness at a cryogenic temperature is very short compared to the interpulse time. Thus, in the long-term accumulation the temperature gradient across the thickness is negligible, or even the next pulse would see the temperature distribution is nearly uniform upon the entire crystal. Thus, it can be assumed that the most crucial issue of the thermal loading is the shift of the reflectivity curve due to the homogeneous expansion through the crystal thickness.
Therefore, the simplified model could focus on the temperature change averaged over the whole crystal or over the effective volume of the irradiation. It is reasonable because the variation of the temperature inside the diamond is negligible when the thermal relaxation time is much shorter than the cooling time. This assumption is based on the fact that the thermal conductivity of diamond at cryogenic temperature is very high. Thus, a simplified model of the averaged temperature can be expressed as [40]: where ∆T = T −T 0 is the relative change of the averaged temperature, h ef f is the heat transfer coefficient related to the design of the cooling system, S is the surface area of heat transfer, m is the mass of diamond and C p is the specific heat [39]. When the temperature does not change significantly, the value of h ef f S mCp can be assumed to be a constant. Then the solution can be written as where α = h ef f S mCp . With a stable repetition rate, the time for cooling is a constant. Thus, we can define η = e −α∆t as the cooling efficiency at the repetition rate of 1/∆t. It describes how much temperature is left behind after a time of ∆t. Its approximation value can be obtained from the transient thermal analysis of single-pulse input.
Using the parameter η, the averaged temperature change of multiple-pulse input can be written as: where n indicates nth pulse. At steady state, the averaged temperature change become stable and we can get a solution, ∆T n+1 = ∆T n = η∆T p /(1 − η). In the singlepulse thermal analysis, η is about 1.8/4 = 0.45. And the steady-state temperature change is 3.27 K, which is a little lower than that produced by multiple-pulse thermal analysis but reasonable for the simplicity. Then, the shift of reflectivity can be written as ∆E ≈ E H β∆T n , where E H is the Bragg energy.
With the help of the simplified thermal loading model, an XFELO operation of SHINE coupled with heat load is conducted. The parameters of SHINE is presented in Tab. I. The main accelerator of SHINE uses the superconducting technology to produce electron bunches with the energy of 8 GeV, a 100 pC charge compressed to the peak current of 700 A, and 1 MHz repetition rate.
The optical cavity is built from four diamond (3, 3, 7) crystal mirrors, seen in Fig. 10. With a thickness of 70 µm, the downstream diamond mirror C 4 at the Bragg energy of 14.33 keV reaches 87% (coupling output is 12%). The other three mirrors is expected to reach the total reflection with thick crystals. To simplify the calculation, the temperature of the three crystals is appropriate to be fixed by 70 K. The heat load is mainly considered in the downstream mirror C 4 due to its thinner thickness. The simulations are conducted by using the combination of a time-dependent FEL code GENESIS [41], a field propagation simulation code OPC [42], and a Bragg reflection simulation code BRIGHT [43]. The model of simplified thermal loading is integrated into the BRIGHT as a part of Bragg reflection. In the following analysis, the X-ray pulse is indicated to be the intra-cavity one to present how much energy is induced to the mirrors. The output value needs to be multiplied by an outputcoupling factor, the transmittance.
In order to investigate the influence of the thermal loading, the simulations were conducted while η is 0.1, 0.5, 0.75, and 0.95. The results are presented in Fig. 11 and Fig. 12. When η = 0.95, the reflectivity shift of the downstream mirror can be observed, shown by a green dash line in the top plot of Fig. 11. In comparison, the Bragg reflectivity of the upstream mirror is shown by an orange dash line. A typical spectrum of the XFELO pulse at η = 0.95 is shown by a red line. As expected, the spectral intensities exist in the overlap of the two reflectivity curves from different mirrors. The total integrated reflectance changes dramatically while the shift of the Bragg reflectivity position exceeds half of the Darwin width. This effect would significantly influence the stability of optical cavity and gain process of the XFELO. Thus, the peak power decreases significantly as η rises, as shown in the bottom plot of Fig. 11. Another fact should be noted that the different reflectivity also results in a different time delay caused by Bragg reflection. The mismatch of time delay decreases the peak power further, for the X-ray pulses overlap on the electron beam imperfectly. It would be a crucial issue for electron beams with a short bunch length. Fig. 12 shows the detailed simulation results, including the spectrum evolution and pulse energy growth with each η. The left panel of Fig. 12 presents the energy evolutions (red line) and its corresponding temperature change (green line). With a small η (0.1), which gives an insignificant temperature change, the energy of the x-ray pulse grows as it is amplified by electron beams again and again. The pulse energy reaches 1 mJ at saturation, and the corresponding temperature change is only 0.8 K. At the steady-state, pulse energy has a very small oscillation, which is an XFELO nature. When η is set to 0.5, the saturated pulse energy decrease to about 900 µJ, and the temperature change increase to 10 K. The oscillation is also enlarged, as the reflectivity curve is more unstable. The decrement of the pulse energy results from the mismatch of time delay and the relatively large cavity 5RXQGWULS (QHUJ\>µM@ 7HPSHUDWXUH>.@ loss due to the shift of the reflectivity curve. While η increases to 0.75, the large remaining heat would lead to a non-trivial temperature change that in-terrupts X-ray pulse growth through the reflectivity shift. When the total reflectivity can not maintain the stability of the system, the pulse energy decreases, and the temperature change follows it. However, when the temperature change is below the critical value, the pulse energy grows exponentially again as a comparatively intense seed is retained. The repetition of this process leads to that the temperature gradually approaches the critical value and that pulse energy has an attenuating oscillation. The steady-state pulse energy is about 350 µJ, only a third of the value of η = 0.1.
A large oscillation occurs when much heat remains in the crystal with a large η. The oscillations are attenuated but stay at a relatively large value eventually. The root of the oscillation is the lingering negative feedback of the temperature change on the pulse energy. Besides, the temperature change is always behind the energy increase. The combined result is that temperature crest gradually approaches the troughs of the pulse energy and reach a balance where oscillation occurs.
The evolution of the spectrum is presented in the right panel of Fig. 12. The simulations share the same initial random number. Thus, the result for each η has nearly the same period of struggling to build up the longitudinal coherence of the X-ray pulse. As energy increases, the spectral width decreases further. The final FWHM spectral width is nearly 13 meV when η is 0.1 or 0.5. This value decreases to about 9 mev at η = 0.75 because the overlap of the reflectivity between the upstream and downstream mirror is narrower. The maximum available shift is about 6 meV, which is nearly half of the Darwin width.

V. CONCLUSION
In this paper, we have performed a numerical study to modelling heat load of diamond mirrors in the XFELO operation, including the light-material interaction, the thermal behavior analysis, and the corresponding heatload coupled XFELO simulation. The light-material interactions were simulated by the particle tracking package GEANT4 with the dedicated Bragg reflection physical process. The advantage of this GEANT4 model is to give accurate and reliable predictions of X-ray absorptions compared to the theoretic attenuation model. The following transient thermal conduction was analyzed by the commercial finite-element analysis software AN-SYS. With the 70 K initial temperature and crystal size of 800µm×800µm×70µm, cases for single and multiple pulse input are analyzed. The single-pulse input result indicates that the temperature gradient along the thickness would be negligible in 20 ns, and the temperature would be nearly homogeneous on the whole crystal in 150 ns. The multiple-pulse result indicates that the accumulated temperature increment would reach 3.5 K, which meets the requirement of mirrors for XFELO operation at 70 K in principle. To further investigate the impact of the heat load on XFELO operation, a simplified thermal loading model focusing on the temperature change averaged over the entire crystal is described. Within the simplified model, the parameter η is defined to describe the ratio between remaining and initial temperature change. Then, heat-load coupled XFELO simulations are conducted with the parameters of SHINE, while η is 0.1, 0.5 0.75, and 0.95. We study the generation of 14.33 keV X-ray pulse, whose energy was found to be 1 mJ in the cavity. With a large η, the pulse energy and peak power were found to decrease markedly due to the reflectivity shifting and time delay mismatching. The simulations also indicate that the heat load would induce the oscillation of the pulse energy and spectrum at a steady state. And the scope of the energy oscillation is determined by the heat load and the value of η.
Beyond this work, the cooling system of mirrors for the XFELO operation of SHINE will be designed and optimized. With the optimized cooling system, the thermal loading analysis with non-Fourier heat conduction will be conducted, so as to provide a better understanding of various heat load effects.