Controlled strong excitation of silicon as a step towards processing materials at sub-nanometer precision

Interaction of a solid material with focused, intense pulses of high-energy photons or other particles (such as electrons and ions) creates a strong electronic excitation state within an ultra-short time and on ultra-small spatial scales. This offers the possibility to control the response of a material on a spatial scale less than a nanometer — crucial for the next generation of nano-devices. Here we create craters on the surface of a silicon substrate by focusing single femtosecond extreme ultraviolet pulse from the SACLA free-electron laser. We investigate the resulting surface modi ﬁ cation in the vicinity of damage thresholds, establishing a connection to microscopic theoretical approaches, and, with their help, illus-trating physical mechanisms for damage creation. The cooling during ablation by means of rapid electron and energy transport can suppress undesired hydrodynamical motions, allowing the silicon material to be directly processed with a precision reaching the observable limitation of an atomic force microscope. Implementation of XTANT code . The hybrid code XTANT (X-ray-induced thermal and non-thermal transitions) consists of a few modules simulating the effects of X-ray FEL irradiation of a target. X-ray photoabsorption, electron cascading for electrons with kinetic energies above a certain threshold, and Auger-decays of core-holes are treated within a Monte-Carlo scheme. We used ~30,000 iterations in the MC module to gather reliable statistics. Low-energy electrons kinetics is simulated within the rate equations, assuming Fermi-Dirac electron distribution and Boltzmann collision integrals for nonadiabatic electron-ion (electron-phonon) coupling. The matrix elements entering the collision integral are obtained within a dynamical-coupling scheme developed in ref. 24 , which delivers transient coupling factor depending on electron and atomic temperatures, as well as on the transient electronic band structure. A transient potential energy surface (PES) and the electronic band structure are modeled within the transferable tight binding method 50 . Atomic dynamics is simulated with a molecular dynamics (MD) method on the evolving PES. Energy transfer from electrons via electron-ion coupling is implemented via velocity scaling at each time step 24 . The Parrinello-Rahman barostat allows XTANT to model material expansion and eventual ablation while preserving periodic boundary conditions (NPH ensemble) 28 . MD time step was set to 0.01 fs to ensure convergence of the Boltzmann collision integrals 24 . All material parameters of silicon used in the simulations can be found in refs. 24,29 . Following experimental conditions, the FEL pulse duration was set to 70 fs FWHM, and photon energies of 92 and 120 eV were used.

W ell-shaped nanostructures are fundamental for nanoscience and technology. A control over nanometer-sized details of material structures is required to enhance electromagnetic field localization of advanced nano-devices [1][2][3] . To date, nano-processing of matter has relied on techniques that utilize low-intensity radiation sources, such as lithography with highenergy beams of ions or electrons, extreme ultraviolet (EUV) photons or focused ion beam etching 4 . However, they suffer from multiple and time-consuming processing steps, and only a limited number of materials for which etchants and resists are available. Laser ablation is a promising alternative in high-intensity regime, offering rapid processing with a feature resolution reaching the submicrometer by a direct removal of material 5,6 .
The main factors limiting the quality of these techniques are: (i) the spot size of the radiation beam, (ii) fast energy transport out from the primary interaction region. The spot size of the beam, which is constrained, e.g., by the diffraction limit at the corresponding wavelength, limits feature resolution of the processing. Energy transport out of the beam focus via secondary electrons leads to undesirable damages, decreasing precision of the techniques (see Fig. 1). It also lowers the 'effective' dose absorbed by the material in the focus of energetic photon (or particle) beam.
The recent advent of fourth generation light sources, X-ray free-electron lasers (XFELs) emitting intense light with wavelength down-to angstrom and a femtosecond pulse duration, offers new possibilities to investigate, and push, the limits of ultrahigh precision material processing. Despite the expense and limited availability of XFELs at present, promising dedicated studies are on-going to improve XFEL performance for practical applications e.g., EUV lithography in high volume manufacture of semiconductors 7,8 . In contrast to conventional light sources including incoherent X-ray sources and optical lasers, XFELs can provide access to a strong electronic excitation state of solid materials within an unprecedented ultra-short time on ultra-small spatial scales. Recently, focusing of XFEL beams to a few tens of nm has been already achieved 9,10 and dedicated efforts to break sub-10 nm size [11][12][13] are on-going.
Although the dynamic response of solid to a nano-excitation has a great application potential for material processing, it is far from being fully understood due to lack of experimental tools (e.g. light sources, detection techniques) precise enough to operate at the relevant temporal and spatial scales. In this study, we benchmark a macroscopic experimental result against microscopic theoretical approaches, revealing features of the nano-excitation in silicon material. We outline the possibility to reach a sub-nanometer precision of the material processing by using nano-spots of ultrashort X-ray pulse.

