research papers
Time and space resolved modelling of the heating induced by synchrotron Xray nanobeams
^{a}ESRF – The European Synchrotron, 71 Avenue des Martyrs, 38000 Grenoble, France, ^{b}Department of Physics, University of Torino, via P. Giuria 1, 10125 Torino, Italy, ^{c}Istituto Nazionale di Fisica Nucleare, Sezione di Torino I, 10125 Torino, Italy, ^{d}Department of Applied Science and Technology, Politecnico di Torino, Torino 10129, Italy, ^{e}Institut Sciences Chimiques de Rennes, UMRCNRS 6226, Campus de Beaulieu, Université de Rennes 1, 35042 Rennes Cedex, France, and ^{f}Department of Chemistry and Interdepartmental Centre NIS, University of Torino, via P. Giuria 7, 10125 Torino, Italy
^{*}Correspondence email: valentina.bonino@esrf.fr
Xray synchrotron sources, possessing high power density, nanometric spot size and short _{2}Sr_{2}CaCu_{2}O_{8+δ} microcrystals by means of nanopatterning experiments with intense hard Xray nanobeams is calculated. It is demonstrated that the temperature evolution is strictly connected to the filling mode of the storage ring. By coupling the MC with the heat equation, Xray pulses that are 48 ps long, possessing an instantaneous of ∼44 × 10^{13} photons s^{−1}, were found to be able to induce a maximum temperature increase of 42 K, after a time of 350 ps. Inversely, by ignoring the energy redistribution calculated with the MC, peaks temperatures up to hundreds of degrees higher were found. These results highlight the importance of the energy redistribution operated by primary and in the theoretical simulation of the Xray heating effects.
are extending their application frontiers up to the exploration of direct matter modification. In this field, the use of atomistic and continuum models is now becoming fundamental in the simulation of the photoinduced excitation states and eventually in the triggered by intense Xrays. In this work, the Xray heating phenomenon is studied by coupling the Monte Carlo method (MC) with the Fourier heat equation, to first calculate the distribution of the energy absorbed by the systems and finally to predict the heating distribution and evolution. The results of the proposed model are also compared with those obtained removing the explicit definition of the energy distribution, as calculated by the MC. A good approximation of experimental thermal measurements produced irradiating a millimetric glass bead is found for both of the proposed models. A further step towards more complex systems is carried out, including in the models the different time patterns of the source, as determined by the filling modes of the synchrotron storage ring. The two models are applied in three prediction cases, in which the heating produced in BiKeywords: oxides; radiation damage; Monte Carlo method; finite element method; Xray nanopatterning.
1. Introduction
Since the firstgeneration synchrotrons, Xray source characteristics have evolved under the strong demand for higher and higher spatial and temporal resolution. Today, peak brilliances up to 10^{26} photons s^{−1} mm^{−2} mrad^{−2} (0.1% bandwidth)^{−1} are achieved with thirdgeneration synchrotron sources, and Xray beams with nanometric spot sizes of ∼50–100 nm and pulses of hundreds of photons s^{−1} with energies in the pJ range are routinely achieved (Mino et al., 2018; MartínezCriado et al., 2016). However, the downside of this ongoing evolution of Xray sources is that such beam characteristics can exceed the threshold where photonflux density can affect the organization of matter (Bras & Stanley, 2016). As a consequence, the displays of synchrotron Xray beams inducing modification phenomena in inorganic materials, such as phase transitions (Adriaens et al., 2013), reduction (Stanley et al., 2014) and crystallization (Feldman et al., 2009; Martis et al., 2011), are now increasing in number. The Xray heating effect is frequently claimed as a possible modification mechanism. So far, however, only a few experimental attempts have been made to clarify this point. For example, in two thermalimaging experiments (Rosenthal et al., 2014; Snell et al., 2007) carried out on glass beads and on an indium microsized particle, Xray heating up to ∼30 and 0.2° have been recorded, respectively.
Although more experimental efforts are needed in this direction, new research fields are emerging in which Xray probing and material modification are strictly correlated (Tu et al., 2017; Hsu et al., 2015; Bonino et al., 2019; Pagliero et al., 2014; Truccato et al., 2016; Mino et al., 2017, 2019). In this background, we have recently demonstrated the feasibility of a new patterning concept, based on a resistfree method, which exploits the structural and electrical modifications induced in condensed matter by highintensity nanometric beams from thirdgeneration synchrotrons (Bonino et al., 2019; Pagliero et al., 2014; Truccato et al., 2016; Mino et al., 2017, 2019). For example, by irradiating the hightemperature superconductors Bi_{2}Sr_{2}CaCu_{2}O_{8+δ} and YBa_{2}Ca_{3}O_{7−δ} with a 17 keV nanobeam with a of ∼10^{11} photons s^{−1}, a change in the electrical behaviour was related to an increase in the crystal mosaicity and a depletion of the doping oxygen content (Bonino et al., 2018, 2019). Following these findings, we exploited these local modifications to nanopattern three Bi_{2}Sr_{2}CaCu_{2}O_{8+δ} microcrystals, finally fabricating electrical devices based on the Josephson effect (Pagliero et al., 2014; Truccato et al., 2016; Mino et al., 2017).
Despite all the experimental displays of photoinduced phenomena, owing to the large variety of Xray–matter interaction mechanisms that can be accessed with highintensity Xrays, a clear picture of the mechanisms involved in material modification is still missing. Depending on the beam intensity and et al., 2001; HauRiege, 2012). This is not surprising; indeed similar concepts have been extensively investigated in the field of highpower laser irradiation, where both thermal and nonthermal effects are currently exploited in many laserbased technological processes (Mirzoev et al., 1996; Liu et al., 1997).
both reversible and irreversible matter states can be induced in a variety of inorganic materials (LondonThe similarities between highintensity synchrotron Xray and optical and IR laser irradiations are several. As in femtosecond lasers, Xrays ionize the matter, but, differently from them, they excite the inner shells of the atoms and penetrate indepth into the material volume. In the femtosecond scale, the absorbed Xray energy is dissipated through the emission of fluorescent photons and Auger and photoelectrons. Then a cascade of excitation processes is initiated by electrons and eventually, in the picosecond timescale, they thermalize firstly among themselves and then with the lattice. Another difference with lasers concerns the surface power densities that can be delivered with a single pulse. Femtosecond lasers can reach 10^{22} W cm^{−2} (Bahk et al., 2004), while the picosecond pulses of thirdgeneration synchrotrons can only achieve values in the order of 10^{14} W cm^{−2}. For all of these reasons, although many concepts valid for highpower lasers can be extended to the Xray regime, differences in thermal distribution are expected and a simple extrapolation from laser data will not work. In this sense, more efforts are necessary from the theoretical point of view to fully understand the effects of highintensity Xray–matter interactions.
Both space and time scales are important in determining the simulation method to numerically approach the problem. et al., 2009; Neutze et al., 2000; London et al., 2001). As an example, by considering the Seitz–Koehler (Seitz & Turnbull, 1956), which accounts for the atom knockon by photoelectrons, in a previous study we applied the MC to simulate the displacement of the oxygen doping content to nonactive doping positions in Bi_{2}Sr_{2}CaCu_{2}O_{8+δ} microcrystals. This mechanism was found to be principally responsible for the photoinduced resistivity change for a fraction of the investigated experimental cases (Torsello et al., 2018). For larger scales, once the electrons have thermalized, continuum models are preferred to simulate effects such as melting and cracking (Nicholson et al., 2001; Mino et al., 2017; Wallander & Wallentin, 2017; London et al., 2001). As an example, by means of the finite element method (FEM) we solved the heat equation for an Xray nanopatterning experiment, calculating the temperature modulation generated by synchrotron pulses (Mino et al., 2017). Similarly, Wallander & Wallentin (2017) modelled the temperature increase induced in an InP nanowire. However, in these cases the derived heating was calculated implementing a heat source implicitly defined as the product between the beam profile and the decaying profile of Xrays as defined by the attenuationlength coefficient, i.e. without taking into account the details of the fraction of absorbed photons and of the spatial distribution of the absorbed energy.
and Monte Carlo methods (MCs), simulating the motion of atoms, ion and molecules, and the transport of photons and electrons, respectively, are applied in systems with sizes up to the micrometric range for fast processes with a maximum time limit in the picosecond scale (GnodtkeIn the present work, by continuing our efforts towards the simulation of the heating effects induced by intense Xrays, we want to assess the temperature increase in both space and time in more detail. For this purpose, by considering for the first time both the spatial distribution of the absorbed energy and the real time structure of synchrotron radiation, a more accurate numerical model has been formulated in which the heat source is modulated in both space and time. This was achieved by coupling the radiationparticle transport problem, solved with the MC, with the timedependent heat equation. With such a detailed model, the temperature increase can be calculated with an unprecedented accuracy. For comparison, the heat equation is also solved implicitly, defining the heat source with the aforementioned method already adopted in the literature.
In the following, the model validation and accuracy are first discussed by modelling the first thermalimaging experiment reported in 2006 by Snell et al. (Snell et al., 2007). The temperature transient, experimentally monitored by the authors, is compared with that predicted by the model. Then, three examples of Xray nanopatterning on Bi_{2}Sr_{2}CaCu_{2}O_{8+δ} microcrystals (Pagliero et al., 2014; Truccato et al., 2016; Mino et al., 2017) are considered to discuss the role of heating on the structural and electrical modifications observed in the fabrication process. For this purpose, the stationary, transient and pulsed regimes of the thermal model are separately treated.
2. Method
2.1. Experimental methods for validation and prediction
Based on their scope, the cases considered can be divided into two categories: validation and prediction. In the validation case, the irradiation of a glass bead (composition in mass percentage: 66% SiO_{2}, 15% Na_{2}O, 7% CaO, 5% Al_{2}O_{3}, 3% B_{2}O_{3}, 2% ZnO, 1% K_{2}O and 1% MgO) (Snell et al., 2007) at ambient temperature is implemented. In Table 1 the main experimental parameters are summarized. The modelling method used for this case is based on the definition of a range of validity in which the solutions can lie. The range is defined by varying the experimental parameters within their experimental error (reported in Table 1). The complete experimental method can be found in the reference work (Snell et al., 2007). Conversely, the three Xray nanopatterning experiments fall within the prediction category. In these cases, no experimental measurements of the temperature are given. Therefore, the related simulations have the unique purpose of estimating the temperatures reached during the patterning process. These experiments were directly performed by the authors at the European Synchrotron Radiation Facility (ESRF, Grenoble, France). The experimental details are reported in the following.
To ensure continuity with the previous works, the Bi_{2}Sr_{2}CaCu_{2}O_{8+δ} samples are named following the same sample nomenclature: WBVB08 (Mino et al., 2017), WBAP13 (Pagliero et al., 2014) and WBAP14 (Truccato et al., 2016). (i) WBVB08 was irradiated at the ID16B beamline using a pink beam (ΔE/E ≃ 10^{−2}) and a 16bunch filling mode of the storage ring with a maximum current of 90 mA. This filling mode consists of a train of 16 equally spaced highly populated electron bunches, each of them generating a current equal to 90 mA 16bunch^{−1} = 5.625 mA bunch^{−1} and a Gaussian time profile for the corresponding Xray pulse with a root mean square (RMS) duration of RMS_{t} ≃ 48 ps. (ii) WBAP13 was irradiated at the ID22 beamline with a monochromatic beam (ΔE/E ≃ 10^{−4}) with 16bunch filling mode. (iii) WBAP14 was irradiated at the ID16B beamline using a pink beam (ΔE/E ≃ 10^{−2}) with a 7/8+1 filling mode of the storage ring. This mode consists of filling 7/8 of the storage ring length with 868 equally spaced bunches of 0.23 mA per bunch, having at their extremes two bunches of 1 mA. The remaining 1/8 of the storage ring is filled in its centre by a single bunch of 2 mA. This results in a typical RMS duration of the pulses of RMS_{t} ≃ 20 ps. The samples dimensions are reported in Table 1 together with the experimental beam parameters: timeaveraged , photon energy E_{0} and beam size.
An overview of a typical Bi_{2}Sr_{2}CaCu_{2}O_{8+δ} chip is reported in Fig. 1(a). The irradiated region is located between the voltage electrodes (framed in black) and Xrays perpendicularly impinge on the ac plane. As determined during scanning electron microscope (SEM) analysis, this region consists of a section of the Bi_{2}Sr_{2}CaCu_{2}O_{8+δ} crystal mounted on top of a sapphire substrate (0.5 mm thick) and separated from it by a thin layer of air ∼1 µm in thickness.
2.2. Simulation strategy and geometry
The selection of the simulation method is based on the analysis of the deexcitation processes. The electron τ_{el} during which nonthermal phenomena can be observed is less than 1 ps for the systems considered (Lindgren et al., 1999; Truccato et al., 2006; Medvedev, 2015; Kaiser et al., 2000), while the timescale of the heating process is defined by the electronlattice thermalization τ_{elph}, which is typically a few picoseconds long (Lindgren et al., 1999; Truccato et al., 2006; Gadermaier et al., 2010; Medvedev, 2015; Kaiser et al., 2000). These timescales must be compared with the duration of the Xray pulse length t_{p} ≃ 2.35 RMS_{t}. Synchrotrons have t_{p} values in the order of tens of picoseconds; therefore, in the first instants of the pulse the sample rapidly reaches the electronlattice equilibrium and then starts dissipating the absorbed energy by following the thermal route. Accordingly, the spatial dispersion of the absorbed energy within the target, triggered by the emission of photoelectrons and Auger electrons, starts well before the thermal deexcitation process. We have already explored the concept of energy dispersion by means of radiationtransport simulations of the Xray nanopatterning process, finding that the photogenerated electrons exist well beyond the exposed volume and can spread around up to 200 nm from the beam centre (Torsello et al., 2018). Based on these considerations, the fraction of absorbed energy and its spatial distribution, owing to the quick electronrelaxation mechanisms, must be considered to calculate, in an accurate way, the slower and consequent temperature increase. Therefore, the new simulation strategy proposed relies on the coupling of the radiation and heattransport physics.
In the present work, when the distribution of the absorbed Xray energy is considered, it was calculated through the MC using the MCNP6 code (Goorley et al., 2012) (https://mcnp.lanl.gov/). Each geometry was limited to the region involved in the transport of radiation and particles. Each physical event was tracked by means of a mesh of elementary cubic volumes called voxels. In order to reach a good compromise between space resolution and computation time, the size of voxels edges was fixed to 50 nm.
The heat transport was calculated with the FEM with the software COMSOL Multiphysics (Version 4.3b; COMSOL AB, Stockholm, Sweden). The heating source was implemented in agreement with the distribution of energy calculated by means of the radiationtransport simulation. The origin of the frame of reference was set to coincide with the point where the beam strikes the samples. Concerning the geometries, to model the heating of the glass bead, two concentric spheres were used. The inner one represents the glass domain of the bead and the outer one represents the surrounding air domain. The volume of the air domain was fixed to be 20 times the bead radius to safely exclude any influence of the boundary conditions on the results. For the prediction cases, the geometry implemented is more complex and reports only the regions of interest [see Fig. 1(b)]. It consists of a 140 µmlong Bi_{2}Sr_{2}CaCu_{2}O_{8+δ} crystal, representing the crystal portion between the voltage electrodes, surrounded by air (100 µm thick) and by the sapphire substrate (100 µm thick). The silver electrodes were not reported in the models because their distance from the irradiated region (∼70 µm) is large enough not to affect the quality of the results. The thermal contact between crystal and substrate, which is guaranteed in the proximity of the silver electrodes, was represented by means of two sapphire pads located at the extremities of the crystal length [see Fig. 1(b)].
2.2.1. MCNP6 code
In MCNP6 code, both the photoatomic interactions (photoelectric effect, Compton scattering, coherent and incoherent scattering) and the electron interactions (electronenergy loss, electron angular deflection, Bremsstrahlung emission, Auger electrons, knockon of electrons) are considered. However, some approximations are made. The of the material is not taken into account, so that the materials are always considered as amorphous. Moreover, a lowenergy threshold of 1 keV is applied for particle tracking, which neglects the fact that the lowenergy particles can move further in space and therefore slightly overconfines energy deposition by the beam in the material.
The photon source was defined specifying the energy of the impinging photon E_{0} and the experimental Gaussian profile of the Xray beam (see Table 1). The latter was approximated by a 2D Gaussian function as follows:
with W_{x} and W_{z} defined as 0.85 FWHM_{i}, and FWHM_{i} corresponding to the beam sizes of Table 1.
The deposited energy density for each voxel is expressed per impinging photon E(x, y, z). In order to reach a valid statistical result, the final spatial distribution was calculated with a total number of 10^{8} trial photons.
Concerning the supporting information), and the micrometric beam size is very small. Therefore, for this case, the contribution of photoelectrons and is expected to play a minor role in the heating effects. Secondly, because both the beam size and the attenuation length are much smaller than the sample diameter, an important fraction of the volume is not affected by the relaxation phenomena considered in the MCNP6 code. On the contrary, in the latter cases, along the direction perpendicular to the beam, the spatial length of the energyredistribution processes is comparable with the nanometric dimension of the spot size. This is represented in Fig. 2, in which the energydensity distribution for the prediction cases WBVB08, WBAP13 and WBAP14 is reported as a function of the beam penetration direction y and of the two perpendicular directions z [Figs. 2(a), 2(b) and 2(c)] and x [Figs. 2(d), 2(e) and 2(f)], respectively. The effects of energy dissipation are evident in the distribution widths, which are wider than the corresponding FWHM of the beam profile. Nonetheless, the energydensity distributions are still correlated with the beam shape. This is highlighted by the fact that the lateral spread of the energy densities rescales coherently with the Xray beam size adopted in each experiment. Finally, on the beam direction, a fraction of the incident energy is transmitted. From these calculations, it is visible how, in these cases, an accurate definition of the energy distribution becomes important.
conditions, an important difference exists between the validation and the prediction cases. In the former, the ratio between the radial redistribution of the deposited energy (a few hundreds of nm, Fig. S1 in theIn Table 2 the peak values of the energydensity distributions (E_{max}) are summarized for each sample.

2.2.2. COMSOL model
Since t_{p} is much greater than both τ_{e} and τ_{elph}, the classical Fourier heat diffusion equation should provide a good description of the temperature evolution. The Fourier heat equation is then used:
where T represents the space and timedependent temperature field, ρ, C_{p} and k are the density, the at constant pressure and the of the materials (see values listed in Table 3), respectively, and Q represents the heat source. The thermal coefficients of the materials were created to be constant in the temperature range of interest. Boundary conditions were assumed by imposing a constant temperature, equal to the ambient temperature, on all the model external boundaries.

In the explicit case, the heat source Q of the FEM model was set equal to the power density absorbed by the crystal P_{absex},
with Φ(t) representing the instantaneous and E(x, y, z) representing the deposited energydensity distribution as calculated by MC simulations. Moreover, in order to highlight the impact of the MC simulation in the Xray heating process, the same model was also solved ignoring the physics of the beam–crystal interaction. In this case, the space distribution of the heating power Q was implicitly considered equal to the space distribution of the Xray absorption, as appears from the experimental conditions. Therefore, along the beam direction (i.e. the y axis) the power absorbed was represented in agreement with the Lambert–Beer law, and along the two directions normal to the beam direction it was defined by means of the experimental intensity profile of the incident beam [equation (1)]. The final form of the implicit definition of the heat source is equal to
where λ represents the attenuation length of the material [51 µm for glass and ∼17 µm for Bi_{2}Sr_{2}CaCu_{2}O_{8+δ} (Henke et al., 1993)]. For simplicity, in the following we will refer to the two approaches, resulting from the use of the two heat sources Q_{ex} and Q_{im}, as the explicit and implicit methods, respectively.
By substituting the instantaneous Φ(t) with the experimental timeaveraged (see Table 1), the total timeaveraged power density absorbed by the samples was obtained. In Table 2, with and the corresponding maximum values are reported for the implicit and explicit methods, respectively.
A few hundreds of thousands of elements were created in the mesh, and a computer with 88 Gb of RAM and two processors with a clock of 2.4 GHz were used for the calculations. The computing time for all the simulations did not exceed 12 h.
3. Results and discussion
3.1. Validation case
We have carried out the simulation corresponding to the validation case by calculating the transient behaviour induced in the experimental conditions in the work of Snell et al. (2007) in the presence of a timeaveraged = 3.24 × 10^{12} photons s^{−1} (see Table 1). In Fig. 3, the evolution of the experimental (dots) and the calculated (lines) heating for the glass bead is reported within the first 5 s of irradiation. The data refer to the average temperature at the sample surface in the area impinged by the beam. The solid lines define the limits of the range (highlighted in grey) of all of the possible behaviours predicted with the explicit method. The upper and lower limits were found by considering the experimental ranges shown in Table 1. A good agreement between experimental and calculated behaviours is shown. Indeed, except for the first instants of irradiation, the experimental data fall inside the simulated range. It is also important to note that if other experimental uncertainties were implemented, the range of existence of the simulated behaviour (which was found considering only the variations in the sphere diameter and in the photon flux) would be even wider. Among the uncertainties, the geometry and materials of the sphere holder, the shape irregularities present in the sample geometry, and the imperfect alignment of the beam direction with the centre of the bead would certainly play a role in further increasing the grey area of Fig. 3. This means that to have an accurate prediction of the Xray heating, all the experimental parameters must be carefully controlled and the experimental uncertainties minimized.
In Fig. 2, the theoretical behaviour predicted by the implicit method is also reported with the dashed lines. As expected, the heating for the implicit method is slightly greater than the one calculated with the explicit method. However, this is coherent with the fact that the ratio between the nanometric length of the diffusive behaviour of the relaxation phenomena and the micrometric dimension of the Xray beam is very small.
3.2. Prediction cases
3.2.1. Steadystate solution
The steadystate solution of the model can be easily obtained by imposing a vanishing time derivative in equation (2) and assuming that the heat source is time independent, then . The corresponding temperature distributions obtained defining the heat source with the explicit method are shown in Fig. 4 for the samples WBVB08 [Fig. 4(a)], WBAP13 [Fig. 4(b)] and WBAP14 [Fig. 4(c)]. It can be observed that the temperature increase is significant only within ∼30 µm from the beam, which confirms a posteriori the possibility of safely neglecting in the model the presence of the Ag contacts. Furthermore, air heating can be appreciated only within a maximum distance of ∼3 µm from the crystal surface. Therefore, the confinement of the temperature increase in these regions also confirms the validity of the boundary conditions assumed for the model. The maximum temperature T_{max} is ∼314, 307 and 303 K for the samples WBVB08, WBAP13 and WBAP14, respectively.
This is in accordance with the fact that sample WBVB08 is the one with the highest maximum value . However, this seems to be not the only parameter that determines T_{max}. Indeed, although has the lowest intensity in WBAP13, the temperature rise calculated in this case is not the lowest one. This indicates that T_{max} is also strongly dependent on the energy spatial distribution of the heating source.
Figs. 5(a), 5(b) and 5(c) display the corresponding temperature profiles calculated inside the crystals along the beam direction for samples WBVB08, WBAP13 and WBAP14, respectively. A temperature decrease can be noticed near the incident surface owing to the cooling effect of air. In this regard, measuring the distance dT_{max} between the incidence surface and the position of T_{max} and normalizing it with respect to the crystal width w, a linear correlation is found with T_{max} [Fig. 5(d)]. This implies that the higher the temperature, the closer (in relative units dT_{max}/w) the T_{max} point will be to the Xray incidence surface. The temperature behaviours in the crystallength direction (x axis of the model) are reported in Fig. 5(e) for samples WBVB08, WBAP13 and WBAP14 in red, blue and green lines, respectively. The temperature profiles decay to room temperature less rapidly with the increase of T_{max}. Moreover, at most 40 µm away from the incidence point all the samples reach the ambient temperature, confirming the confinement of the heating effect.
3.2.2. Transient solution
The evolution of the temperature in the crystals under transient conditions, e.g. after turning the beam on at t = 0, can be studied by restoring the timederivative term in equation (2), while keeping the heatsource term Q equal to a timeindependent constant . The corresponding time behaviours shown in Fig. 5(f) for samples WBVB08, WBAP13 and WBAP14 calculated at the hottest point are represented with the solid red, blue and green curves, respectively. In the same figure, with the same colours, the corresponding behaviours calculated with the implicit method are represented with dashed lines. In general, three regimes can be distinguished: (i) an initial very fast increase of the temperature (typically less than a few microseconds), (ii) a slowdown of the heating rate and (iii) a final steadystate condition. Comparing the results from the implicit and the explicit models, the overestimation of the heating response obtained with the former method is clear.
3.2.3. Timedependent solution
The implementation of the model with the pulsed nature of synchrotron radiation adds a remarkable complication to the problem by concentrating the absorbed power density in time intervals of tens of picoseconds and by modulating the pulse features in frequency and intensity. In order to obtain a closer description of the experimental situation, the realtime dependence of the heatsource term, Q = Q(t), was implemented in the model. The single time pulse of the beam is Gaussian in shape and has an RMS duration RMS_{t} defined according to the relation FWHM_{t} ≅ 2.35 RMS_{t}. Therefore, in both equations (3) and (4), Φ(t) can be rewritten for each pulse as
where W_{t} = FWHM_{t}/(2ln2 )^{1/2} 0.85 FWHM_{t} 1.996 RMS_{t} represents the t_{0,i} represents the time of the ith peak pulse, and Φ_{0}(t_{0,i}) is the corresponding instantaneous Since the time pattern of the beam is periodically repeated with a period T_{rev} = 2.82 µs, Φ_{0}(t_{0,i}) can be determined for each mode of the storage ring by requesting that
The values of the parameters used to describe the time modulation of the synchrotron radiation are listed in Table 4.

In 16bunch filling mode, used for samples WBVB08 and WBAP13, all of the pulses are equal and Fig. 6(a) shows the time behaviour of the power density at the beam incidence point Q(0, t) corresponding to the first eight bunches, normalized to their peak value Q_{tmax} which corresponds to the most intense power density at the incidence point obtained at t_{0,i}, i.e. when the absolute maximum Φ_{0}(t_{0,i}) is reached. On the other hand, in the 7/8+1 mode (sample WBAP14) three different kinds of pulses are present, with intensities of 0.23, 1 and 2 mA bunch^{−1}. The instantaneous corresponding to the smallest peak can be obtained from equation (5) as . The time evolution of the normalized power density Q(0, t)/Q_{tmax} for sample WBAP14 is represented in Fig. 6(b) for two representative portions of the time pattern.
16bunch irradiation of WBVB08 and WBAP13. The comparison between the temperature time behaviours of the transient and the timedependent solutions, as calculated at the hottest point with the explicit method, is reported in Fig. 7(a) for WBVB08 and in Fig. 7(b) for WBAP13. Because of the dense time packing of the Xray energy, the temperature rapidly increases in correspondence with the pulse arrival. Indeed, the instantaneous heating power density is four orders of magnitude larger than the respective steadystate cases. In the first pulse, the temperature rises to a maximum of ∼340 and 312 K within ∼350 ps in samples WBVB08 and WBAP13, respectively. Then, it slowly decreases in the interval between two pulses, reaching a minimum temperature of ∼300 K immediately before the second pulse arrives. In Figs. 7(a) and 7(b), it is also possible to appreciate how the difference between the transient (red and blue curves) and the pulsed solutions (corresponding black curves) becomes more evident with increasing the maximum temperature achieved with the pulsed heat source. This fact clarifies how considering synchrotron experiment's measurements and simulations, that average the sample temperature over time scales of the order of 100 ns or more, can result in an underestimation of the real instantaneous sample temperature.
In Figs. 7(d) and 7(e), a blowup of the first temperature peak, as calculated with the explicit (solid line) and implicit (dashed lines) methods, is reported for samples WBVB08 and WBAP13, respectively. It is possible to observe that the temperature increases simulated with the implicit method are about one order of magnitude higher than those obtained with the explicit method. Precisely, an overestimation of ∼500 and 90 K for the samples WBVB08 and WBVB13 is obtained with the implicit method, respectively. By comparing these results with those shown in Fig. 5(f), for the corresponding transient cases, it can be observed how the temperature difference between the implicit and explicit methods becomes larger when the time modulation of the heat source is considered.
7/8+1 irradiation of WBAP14. The heating evolution calculated with the explicit method at the maximum temperature point of WBAP14 is reported in Fig. 7(c). The more complex structure of the irradiation time pattern is reflected in the temperature evolution of the sample: the two temperature peaks of ∼301.5 and 303 K (at t ≃ 0 and t ≃ 2500 ns, respectively), corresponding to the two 1 mA bunches of the storage ring, delimit the train of 868 smaller bunches, generating lower temperature peaks whose maxima lie between 300 and 301 K. The highest temperature peak at ∼304 K corresponds to the single bunch of 2 mA (at t ≃ 2600 ns). It can be noted that the temperature increase is lower than in the cases in 16bunch mode. This is consistent with the peak value of the heat source of 1.9 × 10^{17} W m^{−3} (2 mA pulse), which is considerably less than those calculated for samples WBVB08 and WBAP13 of 1.3 × 10^{18} W m^{−3} and 3.2 × 10^{17} W m^{−3}, respectively. By comparing the timedependent solution with the transient solution (green line), it can be noted that in the region of the closely spaced pulses at t ≃ 0 ns the first three peaks have a minimum temperature which is higher than the transient behaviour. This is caused by the presence of the 1 mA peak at t ≃ 0 in the power density Q(t). Indeed, after the repetition of the 868 pulses, the temperature of the transient case is stabilized and the baseline of the pulsed profile falls below the transient behaviour as expected, because the transient solution should represent a sort of time average of the real pulsed regime. A maximum difference of ∼2 K from the transient behaviour can be detected in this region around t ≃ 2440–2800 ns for the last two pulses.
In Fig. 7(f), the temperature evolution calculated from the explicit method is compared with the one derived from the implicit method for the most intense peaks observed in sample WBAP14. A maximum overestimation of the heating of ∼40 K is obtained. Therefore, the lower heating predicted with the explicit method can be considered as an improvement in the accuracy of the model. This fact not only clarifies the role of the absorbedenergy spatial distribution in the thermal phenomenon but it also allows us to effectively determine the impact of the thermal effects in the observedmaterial property changes. From this point of view, the hypothesis that very high temperature peaks could have been induced by Xrays during these experiments must be rejected. Anyway, the possibility that synchrotron pulses may induce thermal stresses in the crystals cannot be excluded. Indeed, plastic deformation and cracks formation can occur because of the strong thermal gradient (Moshe et al., 2000; Nicholson et al., 2001; HauRiege, 2012). Moreover, even in the case that thermal and mechanical phenomena are not important, the repetitive pulsed irradiation can cause thermal fatigue and fracture (de Castro et al., 2010; Ryutov, 2003). The numerical assessment of these possibilities requires more dedicated investigations.
The temperature distributions calculated at the instant of maximum temperature increase with the explicit method are reported in Fig. 8. For all of the samples, the heating effects are almost completely confined within the beam spot size.
3.2.4. Adiabatic approximation
It is important to note that without considering the photoatomic interactions and the electron interactions, it is not possible to correctly estimate the instantaneous temperature increase for the experiments. On the contrary, if the spatial extension of these interaction is given by any simulation code, it is possible to exploit the adiabatic approximation to obtain a preliminary estimation of the results of the numerical solution of the heattransfer equation. Indeed, if we consider the predictive cases corresponding to the Bi_{2}Sr_{2}CaCu_{2}O_{8+δ} Xray nanopatterning experiments, the thermal diffusivity of the material is given by . This implies that for any typical radius Δr of the cylindrical volume where energy deposition takes place around the beam axis, the time Δt needed by the heat to diffuse out of this region is of the order of Δr^{2}/α. It can be easily checked, even in the case of a very small Δr value in the order of magnitude of 100 nm, then Δt ≃ 4.8 ns, which is much longer than any Xray reproduced by our numerical simulations (20–50 ps). This means that the energy deposition is so fast that during each pulse there is no time for the heat to diffuse out of the energydeposition region, which makes the adiabatic approach meaningful. Therefore, during such a short timescale, the energydeposition volume can be considered as isolated from the rest of the sample and the corresponding temperature increase can be calculated, which is expected to correspond to the maximum temperature reached in the sample.
Therefore, the definition of the energydeposition volume is an important step. Because of the irradiation geometry, this volume can be identified with a cylindrical volume with radius Δr and height Δl. Keeping in mind that we are interested in estimating the maximum temperature to be observed in the sample, it is expected that the point where this temperature is reached has to be close to the sample surface, which means that this point lies in the sample within a depth that is much less than the Xray attenuation depth λ. Therefore, the choice Δl = λ/10 is expected both to contain the maximum temperature point and to be reasonably uniform in energy deposition and in temperature increase. Concerning the radius Δr, this can be estimated only on the basis of the MC results. By fitting to Gaussian curves (see supporting information) the profiles of the energydeposition curves plotted in Fig. 2, the FWHM_{xz} values reported in Table 5 can be obtained.

Therefore, a natural choice in this respect involves assuming Δr = FWHM_{xz}/2. Then, the maximum temperature increase ΔT can be estimated by considering all of the energy absorbed by the cylinder during a single bunch as heat adiabatically delivered to it,
where E_{bunch} is the energy corresponding to all of the photons of a single bunch, and m = ρπΔr^{2}Δl. The corresponding values are listed in Table 5 and compare reasonably well with the results from the numerical solutions.
4. Conclusions
In this work, the Fourier heat equation was solved with the finite element method to simulate the Xray heating effects induced by thirdgeneration synchrotrons in inorganic materials. Two approaches were used to define the heatsource term. In addition to the already reported approach, in which the heat source is implicitly approximated by using the experimental parameters coming from the beam profile and the material attenuation length, we proposed a second method in which the heat source is defined explicitly. This is carried out using the Monte Carlo method which explicitly considers the photon and particle interactions involved in the absorption and relaxation processes of the target. Both of the models were applied to simulate the heating measurements carried out by Snell et al. (2007) in a millimetric glass bead, finding a good agreement. However, it also emerges how a control of the experimental parameter is fundamental to accurately model the system and reduce the variation of the predicted behaviour.
In the other scenario explored, the role of heating in three examples of Xray nanopatterning experiments, performed in Bi_{2}Sr_{2}CaCu_{2}O_{8+δ} microcrystals (Mino et al., 2017; Pagliero et al., 2014; Truccato et al., 2016), was predicted. For these cases, the theoretical prediction of the heating behaviours was studied modelling the heat source in stationary and timedependent conditions, i.e. the filling pattern of the storage ring was considered when modulating in time the heat source.
In the stationary case, the maximum temperature achieved is found to be dependent not only on the power density but also on the spatial distribution of the heat source in the plane normal to the beam. Moreover, an overestimation of up to tens of degrees is obtained with the implicit method. This strong difference was attributed to the energy redistribution that takes place because of photoelectron emission and transport. This implies that the more accurate definition of the powerdensity distribution in the explicit method is fundamental when Xray nanobeams are considered, i.e. when the irradiated region is in the same order of the diffusion length of the photoelectrons. In the timedependent problem, in which the pulsed behaviour of synchrotrons is considered, the heating phenomenon was demonstrated to be roughly proportional to the instantaneous peak value of the absorbed power density. Finally, in this more complex regime, the comparison between the implicit and explicit methods highlighted the crucial impact of the powerdensity distributions. Indeed, in the implicit method, which disregards the energy distribution taking place with photoabsorption events, an overestimation of the thermal effects of the order of hundreds of degrees is predicted. Conversely, if the spatial redistribution of the absorbed energy is explicitly considered, an adiabatic approach can estimate the temperature increase reasonably well.
Ultimately, from the point of view of the Xray nanopatterning process, the calculated instantaneous temperature increase, for the three considered samples, varies from a few degrees to tens of degrees. Therefore, by considering that Bi_{2}Sr_{2}CaCu_{2}O_{8+δ} has a melting temperature of ∼860°C, an ordinary melting process can be excluded.
Supporting information
Figures S1 and S2. DOI: https://doi.org/10.1107/S1600577520010553/fv5124sup1.pdf
Acknowledgements
The authors thank the ESRF for the beam time allocated.
Funding information
MT and VB acknowledge partial support from the `Department of Excellence' (L.232/2016) grant, funded by the Italian Ministry of Education, University and Research (MIUR). This work has been partly carried out under project NANOX, jointly approved and funded by University of Turin and Compagnia di San Paolo.
References
Adriaens, A., Quinn, P., Nikitenko, S. & Dowsett, M. G. (2013). Anal. Chem. 85, 9556–9563. Web of Science CrossRef CAS PubMed Google Scholar
Archer, D. G. (1993). J. Phys. Chem. Ref. Data, 22, 1441–1453. CrossRef CAS Web of Science Google Scholar
Bahk, S. W., Rousseau, P., Planchon, T. A., Chvykov, V., Kalintchenko, G., Maksimchuk, A., Mourou, G. A. & Yanovsky, V. (2004). Opt. Lett. 29, 2837–2839. Web of Science CrossRef PubMed CAS Google Scholar
Bonino, V., Agostino, A., Prestipino, C., Hernandez, O., Fretto, M., Mino, L. & Truccato, M. (2018). CrystEngComm, 20, 6667–6676. Web of Science CrossRef CAS Google Scholar
Bonino, V., Mino, L., Agostino, A., Prestipino, C., Fretto, M. & Truccato, M. (2019). Proc. SPIE, 11035, 110350I. Google Scholar
Bras, W. & Stanley, H. (2016). J. NonCryst. Solids, 451, 153–160. Web of Science CrossRef CAS Google Scholar
Burghartz, S. & Schulz, B. (1994). J. Nucl. Mater. 212–215, 1065–1068. CrossRef CAS Web of Science Google Scholar
Castro, A. R. B. de, Vasconcellos, A. R. & Luzzi, R. (2010). Rev. Sci. Instrum. 81, 073102. Web of Science PubMed Google Scholar
Crommie, M. F. & Zettl, A. (1990). Phys. Rev. B, 41, 10978–10982. CrossRef CAS Web of Science Google Scholar
Feldman, Y., Lyahovitskaya, V., Leitus, G., Lubomirsky, I., Wachtel, E., Bushuev, V. A., Vaughan, G., Barkay, Z. & Rosenberg, Y. (2009). Appl. Phys. Lett. 95, 051919. Web of Science CrossRef Google Scholar
Gadermaier, C., Alexandrov, A. S., Kabanov, V. V., Kusar, P., Mertelj, T., Yao, X., Manzoni, C., Brida, D., Cerullo, G. & Mihailovic, D. (2010). Phys. Rev. Lett. 105, 257001. Web of Science CrossRef PubMed Google Scholar
Gnodtke, C., Saalmann, U. & Rost, J. M. (2009). Phys. Rev. A, 79, 041201. Web of Science CrossRef Google Scholar
Goorley, T., James, M., Booth, T., Brown, F., Bull, J., Cox, J., Durkee, J., Elson, J., Fensin, M., Forster, R. A., Hendricks, J., Hughes, H. G., Johns, R., Kiedrowski, B., Martz, R., Mashnik, S., McKinney, G., Pelowitz, D., Prael, R., Sweezy, J., Waters, L., Wilcox, T. & Zukaitis, T. (2012). Nucl. Technol. 180, 298–315. Web of Science CrossRef CAS Google Scholar
HauRiege, S. P. (2012). HighIntensity Xrays – Interaction with Matter: Processes in Plasmas, Clusters, Molecules and Solids. Hoboken, New Jersey: Wiley. Google Scholar
Henke, B. L. G. E. M., Gullikson, E. M. & Davis, J. C. (1993). At. Data Nucl. Data Tables, 54, 181–342. CrossRef CAS Web of Science Google Scholar
Hsu, P.C., Chen, Y.S., Hwu, Y., Je, J. H., Margaritondo, G. & Tok, E. S. (2015). J. Synchrotron Rad. 22, 1524–1527. Web of Science CrossRef IUCr Journals Google Scholar
Kaiser, A., Rethfeld, B., Vicanek, M. & Simon, G. (2000). Phys. Rev. B, 61, 11437–11450. Web of Science CrossRef CAS Google Scholar
Lide, D. R. (2003). CRC Handbook of Chemistry and Physics. CRC Press: Boca Raton. Google Scholar
Lindgren, M., Currie, M., Williams, C., Hsiang, T. Y., Fauchet, P. M., Sobolewski, R., Moffat, S. H., Hughes, R. A., Preston, J. S. & Hegmann, F. A. (1999). Appl. Phys. Lett. 74, 853–855. Web of Science CrossRef CAS Google Scholar
Liu, X., Du, D. & Mourou, G. (1997). IEEE J. Quantum Electron. 33, 1706–1716. CrossRef CAS Web of Science Google Scholar
London, R. A., Bionta, R. M., Tatchyn, R. O. & Roesler, S. (2001). Proc. SPIE, 4500, 51–62. CrossRef CAS Google Scholar
MartínezCriado, G., Villanova, J., Tucoulou, R., Salomon, D., Suuronen, J.P., Labouré, S., Guilloud, C., Valls, V., Barrett, R., Gagliardini, E., Dabin, Y., Baker, R., Bohic, S., Cohen, C. & Morse, J. (2016). J. Synchrotron Rad. 23, 344–352. Web of Science CrossRef IUCr Journals Google Scholar
Martis, V., Nikitenko, S., Sen, S., Sankar, G., van Beek, W., Filinchuk, Y., Snigireva, I. & Bras, W. (2011). Cryst. Growth Des. 11, 2858–2865. Web of Science CrossRef CAS Google Scholar
Medvedev, N. (2015). Appl. Phys. B, 118, 417–429. Web of Science CrossRef CAS Google Scholar
Mino, L., Bonino, V., Agostino, A., Prestipino, C., Borfecchia, E., Lamberti, C., Operti, L., Fretto, M., De Leo, N. & Truccato, M. (2017). Sci. Rep. 7, 9066. Web of Science CrossRef PubMed Google Scholar
Mino, L., Bonino, V., Picollo, F., Fretto, M., Agostino, A. & Truccato, M. (2019). Adv. Electron. Mater. 5. Google Scholar
Mino, L., Borfecchia, E., SeguraRuiz, J., Giannini, C., MartinezCriado, G. & Lamberti, C. (2018). Rev. Mod. Phys. 90, 65. Web of Science CrossRef Google Scholar
Mirzoev, F. K., Panchenko, V. Y. & Shelepin, L. A. (1996). UFN, 166, 3–32. CrossRef CAS Web of Science Google Scholar
Moshe, E., Eliezer, S., Henis, Z., Werdiger, M., Dekel, E., Horovitz, Y., Maman, S., Goldberg, I. B. & Eliezer, D. (2000). Appl. Phys. Lett. 76, 1555–1557. Web of Science CrossRef CAS Google Scholar
Natividad, E., Castro, M., Burriel, R. & Angurel, L. A. (2006). J. Therm. Anal. Calorim. 84, 307–316. Web of Science CrossRef CAS Google Scholar
Neutze, R., Wouts, R., van der Spoel, D., Weckert, E. & Hajdu, J. (2000). Nature, 406, 752–757. Web of Science CrossRef PubMed CAS Google Scholar
Nicholson, J., Nave, C., Fayz, K., Fell, B. & Garman, E. (2001). Nucl. Instrum. Methods Phys. Res. A, 467–468, 1380–1383. Web of Science CrossRef CAS Google Scholar
Pagliero, A., Mino, L., Borfecchia, E., Truccato, M., Agostino, A., Pascale, L., Enrico, E., Leo, N. D., Lamberti, C. & MartínezCriado, G. (2014). Nano Lett. 14, 1583–1589. Web of Science CrossRef CAS PubMed Google Scholar
Rosenthal, M., Doblas, D., Hernandez, J. J., Odarchenko, Y. I., Burghammer, M., Di Cola, E., Spitzer, D., Antipov, A. E., Aldoshin, L. S. & Ivanov, D. A. (2014). J. Synchrotron Rad. 21, 223–228. Web of Science CrossRef CAS IUCr Journals Google Scholar
Ryutov, D. D. (2003). Rev. Sci. Instrum. 74, 3722–3725. Web of Science CrossRef CAS Google Scholar
Seitz, F. & Turnbull, D. (1956). Solid State Physics, pp. 307–449. New York: Academic Press. Google Scholar
Snell, E. H., Bellamy, H. D., Rosenbaum, G. & van der Woerd, M. J. (2007). J. Synchrotron Rad. 14, 109–115. Web of Science CrossRef CAS IUCr Journals Google Scholar
Stanley, H. B., Banerjee, D., van Breemen, L., Ciston, J., Liebscher, C. H., Martis, V., Merino, D. H., Longo, A., Pattison, P., Peters, G. W. M., Portale, G., Sen, S. & Bras, W. (2014). CrystEngComm, 16, 9331–9339. Web of Science CrossRef CAS Google Scholar
Subramanian, M. A., Torardi, C. C., Calabrese, J. C., Gopalakrishnan, J., Morrissey, K. J., Askew, T. R., Flippen, R. B., Chowdhry, U. & Sleight, A. W. (1988). Science, 239, 1015–1017. CrossRef ICSD PubMed CAS Web of Science Google Scholar
Torsello, D., Mino, L., Bonino, V., Agostino, A., Operti, L., Borfecchia, E., Vittone, E., Lamberti, C. & Truccato, M. (2018). Phys. Rev. Mater. 2, 014801. Web of Science CrossRef Google Scholar
Truccato, M., Agostino, A., Borfecchia, E., Mino, L., Cara, E., Pagliero, A., Adhlakha, N., Pascale, L., Operti, L., Enrico, E., De Leo, N., Fretto, M., MartinezCriado, G. & Lamberti, C. (2016). Nano Lett. 16, 1669–1674. Web of Science CrossRef CAS PubMed Google Scholar
Truccato, M., Agostino, A., Rinaudo, G., Cagliero, S. & Panetta, M. (2006). J. Phys. Condens. Matter, 18, 8295–8312. Web of Science CrossRef CAS Google Scholar
Tu, F., Spath, A., Drost, M., Vollnhals, F., Calderon, S. K., Fink, R. H. & Marbach, H. (2017). J. Vac. Sci. Technol. B, 35, 031601. Web of Science CrossRef Google Scholar
US (1976). US Government Printing Office, Washington, DC. Google Scholar
Wallander, H. & Wallentin, J. (2017). J. Synchrotron Rad. 24, 925–933. Web of Science CrossRef IUCr Journals Google Scholar
© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.