Results
Contributing physical processes. Figure 2 shows typical timescales for physical phenomena associated with material processing using a femtosecond X-ray pulse. X-ray induced modification of a solid material starts with photoabsorption. In contrast to a visible-light irradiation, such a high-energy photon excites electrons from the valence band (for insulator and semiconductor) or, predominantly from deep atomic shells to highenergy states of the conduction band 14,15 . The atom with a core vacancy then readjusts via fluorescence emission or Auger process. The Auger process is the dominant relaxation channel for light elements 16 . This decay leads to an excitation of an electron from a higher shell to the conduction band at a few-femtosecond scales. The ejected photo-and Auger electrons scatter further via inelastic channels (electron-electron interaction), or elastic channels (electron-lattice interaction). Impact ionization cascading is the predominant relaxation process for high-energy electrons, typically occurring on a femtosecond timescale 17,18 . It finishes when the electron loses its energy below an impact ionization threshold. In contrast, relaxation of low-energy electrons is dominated by the electron-phonon scattering in which the electron significantly loses its energy only at longer (typically picosecond) timescales 19,20 .
Apart from that, in covalently bonded materials, transiently excited states of the electronic subsystem can lead to a strong modification of the inter-atomic forces resulting into non-thermal phase transition, e.g., melting 21 . In finite systems, high fluences could also lead to a Coulomb explosion due to charge imbalance after electron emission.
Non-thermal melting in semiconductors is a result of the promotion of a large fraction (above a few percent) of the valence electrons to the conduction band 22,23 . Such electronic excitation modifies an atomic potential energy surface (PES), leading to destabilization of the atomic lattice without an increase of its kinetic energy. Thus, atoms immediately begin to move to possibly new equilibrium positions in the new phase-on timescales much shorter than the several picoseconds required to convert the electronic energy into thermal motions via electron-phonon coupling. Only later, electron-phonon coupling brings the system into thermodynamic equilibrium 19,24 . On the picosecond timescale, the shock-wave emission and thermal diffusion are dominant effects in the formation of a crater structure. Final material relaxation completes on a timescale from nano-to micro-seconds. The morphological changes of the irradiated surface can be used to identify some of the processes and to determine their threshold fluences 25 , including nonlinear effects in soft X-ray ablation 26 .
The recently developed hybrid code XTANT (X-ray-induced Thermal And Non-thermal Transitions) models the above mentioned processes on timescale up to some few tens picosecond 27 . The model accounts for: (a) X-ray photoabsorption, creation of photo-electrons and secondary electron cascading as well as for Auger-decays of core-holes, within a nonequilibrium framework of a Monte-Carlo scheme; (b) low-energy electrons kinetics traced with the rate equations assuming Fermi-Dirac distribution and self-consistent Boltzmann collision integrals for nonadiabatic electron-ion (electron-phonon) coupling 24 ; (c) transient PES and electronic band structure modeled within the transferable tight binding method; (d) atomic dynamics is simulated with a molecular dynamics (MD) method on the evolving PES, including electron-lattice energy exchange via electron-ion coupling at each time step 24 . The Parrinello-Rahman MD scheme allows XTANT to model material expansion and eventual ablation/spallation, while preserving periodic boundary conditions that exclude particle and heat flows 28 . All parameters of silicon used in the simulations can be found in refs. 24,29 .
Experiment and analysis. Figure 3 presents a schematic view of our experimental setup (see details in the "Methods" section). We focused single X-ray pulses from the soft X-ray free-electron laser (SXFEL) of the SPring-8 Angstrom Compact Free Electron Laser (SACLA) 30 onto the surface of single crystalline silicon substrates (p-doped, orientation (100), specific resistivity 0.02 ∼ 30 Ω cm) by using a Kirkpatrick-Baez (K-B) mirror. The target surface was set normal to the incident laser beam. The laser spot measured by knife-edge scanning, has a Gaussian-like energy profile with typically horizontal and vertical diameters of 8.5 and 10.5 μm at full width at half maximum (FWHM), respectively. The pulse duration of the SXFEL was measured to be about 70 fs FWHM by a correlation monitor 31 . To regulate the laser fluence, we used thin-film filters of zirconium and silicon to attenuate the pulse energy. The pulse energy at the target surface was evaluated by using a calibrated X-ray charge-coupled device (CCD) camera. During the experiment, the pulse energy was monitored on a shot-to-shot basis with a relative uncertainty of 10% by means of a gas intensity monitor (GIM). We used the photon energy of 92 and 120 eV around the L-edge of silicon. The produced silicon craters were investigated using an imaging platform including differential interference contrast microscope (DIC), laser scanning microscope (LSM) and atomic force microscope (AFM). The analyzed results were then compared to those of the XTANT simulations (see details in "Methods" section).
We summarize observed features of the modified structures in Fig. 4a-i. Figure 4a-d show the depth of a crater as a function of the SXFEL fluence. The data are fitted with the well-known equation 5 : where the coefficient α eff fitted to the data is the so-called effective attenuation length which describes the effective energy deposition depth, L is the crater depth, F L is the FEL fluence, and F th is the damage threshold fluence.
We converted the FEL fluence into the absorbed dose via the following relation: where, D abs (x) is the absorbed dose (in eV/atom) at the depth x, F L is the incident laser fluence at the material surface, α att is the attenuation length of the considered X-ray photons, R is the reflectivity of the sample at the wavelength of the incoming pulse, and n at is the atomic number density. We assumed here normal incidence of X-rays with respect to the material surface. The experimental crater depth as a function of an incident fluence for the 92-eV irradiation is well-fitted by a single line   Table 1. One can see that the effective attenuation lengths determined for silicon are significantly larger than the photon attenuation lengths. This indicates a significant contribution of the particles and energy transport in the sample, Fig. 4 Measured crater depth in silicon as a function of the laser fluence for 92-eV and 120-eV irradiation (a-d). The symbols in a-d depict distances from the initial surface to the deepest point of a crater. The error of the laser fluence is ±10% due to the uncertainty of the energy monitors, whereas the error of the depth is ±2 nm due to the uncertainty of the microscopy. c, d are the zoom-ins in the dashed-line marked region around the fluence thresholds for the hole formation in (a, b). The colored numerals in (c) and (d) show the laser fluences in mJ/cm 2 for the experimental points discussed in the below section. e-g are atomic force microscopy (AFM) images of three typical damage structures of silicon under the 92-eV irradiation with fluences of 220, 294, and 460 mJ/cm 2 (corresponding to absorption doses of 0.47, 0.62, and 0.97 eV/atom, respectively). h, i are a contrast image in a single exposure in a laser scanning microscope (LSM), and an AFM image of a crater produced at the 120-eV laser fluence of 148 mJ/cm 2 (4.4 eV/atom), respectively. The inserted circles in (h) and (i) indicate the same hillock structure. For 120-eV irradiation, there are two main types of modifications observed. The first type was found by the LSM, but could not be seen by the AFM (fluence regime below 106 mJ/cm 2 ). The other type comes from the region where the modifications could be observed by both AFM and LSM (fluence regime above 106 mJ/cm 2 ). The inserted scale bars in (h, i) indicate the length of 2 µm. which extends the energy deposition beyond the attenuation length of X-ray photons via the physical mechanisms discussed in Fig. 2. Figure 4c, d are the zoom-ins in the region around the fluence thresholds for the hole formation for 92-and 120-eV irradiation. In such an intermediate fluence regime, the hydrodynamic motion at the material surface or the bottom of the ablation crater can lead to (i) expansion which shallows the crater, or (ii) hole formation which additionally deepens the crater. These effects seem to contribute to the fluctuation of the measured depth of crater, as shown in Fig. 4c, d.
Note that the effective attenuation length of the 92-eV irradiation (see Table 1) is similar to that of 120 eV at high laser fluences (>2 J/cm 2 ). This indicates that, despite the fact that the absorbed dose distribution along the depth is different for the two photon energies, the formation of a crater in the high-fluence regime above the ablation threshold seems to be triggered by the same physical mechanisms. At such high laser fluence, a plasma formation can take place, triggering the fast removal of material surface. Apart from that, transport of a large energy flux into the deeper layer of the material via processes such as, e.g., plasma etching, thermo-hydrodynamic can additionally contribute to the damage. These mechanisms deepen the crater additionally to the initial damage, contributing to the change in slope of the curve in Fig. 4b at fluence ca. 2 J/cm 2 . However, a detailed investigation of such a high fluence regime needs dedicated experiment and simulation techniques which are beyond the scope of the present paper.
Mechanisms of hole formation around the threshold. In the case of 92 eV, there are three typical surface modifications produced, as shown in Fig. 4e-g. At the X-ray laser fluence around 150-250 mJ/cm 2 (corresponding to the surface absorbed dose of 0.32-0.53 eV/atom, only an increase of the roughness of the silicon surface with formation of nano-cone structures could be observed (see Fig. 4e). Although the averaged dose is lower than the melting thresholds (shown below), a local phase transition could be triggered due to a spatially inhomogeneous energy deposition. The inhomogeneity of material, intensity profile of the laser, and involved physical processes are factors leading to such deposition.
At the laser fluence around 300 mJ/cm 2 (~0.64 eV/atom), dome-shaped structures surrounded by the nano-cones on Si surface are formed (Fig. 4f), indicating material expansion. The absorbed dose of 0.64 eV/atom coincides well with the predicted threshold of the thermal melting into low-density liquid (LDL) phase estimated with the XTANT code 27 . Note that the attenuation length of 92-eV photon in silicon is relatively long (~588 nm). The melting to low-density liquid phase of silicon takes place within material layers down to the attenuation length. It increases the pressure (as the simulation result shows in Fig. 5a), leading to an expansion of the irradiated material.
The threshold for forming a crater, determined by fitting the experimental result with Eq. (1), is 418 mJ/cm 2 (~0.89 eV/atom). This threshold coincides well with the XTANT-predicted threshold for non-thermal melting into high-density liquid (HDL) phase 27 . During this process, an increase of pressure (as the simulation result in Fig. 5b shows), includes a rapid switching from expansion to shrinkage at sub-ps timescale, followed by a set of 'relaxation' oscillations in the compression region. The nonthermally induced increase of density reduces the volume of material, thereby forming a dimple on the surface, without material removal. The attenuation length of silicon was measured by using synchrotron radiation from the beamline 11D of the Photon Factory at the High Energy Accelerator Research Organization (KEK) We note a co-existence of a dimple with a multi-rim structure at the edges in Fig. 4g. It is widely accepted in the literature that the formation of a rim structure is caused by the heat and pressure waves in the thermodynamics regime 5,32,33 . The doublerim observed in Fig. 4g can be qualitatively understood as an interplay between the two phase transitions: the HDL state formation at the surface in the middle of the laser spot, and the LDL state in the deeper layer and at the periphery. At first, the interaction region is rapidly heated into an LDL phase at 100-fs timescale as shown in Fig. 5a, b. After that, the transition from LDL to HDL phase leads to a shrinkage at the center of the formed liquid, which relaxes through a hydrodynamical motion. The melting surface becomes then 'frozen' via rapid heat diffusion in the thermodynamical regime, and forms a multi-rim structure with a dimple. Despite a relatively low fluence, the local stress enhancement induced by the reflection of the traveling pressure at the phase boundary, including that of HDL-LDL and the liquid-solid, may lead to a removal of the material through a process known as spallation on a timescale of few tens of ps 34,35 .
In the case of 120 eV, the threshold of hole formation, estimated with Eq. (1), is of 106 mJ/cm 2 in fluence, or 3.15 eV/ atom in the surface absorbed dose. This dose is above the threshold of 0.89 eV/atom observed for the 92-eV laser pulse. A presence of new sub-µm hillock structures far from the crater, as shown in the AFM photographs of Fig. 4i, implies that ablation had occurred. In the region below the threshold of 3.15 eV/atom, there is only a contrast enhancement which could be observed by the LSM, but was not detectable with the AFM (see the part surrounding the crater in Fig. 4h, i). Such contrast enhancement is attributed to the refractive index changing at the Si surface. The corresponding realignment of Si atoms has occurred without any observable thermal cracks, or surface roughness changing, or dome feature as in the case of 92 eV (see Fig. 4d). Analogous phenomena have been observed as a result of irradiation with optical laser pulses on femtosecond timescale 25,36 as well as with a focused ion beam 37,38 .
Energy transport effect. The refractive index changing at the Si surface induced by the 120 eV radiation indicates that the Si atoms could not gain enough kinetic energy to induce disintegration, or even significant volume expansion. It is attributed to a rapid diffusion of the excited electrons and their energy, which leads to efficient cooling of the system (note that the photon energy of 120 eV is above the L-edge of Si, which makes the attenuation length significantly shorter than in the case of 92 eV). This, in turn, can increase the dose threshold for the damage, and also provide a kinetic pathway for fast reorganization of the material. The XTANT code without considering the transport effects predicts a removal of large fragments of silicon via ablation or fragmentation, which should take place within a picosecond for the doses above~2.6 eV/atom 27 , in reasonable agreement with the observed threshold; slightly larger experimental value indicates that some part of the energy is indeed brought out of the system due to energy transport effects.
In Fig. 6, we compare the experimental thresholds of hole formation in silicon with thresholds for ablation, thermal melting into the low-density liquid, and non-thermal melting into the high density liquid estimated with the XTANT simulation. For femtosecond X-ray pulses with the photon energy of 39, 92, 120 eV, and 10 keV, a crater was formed above the thresholds of 380 (ref. 39 ), 416, 106, and 78 × 103 mJ/cm 2 (ref. 40 ) (corresponding to the absorbed doses of 1.84, 0.89, 3.15, and 0.71 eV/atom), respectively. This result shows an increase of the dose threshold with increasing absorption coefficient. Attenuation of X-ray radiation with a small absorption coefficient in the target results in a relatively uniform absorbed dose profile along the penetration depth, and consequently generates a small gradient in the excited-electron distribution. Electron transport in a system with such a small gradient is slow, preserving the energy within the irradiated region for a long time. Therefore, the modification of silicon induced by a femtosecond pulse of 92 eV and 10 keV photons takes place under the conditions close to those assumed in the simulation (no heat sinks out of the simulation box), and the damage thresholds are close to the calculated thermal-melting threshold. In contrast, the energy flow out of the affected region is relatively large in the case of 39 and 120 eV photons (due to short photon absorption length and thus large gradients in the deposited energy). Consequently, the surface modification requires higher deposited doses to compensate for this energy sink. Therefore, the effect of specific laser parameters: wavelength, pulse duration and focusing spot size on the production of damage or ablation, should be accounted for, when estimating the transport of particles and their energy within the interaction region.

Discussion
After soft X-ray excitation, the emitted high-energy photo-and Auger-electrons can travel extremely fast, but they lose their energy via the inelastic channels at a femtosecond timescale. The mean free path of the electrons at such impact energies can be of a nanometer order 41,42 . Later, energy diffusion stems from the transport of low-energy electrons scattering via the elastic channels (electron-phonon coupling). High precision processing by femtosecond X-ray pulses, therefore, requires to minimize the damage zone induced by the low-energy electrons scattering during ablation (up to picoseconds timescale). The energy sink effect would suppress undesirable phase transitions and allow to keep a sharp boundary between the interaction region and the remaining part of the material. This proof of principle is shown below. Figure 7a-d show the cross section of craters induced by the 120-eV radiation. An evolution of multi-rim structures as the laser fluence increased is marked by colored-arrows. Due to the cooling effect by rapid electron and energy diffusion during the ablation, generation of the outermost rims, as marked by the green arrows, was not observed by the AFM with a typical resolution of~2 nm. The rims, as shown by the purple and blue arrows in Fig. 7b-d, are only formed on the inside of the crater. This, as discussed above, is an indication of interplay of different phase transitions (thermal one into LDL state and non-thermal one into HDL state). The presence of "rim and hole" features on the inside of the primary crater, as shown in Fig. 7b-d indicates a possibility that a transient dome feature might be formed; however, there is no permanent dome feature observed in the 120 eV data. The transient dome may be destroyed during the solidification or due to breaking the vacuum within the interaction chamber. By optimizing the laser fluence, a rimless hole structure as shown in Fig. 7a was machined with a precision reaching the observable limitation. Here, the energy sink effect helped to avoid unwanted hydrodynamical motions of the material at the edge and the bottom of the crater.
It is expected that an even greater energy sink effect can be achieved by reducing the photo-excited volume, after optimizing the incident photon energy and minimizing beam spot size. Note that the focal spots of micrometer order in our experiment were two orders of magnitude greater than the attenuation length of 120 eV photon. In contrast, the irradiation with a laser spot as small as the attenuation length offers the gradient of deposited energy increasing as a function of the beam radius, by two orders of magnitude, along the radial axis. Consequently, the energy sink effect in 2D system (axial and radial axis) becomes much greater, allowing to significantly reduce the spatial scale of the damaged zone. By scaling down the irradiated laser spot from micro-to nano-size using advanced X-ray optics 11,12 , we expect that the direct machining with a sub-nanometer precision can be achieved.
Since the laser fluences needed to create the modifications are relatively low (see Fig. 6), large areas of material could be machined simultaneously by using the X-ray interferometer technique 43 . One could split a single FEL to several beams for a high throughput processing. Compact coherent X-ray sources, including plasma-based X-ray laser 44,45 and high harmonic generation 46,47 , could also be applied for the processing. Other approaches towards high precision processing could be through steering the electron transport within the processed material by an intense electromagnetic field. Pioneering works exploiting properties of intense optical lasers for electron transport control 48,49 indicate such a possibility in a near future.
In conclusion, we investigated the formation of damage and a crater in silicon substrates by focusing femtosecond EUV pulse from the SACLA SXFEL as a proof of concept for high control in material processing for the next generation of nano-devices. Benefits for the radiobiological research on radio-and particletherapy can also be anticipated. In parallel to experimental efforts, further development of simulation tools enabling long-timescale and large-spatial scale simulations of the processed materials is also necessary.

Methods
Monitor and regulation of the laser pulse energy. The experiment was performed by using the SACLA SXFEL operating in the self-amplified spontaneous emission mode at a repetition rate of 60 Hz. A carbon-coated multilayer K-B mirror was used to focus the SXFEL pulse onto the target. It located~85 m far from the light source point. The vertical and horizontal focal length of the K-B mirror are 2.00 and 2.65 m, respectively. The incident beam size was 10 mm, corresponding to the numerical apertures of 0.0025 (vertical) and 0.0018 (horizontal). The thin-film filters of zirconium and silicon used to regulate the laser energy are also useful for eliminating higher harmonics from the SXFEL pulse. They were placed~0.8 m far from the target to prevent damages or nonlinear effects induced by the SXFEL. The 16 different energy-attenuating filters can change the relative pulse energy from 0.0044 to 100% of the maximum value. The CCD camera was set at~0.65 m after the focal point to evaluate the pulse energy on the target surface. During the irradiation, the pulse energy was monitored for each shot using the gas intensity monitor installed in front of the K-B mirror. Calibration of the responsivity of the X-ray CCD camera and the attenuation of the filters had been performed in advance at the beamline 11D of the High Energy Accelerator Research Organization (KEK)-photon factory (PF). Implementation of XTANT code. The hybrid code XTANT (X-ray-induced thermal and non-thermal transitions) consists of a few modules simulating the effects of X-ray FEL irradiation of a target. X-ray photoabsorption, electron cascading for electrons with kinetic energies above a certain threshold, and Auger-decays of core-holes are treated within a Monte-Carlo scheme. We used 30,000 iterations in the MC module to gather reliable statistics. Low-energy electrons kinetics is simulated within the rate equations, assuming Fermi-Dirac electron distribution and Boltzmann collision integrals for nonadiabatic electron-ion (electron-phonon) coupling. The matrix elements entering the collision integral are obtained within a dynamical-coupling scheme developed in ref. 24 , which delivers transient coupling factor depending on electron and atomic temperatures, as well as on the transient electronic band structure. A transient potential energy surface (PES) and the electronic band structure are modeled within the transferable tight binding method 50 . Atomic dynamics is simulated with a molecular dynamics (MD) method on the evolving PES. Energy transfer from electrons via electron-ion coupling is implemented via velocity scaling at each time step 24 . The Parrinello-Rahman barostat allows XTANT to model material expansion and eventual ablation while preserving periodic boundary conditions (NPH ensemble) 28 . MD time step was set to 0.01 fs to ensure convergence of the Boltzmann collision integrals 24 . All material parameters of silicon used in the simulations can be found in refs. 24,29 . Following experimental conditions, the FEL pulse duration was set to 70 fs FWHM, and photon energies of 92 and 120 eV were used.

Data availability
The authors declare that all relevant data are included in the paper.

Code availability
The computer codes used in the current study are currently unavailable, however, reasonable request will be considered.