(IUCr) Time dependence of X-ray polarizability of a crystal induced by an intense femtosecond X-ray pulse

research papers

IUCrJ

Volume 1| Part 6| November 2014| Pages 402-417

ISSN: 2052-2525

doi:10.1107/S2052252514018156

PHYSICS | FELS

Open

access

Time dependence of X-ray polarizability of a crystal induced by an intense femtosecond X-ray pulse

A. Leonov,^a ^* D. Ksenzov,^b A. Benediktovitch,^a I. Feranchuk ^a and U. Pietsch ^b

^aDepartment of Theoretical Physics, Belarusian State University, 220030 Nezavisimosti Avenue 4, Minsk, Belarus, and ^bFestkörperphysik, Universität Siegen, 57072 Walter-Flex-Straße 3, Siegen, Germany
^*Correspondence e-mail: leonov.bsu@gmail.com

Edited by T. Ishikawa, Harima Institute, Japan (Received 4 June 2014; accepted 7 August 2014; online 23 October 2014)

The time evolution of the electron density and the resulting time dependence of Fourier components of the X-ray polarizability of a crystal irradiated by highly intense femtosecond pulses of an X-ray free-electron laser (XFEL) is investigated theoretically on the basis of rate equations for bound electrons and the Boltzmann equation for the kinetics of the unbound electron gas. The photoionization, Auger process, electron-impact ionization, electron–electron scattering and three-body recombination have been implemented in the system of rate equations. An algorithm for the numerical solution of the rate equations was simplified by incorporating analytical expressions for the cross sections of all the electron configurations in ions within the framework of the effective charge model. Using this approach, the time dependence of the inner shell populations during the time of XFEL pulse propagation through the crystal was evaluated for photon energies between 4 and 12 keV and a pulse width of 40 fs considering a flux of 10¹² photons pulse⁻¹ (focusing on a spot size of ∼1 µm). This flux corresponds to a fluence ranging between 0.8 and 2.4 mJ µm⁻². The time evolution of the X-ray polarizability caused by the change of the atomic scattering factor during the pulse propagation is numerically analyzed for the case of a silicon crystal. The time-integrated polarizability drops dramatically if the fluence of the X-ray pulse exceeds 1.6 mJ µm⁻².

Keywords: femtosecond pulse; X-ray diffraction; polarizability; electron density; rate equations.

1. Introduction

The first hard X-ray free-electron lasers (XFELs) (Emma et al., 2010 ; Ishikawa et al., 2012 ; Feldhaus et al., 2005 ; Pellegrini & Reiche, 2004 ; Chapman, 2009 ) are already in operation at SLAC (USA) and SPring-8 (Japan); other XFEL facilities are under construction, including the European XFEL at DESY (Altarelli et al., 2007 ). These facilities will provide ultra-bright femtosecond X-ray radiation with unique possibilities to study the structure of matter with angström resolution on a time scale of femtoseconds. Most of the current experiments using FEL radiation focus on single-shot exposure of molecules and clusters, assuming that structure data can be taken before sample destruction takes place (Neutze et al., 2000 ) on a time scale much larger than the FEL pulse length. Having this sample destruction in mind, FEL experiments on crystals are rare at present (Shastri et al., 2001 ; PSI, 2009 ). Specific experimental conditions for FEL experiments have to be defined in order to solve specific questions of solid state physics.

At present, crystal diffraction is used for monochromators or other optical elements. During the first experiments with XFEL sources it was discovered that the crystal response known from conventional experiments at synchrotron sources is maintained as long as the fluence, i.e. the deposited photon energy per sample area, is below a certain threshold (Hau-Reige et al., 2007 , 2010 ). Therefore for current experiments the crystal is illuminated by a wide beam and the focusing takes place after monochromatization. However, other experimental scenarios might be realised in future experiments. One is the photon–photon pump–probe experiment where the sample is excited by one FEL pulse followed by a second one after a time span much shorter than the repetition time of the FEL source. A respective time delay set-up has been proposed recently, equipped with four crystal reflections (Roseker et al., 2009 ). For this experiment it is important to know how both the pulse shape and the intensity of the delayed pulse differ from those of the first pulse if a highly intense FEL femtosecond-pulse propagates throughout the crystal (Bushnev et al., 2011 ; Shvyd'ko & Lindberg, 2012 ).

Additional motivation for analysis of the time evolution of Fourier components of polarizability of the crystal induced by femtosecond intense X-ray pulse is connected with investigations of compact XFEL sources on the basis of the femtosecond relativistic electron bunches produced by the laser-driven accelerators (Nakajima, 2008 ; Corde et al., 2013 ). These bunches can be used to generate coherent X-ray pulses on the basis of the effect of parametric X-ray beam instability (PXBI) in crystals (Baryshevsky & Feranchuk, 1984 ; Leonov et al., 2013 ; Baryshevsky et al., 2005 ). It is well known (e.g. Baryshevsky et al., 2005; Akhiezer & Berestetzkii, 1969 ; Ter-Mikaelian, 1972 ) that electromagnetic interaction of the relativistic electron bunch with the crystal is analogous to interaction between the crystal and the X-ray pulse with the same duration and intensity being proportional to the electron current density. Therefore realisation of the PXBI effect depends significantly on the evolution of crystal polarizability during the passage of the electron bunch.

Up to now the interaction of FEL pulses with a crystal has been described by many authors in terms of X-ray dynamical theory considering the time delay of the X-ray beam while propagating through the crystal (Shastri et al., 2001; Shvyd'ko & Lindberg, 2012; Malgrange & Graeff, 2003 ) but using time-independent atomic scattering factors (ASF). However, it was shown by Hau-Riege (2011 ) that such an approach remains valid only in the case of relatively small fluences. In our paper we will show that the major variation in the crystal polarizability being proportional to ASF originates from the alteration of the ASF as a function of the pulse duration and fluence.

In the femtosecond time range the atomic positions in a crystal are fixed and the main source of variation is the electronic excitation and Auger recombination of bound electrons induced by the X-ray beam. Because the time scale of these processes is in the same time range as the FEL pulse length, the population of electronic states of an atom and subsequently the atomic form factor become time-dependent. Under these conditions, conventional theories of X-ray diffraction that are based on the stationary X-ray susceptibility of the crystal (Authier, 2003 $[Authier, A. (2003). Dynamical Theory of X-ray Diffraction, International Union of Crystallography Monographs on Crystallography. Oxford University Press.]$ ) are no longer valid because of the fast evolution of the electron density in the crystal. Since the duration $[\tau_{\rm{D}}]$ of the formation of a diffraction peak, defined by the extinction length L_ext ( $[\tau_{\rm{D}}]$ $[\approx]$ L_ext/c ≃ 10 fs, where c is the speed of light), is comparable with the duration of the XFEL pulse it is necessary to take into account the dynamics of electronic redistribution within the atomic shells. These processes finally result in the time-dependence of the ASF and the integrated Bragg peak intensity that is proportional to the square of the Fourier components of the crystal X-ray polarizability.

The evolution of electron density of an object irradiated by an XFEL pulse can be described by the solution of rate equations for the atomic state populations (e.g. Son et al., 2011 ; Santra, 2009 , and references therein) or by the simulation of microscopic processes in terms of the Monte Carlo method (Hau-Riege, 2011). An alternative approach is focusing on the description of the evolution of the electron plasma that is created in the process of ionizing the atoms (e.g. Ziaja et al., 2002 ; Hau-Riege, 2013 , and references therein). Moreover, it was also shown (Gnodtke et al., 2012 ; Iwayama et al., 2009 ; Bostedt et al., 2010; Schorb et al., 2012 ) that the ionization dynamics of individual atoms changes substantially considering the influence of the electron plasma on the time-dependent evolution of the population probabilities. As a result, the population of the atomic configurations depends on the relation between pulse duration and the size of the cluster on the one hand and the energy distribution of plasma electrons on other (Schorb et al., 2012). Evidently, the latter effect becomes essential in the case of crystals where the electronic band spectrum differs substantially from the energy spectrum of electrons in isolated atoms and molecules.

The specific feature of our approach is based on the numerical solution of a self-consistent system of master equations that includes both the rate equations for the population of bound electrons and the Boltzmann kinetic equation for the distribution function of unbound (plasma) electrons generated by the ionization of the atoms during the pulse propagation in the medium. Such an approach allows one (i) to trace explicitly the evolution of all possible atomic/ionic configurations as it is vital for further estimation of the X-ray diffraction intensities (this means that if one considers an ion with total charge +1 the diffraction signal is different for the cases of inner and outer vacancies) and (ii) to take into account secondary ionization processes and the role of free-electron plasma in the problem of evolution of atomic states of the system. The latter part considers the band spectrum of unbound electrons and additional relaxation, such as the ionization of the atoms by the electrons, electron–electron collisions and three-body recombination. The numerical treatment of these additional processes makes the solution of the master equations very expensive. Therefore analytical expressions for the cross sections of all the electron configurations in the ions have been derived on the basis of the effective charge model (ECM) for single-particle atomic wavefunctions (Feranchuk et al., 2002 ; Triguk & Fernachuk, 2011). They have been implemented in the numerical algorithm of the solution of the master equations. The developed software, crystal evolution induced by X-ray (CEIX), is applicable to various atoms. Its possibilities are demonstrated for a Si crystal as an example.

The present paper deals with theoretical investigation of the electron density evolution of atoms arranged in a crystal and the estimation of the time-dependence of ASF during the propagation of an intense XFEL femtosecond-pulse through the crystal. As shown by Ziaja et al. (2012 ), the ASF decreases remarkably during the time of the pulse propagation through the sample. This means that the conventional linear theory of diffraction (Shvyd'ko & Lindberg, 2012), assuming a constant crystal susceptibility, is no longer valid.

We concentrate on the calculation of the population dynamics of the atomic electronic states considering bound and unbound electronic states and the resulting time-dependence of the ASF. The time-dependence of the Bragg peak intensities is estimated from the square of the structure factors making up the time-dependent ASFs. First of all it is important to analyze the role of electron density evolution during the initial stage of X-ray diffraction; therefore we describe Bragg peak intensities in terms of kinematical theory which is valid as long as the crystal thickness is smaller than the extinction length L < L_ext $[\approx]$ $[c/\omega|\chi_g|]$ (where $[\chi_g]$ is the Fourier component of the X-ray polarizability of the crystal and ω is the frequency of X-ray radiation) so that the dynamical effects are negligibly small. For silicon at 8 keV photon energy L_ext = 18.5 µm at the (111) reflection in Laue geometry, for instance (Stepanov, undated $[Stepanov, S. A. (undated). X-ray dynamical diffraction web server, https://sergey.gmca.aps.anl.gov .]$ ).

Considering its femtosecond time range, the FEL pulse will probe a snapshot of the atomic arrangement in the crystal affected by random displacements of the atoms due to thermal displacements. We suppose that the respective damping of the diffraction intensity can be effectively described in terms of the static Debye–Waller approach, causing a certain reduction in the Bragg peak intensity. Whereas this part is not considered in our approach for now, we effectively describe the evolution of the Bragg peak intensity by considering five different processes of electron redistribution in the atoms and their contributions to the ASF. The degree of electron redistribution depends on the pulse length and the pulse intensity, and becomes essential if the time necessary for complete ionization of the atoms is of the order of the time necessary to form the diffraction peak. We show results of numerical investigations at photon energies of 4 keV and 8 keV, i.e. close and apart from the Si K-edge, using a pulse length of 40 fs and a flux of 10¹² photons pulse⁻¹ (the fluence being 0.8 and 1.6 mJ µm⁻², respectively).

The present paper is organized as follows. §2 and §3 motivate the approximations and introduce the processes considered for the description of the evolution of the electron density during the propagation of an XFEL pulse through a crystal. The complete system of master equations that describes the ionization dynamics in the crystal and the algorithm of the numerical solution are described in §4. The numerical results for the evolution of the electron density are discussed in §5 followed by a description of the time-dependence of the diffraction intensities from a Si crystal described in §6. The influence of the electron density evolution on formation of the PXBI effect is estimated briefly in §7.

2. Qualitative analysis

In general, the problem of the propagation of an X-ray pulse through matter is based on the solution of the system of Maxwell equations for the X-ray wavefield coupled to the Schrödinger equation for the quantum states of the electron subsystem of the crystal. In contrast to the widespread approximation of linear X-ray optics that treats the electrons as classical oscillators (Authier, 2003 $[Authier, A. (2003). Dynamical Theory of X-ray Diffraction, International Union of Crystallography Monographs on Crystallography. Oxford University Press.]$ ), a quantum theory approach for the electron density response is required in order to take into account the variations of the atomic state populations during the interaction between the X-ray field and the crystal (Benediktovich et al., 2014 ).

First of all, let us estimate the effect of an intense X-ray laser field on a single atom using the parameters of the XFEL pulse introduced in the EuroXFEL technical design report (Altarelli et al., 2007).

The electric field strength in the photon pulse can be evaluated as (Landau & Lifshitz, 1989 )

$[{\cal E}\approx \left( {{4 \hbar \omega N_{\rm{ph}}} \over {\epsilon_0 \pi d^2 c T}} \right)^{1/2} \simeq 4 \times 10^{9}\,{{{\rm{V}}}\,{{\rm{m}}^{-1}}} \,\lt\,\, {\cal E}_a \approx 5 \times 10^{11}\,{{{\rm{V}}}\,{{\rm{m}}^{-1}}}, \eqno(1)]$

where $[\epsilon_0]$ is the dielectric constant, $[{\cal E}_a]$ = $[m^2 c^3 \alpha^3/ (\hbar e_0)]$ is the characteristic strength of the atomic field with α being the fine-structure constant; e₀ is the electron charge; d is the photon beam size; T is the pulse duration; N_ph is the number of photons per pulse; $[\hbar\omega]$ is the photon energy.

The effect of an alternating laser field on the non-resonant atomic states is defined by the ponderomotive energy (Popov, 2004 ),

$[U_{\rm{p}}\approx {{e_0^2 {\cal E}^2} \over {2 m \omega^2 }} \simeq 2 \times 10^{-9}\,{\rm{eV}}, \eqno(2)]$

that is essentially smaller than the average atomic ionization potential U_i. The probability of non-resonant ionization of atoms by a laser field can be calculated on the basis of Popov (2004). In the considered case the Keldysh parameter

$[\gamma=\left(U_{\rm{i}}/2U_{\rm{p}}\right)^{1/2}\,\, \gg 1, \eqno(3)]$

which means that the probability of under-barrier tunneling is extremely small.

Taking into account (1)–(3), one can conclude that the atomic wavefunctions represent a good basis set to describe the atom–field interaction in terms of perturbation theory.

Let us compare the typical structure of the energy spectrum of electron states in a crystal (Ziman, 1972 ) with the energy spectrum of a single atom (Fig. 1a). The overlap of the electron shells of the atoms in the crystal leads to the formation of the energy bands $[\tilde{E}_n({\bf{p}})]$ (n is the zone number, $[{\bf{p}}]$ is the quasi-momentum). The electron states with n $[\leq]$ n₀ (n₀ is the quantum number of the highest populated energy level for bound electrons) correspond to the ground state of the system, the widths of the allowed bands are defined by the exponentially small overlap integrals between neighboring atomic states (Ziman, 1972), so that the energy levels in every unit cell are approximately equal to $[\tilde{E}_n({\bf{p}})]$ $[\approx]$ E_n, found at an isolated atom. At the same time, the excited states with n $[\gt]$ n₀ correspond to the conduction band. For these states the overlap integral is large and the energy spectrum is described in the framework of the `free-electron approximation' (Ziman, 1972) by $[\tilde{E}_n({\bf{p}})]$ $[\approx]$ p²/2m. This behavior is opposite to the case of an atom in a molecule or a small cluster, where the energy of the unoccupied states is still sharp. Due to the formation of the band structure the effective ionization energy that defines the transition of the electrons from the discrete to the continuous spectrum becomes a little bit smaller in a crystal than in a molecular system.

Figure 1
(a) Comparison of energy spectra of the electron states in isolated atoms and atoms in a crystal; (b) schematic estimation of the role of free electrons in the ionization dynamics. Here $[f({\bf p})]$ is the distribution function of the free electrons, R_mfp is the electron mean free path, and L_ext is the extinction length.

Another important feature of the ionization dynamics in crystals is the role of the free electrons, which are described by the distribution function $[f({\bf{p}})]$ (Fig. 1b). The characteristic energy of the free electrons that appear due to the photoionization is defined by the photon energy p²/2m $[\simeq]$ $[\hbar\omega]$ ∼ 10 keV. The mean free path R_mfp of the electrons of such an energy in media is defined by the energy loss due to secondary ionization processes, and according to the NIST Electron Inelastic Mean-Free-Path database (https://www.nist.gov/srd/nist71.cfm ) it can be estimated as R_mfp ∼ 10 nm. At the same time, in a crystal with a thickness of the same order of magnitude as the extinction length L ∼ L_ext the percentage of ionized electrons that remain within the crystal can be approximately estimated as

$[\lambda\approx \left(1-{{R_{\rm{mfp}}} \over {L_{\rm{ext}}}}\right) \simeq 1, \eqno(4)]$

that is almost a unity. This means that, in a crystal, a considerable part of the free electrons contributes to the evolution of the electron density.

3. Basic assumptions and justifications

The contribution of free electrons to the redistribution of the electron density is essential and needs to take into account additional elementary processes in order to define the ionization dynamics during the interaction of the XFEL pulse with the crystal (Fig. 2). The interaction of the XFEL pulse with a single atom is described by photoionization and Auger processes (Son et al., 2011) (processes 1 and 2, respectively). In a crystal, the large number of electrons excited into the conduction band leads to electron–electron collisions, electron impact ionization of other atoms, and the reverse process of a three-body recombination (processes 3 to 5, respectively). A sixth process is the possibility of induced photorecombination (not shown in Fig. 2). This process is reverse to photoionization, and takes place if the free electrons of the continuous spectrum become excited into unoccupied atomic states under the influence of the electromagnetic field pulse. This process is substantially resonant and involves free electrons with momenta p_r $[\approx]$ $[[2m(\hbar\omega-E_n)]^{1/2}]$ . However, numerical results show (see §5 below) that due to the collisions with electrons and atoms the photoelectrons quickly fill the entire range of the continuous states (Fig. 4) and, hence, the contribution of the resonant photorecombination to the kinetic equation for the distribution function $[f({\bf p})]$ can be neglected.

Figure 2
Elementary processes that define ionization dynamics in the crystal.

In order to find the intensity of a Bragg peak formed by the XFEL pulse, one has to calculate the crystal X-ray polarizability taking into account the evolution of the electron density. Following textbooks (e.g. Landau & Lifshitz, 1982 ) one has to solve the Maxwell equations for the photon field (here the Coulomb gauge is used),

$[\eqalign{\nabla^2{\bf{A}}({\bf{r}},t)-{{1}\over{c^2}}{{\partial^2}\over{\partial t^2}}{\bf{A}}({\bf{r}},t)&= -{{4\pi}\over{c}}{{\partial{\bf{j}}({\bf{r}},t)}\over{\partial t}},\cr \nabla{\bf{A}}&=0,\varphi=0,}\eqno(5)]$

with $[{\bf{A}}]$ and φ being the vector and scalar potentials, respectively, coupled to the Schrödinger equation for the wavefunctions $[\Psi_a({\bf{r}},t)]$ = $[\Psi({\bf{r}}-{\bf{R}}_a,t)]$ of the electron subsystem of the atom in the crystal unit cell localized near the point $[{\bf{R}}_a]$ ,

$[\eqalign{i\hbar{{\partial\Psi_a({\bf{r}},t)}\over{\partial t}}&= \hat{H}\Psi_a({\bf{r}},t),\cr \hat{H}&= {{1}\over{2m}} \left[\hat{{\bf{p}}}-{{e_0}\over{c}}{\bf{A}}({\bf{r}},t)\right]^2\,\,+\,\,e_0V({\bf{r}}),}\eqno(6)]$

where $[V({\bf{r}})]$ is the part of periodic potential of the crystal within the considered unit cell.

The induced current density in the matter can be calculated as the sum over all cells,

$[\eqalignno{{\bf{j}}({\bf{r}},t) ={}& \sum_a \Big\{ {{ie_0 \hbar} \over {2 m}}\Big[\nabla \Psi_a^* ({\bf{r}},t)\Psi_a({\bf{r}},t)-\Psi_a^*({\bf{r}},t)\nabla\Psi_a({\bf{r}},t)\Big] \cr& -{{e^2}\over{m c}} {\bf{A}}({\bf{r}},t)\Psi_a^*({\bf{r}},t)\Psi_a({\bf{r}},t) \Big\}.&(7)}]$

According to the analysis mentioned above, the stationary single-electron wavefunctions $[\psi_n({\bf{r}})]$ of the electrons in the crystal can be used as a basis set for the solution of equation (6). Let us consider the evolution of the electron state with the quantum number l and expand the wavefunction as follows,

$[\eqalign{ \Psi_a({\bf{r}},t)&= C_l \psi_l({\bf{r}}-{\bf{R}}_a) + {\textstyle\sum\limits_{n\,\neq\,l}} C_n \psi_n({\bf{r}}-{\bf{R}}_a),\cr C_{l,n}&\equiv C_{l,n}({\bf{R}}_a,t) = a_{l,n}({\bf{R}}_a,t)\exp\left[-({{i}/{\hbar}})E_{l,n}t\right]_{\vphantom{\big|}}, \cr {\bf{A}}({\bf{r}},t)&= {\bf{A}}_{s}({\bf{r}},t)\exp\left[i({\bf{k}}{\bf{r}}-\omega t)\right]+c.c.}\eqno(8)]$

The quantum number n corresponds to the entire set of the single-electron quantum states including the wavefunctions of the continuous spectrum. The coefficients $[a_{l,n}({\bf{R}}_a,t)]$ and the slope functions $[{\bf{A}}_{s}({\bf{r}},t)]$ (temporal envelope of the pulse) (Ziaja et al., 2012) are varying due to the atom–field interaction rather slowly in comparison with the atomic frequencies.

In the numerical calculations below, the analytical single-electron approximation, ECM (Feranchuk et al., 2002; Triguk & Feranchuk, 2011 ), is used both for the functions $[\psi_{l,n}({\bf{r}})]$ and the energies E_l,n of the atomic stationary states. This approximation is based on the use of hydrogen-like wavefunctions with an effective charge for each orbital so that it provides an accuracy comparable with the results obtained by the Hartree–Fock approximation (LANL Atomic Physics Codes, https://aphysics2.lanl.gov ).

The conventional approach of calculating the linear response of a system (susceptibility) (Batterman & Cole, 1964 ) is based on the approximation a_l = 1 and a_n(t) being calculated by means of the perturbative solution of equation (6). In the present case, a lot of atomic transitions are excited at the same time due to the very strong field. This results in a significant depopulation of the initial state, which must be taken into account when calculating the non-linear and time-dependent response. If one neglects the transitions between different excited states during the pulse propagation (we assume these states to be located in the continuous spectrum), a compact equation for the function a_n(t) can be derived,

$[\eqalignno{ a_n({\bf{R}}_a,t)&= -i\,{{e_0}\over{2mc}} \int\limits_{-\infty}^{t}{\rm{d}}t'\,{\bf{A}}_{s}({\bf{R}}_a,t')\,a_0({\bf{R}}_a,t')\cr&\quad\times \left\langle \Psi_n({\boldrho})\left| \hat{{\bf{p}}} \exp\left(-i{\bf{k}}{\boldrho}\right) \right| \Psi_0({\boldrho})\right\rangle \exp\left[i(\omega_{n0}-\omega)t'\right]\semi \cr {\boldrho}&={\bf{r}}-{\bf{R}}_a\semi\quad n \neq l\semi\cr \dot{a_l}({\bf{R}}_a,t)&=-{{e_0^2}\over{4m^2c^2}}\,{\bf{A}}^*_{s}({\bf{R}}_a,t)&(9)\cr&\quad\times\int\limits_{-\infty}^{t}{\rm{d}}t'\exp\left[i\omega(t-t')\right]Y(t-t'){\bf{A}}_{s}({\bf{R}}_a,t'){a_l}({\bf{R}}_a,t'),}]$

with the response function [Y(t-t')] , which allows one to take into account the effects of memory and coherence in the atom-field interaction,

$[\eqalignno{Y(t-t')={}& \sum\limits_n \left\langle \Psi_0({\boldrho}) \left| \hat{{\bf{p}}} \exp\left(i{\bf{k}}{\boldrho}\right)\right| \Psi_n({\boldrho})\right\rangle \cr& \times \left\langle \Psi_n({\boldrho}') \left|\hat{{\bf{p}}} \exp\left(-i{\bf{k}}{\boldrho}'\right)\right| \Psi_0({\boldrho}')\right\rangle \cr&\times \exp\left[-i\omega_{n0}\left(t-t'\right)\right].&(10}]$

The resonant and non-resonant parts should be treated separately when solving equation (9) for a_l(t). It can be shown that in the non-resonant case ( $[\omega_{nl}]$ $[\neq]$ $[\omega]$ ) the kernel of the integral operator (10) is almost local in time because of the condition $[\omega T]$ $[\gg]$ 1. Then the decrease of population of the atomic ground state reduces to the rate equation

$[\dot{a}_l({\bf{R}}_a,t)= -I({\bf{R}}_a,t) \sigma^{({\rm{tot}})}(\omega)\,a_l({\bf{R}}_a,t),\eqno(11)]$

where $[I({\bf{R}}_a,t)]$ is the XFEL field intensity at the point $[{\bf{R}}_a]$ of the considered atom and $[\sigma^{({\rm{tot}})}(\omega)]$ is the total cross section of inelastic scattering of the radiation by the atom. This value can be found experimentally by measuring the intensity-dependent absorption coefficients $[\mu]$ = $[n_{\rm{res}}\sigma^{({\rm{tot}})}(\omega)]$ (n_res is the resonant atom density).

Another approximation is used in the resonant case when $[\omega_{ln_r}]$ $[\approx]$ $[\omega]$ for one of the transitions. Then the coupled equations define the populations of the resonant levels,

$[\eqalign{ i\dot{a}_l({\bf{R}}_a,t)&= -\Delta\omega{a_l}({\bf{R}}_a,t)+U({\bf{R}}_a,t)a_{n_r}({\bf{R}}_a,t), \cr i\dot{a}_{n_r}({\bf{R}}_a,t)&= -i({{\Gamma}/{2}})a_{n_r}({\bf{R}}_a,t)+U({\bf{R}}_a,t)a_l({\bf{R}}_a,t),}\eqno(12)]$

where $[\Delta\omega]$ = $[\omega-\omega_{ln_r}]$ , Γ is the width of the excited level and $[U({\bf{R}}_a,t)]$ is the coupling function defined as follows,

$[U({\bf{R}}_a,t)=-{{e_0}\over{2mc}} {\bf{A}}_{s}({\bf{R}}_a,t) \left\langle\Psi_{n_r}({\boldrho})\left| \hat{{\bf{p}}} \exp\left(-i{\bf{k}}{\boldrho}\right)\right| \Psi_l({\boldrho})\right\rangle.]$

Substituting equation (8) into equation (7) and summing over the periodic coordinates $[{\bf{R}}_a]$ of the crystal, one can find that the induced current includes only the Fourier component corresponding to the set of the reciprocal lattice vectors $[{\bf{h}}]$ :

$[\eqalignno{ {\bf{j}}({\bf{r}},t)&= \sum\limits_{{\bf{h}}} {\bf{j}}^{{\bf{h}}}_{s}({\bf{r}},t) \exp\left[i({\bf{k}}+{\bf{h}}){\bf{r}}-\omega t\right], \cr {\bf{j}}^{{\bf{h}}}_{s}({\bf{r}},t)&= {{e_0}\over{m\Omega}} \Big\{ a_0^*({\bf{r}},t)a_l({\bf{r}},t) \left\langle\Psi_l({\boldrho})\left|\hat{{\bf{p}}} \exp\left[i({\bf{k}}+{\bf{h}}){\boldrho}\right]\right| \Psi_l({\boldrho})\right\rangle \cr&\quad- {{e_0}\over{c}} {\bf{A}}_{s}({\bf{r}},t) \sum\limits_{l\,\leq\,l_m,\,j} \left|a_l({\bf{r}},t)\right|^2 F_{l,j}({\bf{h}}) \exp\left(i{\bf{h}}{\bf{R}}_j\right) &(13)\cr&\quad\times \exp\left[-W(h)\right] \Big\},}]$

where Ω is the unit cell volume and $[F_{l,j}({\bf{h}})]$ is the partial atomic scattering factor that corresponds to the transferred scattering vector $[{\bf{q}}]$ = $[{\bf{h}}]$ . It is calculated for the state $[\psi_l({\bf{r}})]$ and the coordinates $[{\bf{R}}_j]$ correspond to various atoms in the unit cell; exp[-W(h)] is the Debye–Waller factor (Batterman & Cole, 1964).

The sum is calculated over all atoms within the crystal unit cell and all bound electron states with quantum numbers l $[\leq]$ l_m that were occupied in the initial state of the system. The characteristic time for a change of the atom positions is defined by the value $[\tau_{\rm{ph}}]$ ∼ $[1/\omega_{\rm{D}}]$ ∼ 10⁻¹³ s ( $[\omega_{\rm{D}}]$ is the Debye frequency) (Landau & Lifshitz, 1982). In the considered case this time is larger than the pulse duration $[\tau_{\rm{ph}}]$ $[\gg]$ $[\tau_p]$ ∼ 10^-14 and the Debye–Waller factor has the same value as for the static crystal. So the total scattering factor of the crystal unit cell is defined as follows,

$[F({\bf{h}},t)=\sum_{l\,\leq\,l_m,\,j} \left|a_l({\bf{r}},t)\right|^2 F_{l,\,j}({\bf{h}})\exp\left(i{\bf{h}}{\bf{R}}_j\right)\exp\left[-W(h)\right].\eqno(14)]$

Far from resonance, i.e. far from the K- or L-absorption edges, the anomalous dispersion term in (13) can be neglected (Kissel et al., 1995 ). This means that only the last term in the induced current density (13) defines the diffraction intensity, and the time-dependent Fourier-component of the crystal X-ray polarizability is defined as follows:

$[\chi_{{\bf{h}}}(\omega,t)= -{{4\pi e_0^2}\over{m\omega^2\Omega}}\, F({\bf{h}},t).\eqno(15)]$

The main processes that determine the dynamics of the occupation probabilities and the time-dependence of the current density via equation (13) are the photoionization and the Auger effect. Here we assume that the ionized electrons are described by plane waves and do not contribute to the periodic susceptibility. However, they can strongly affect the bound electron population. During the pulse propagation the inner shells become depleted due to both photon-induced processes and electron–atom impact ionization.

4. Application of the rate equations for ionization dynamics in the crystal

In order to solve the evolution problem for the electron density in the crystal it is convenient to separate the whole system into three subsystems: the bound electrons (discrete spectrum), the free-electron gas (continuous spectrum) and the electromagnetic field.

(i) It has been shown in many papers [for example, Son et al. (2011), and citations therein] that the most efficient way to describe the dynamics of the bound electrons is obtained by studying the time-dependence of any electron configuration of the atom. Since a set of bound electrons at a given time represents a certain atomic configuration, their evolution can be described as time-dependent changes between different possible configurations. It may start from the neutral atom and may finish with a fully ionized atom. If one writes $[P_\lambda(t)]$ for the probability of the λ configuration at an arbitrary moment of time, then the initial condition for this function corresponds to the case where all atoms are in the ground (neutral) state,

$[P_\lambda(0)=\delta_{\lambda,0}.\eqno(16)]$

One should also stress the normalization condition for the whole set of atomic configuration probabilities that should be fulfilled for any arbitrary moment of time,

$[\textstyle\sum\limits_\lambda P_\lambda(t)=1.\eqno(17)]$

With this definition, the population of the atomic level Q_l(t) = $[\langle|a_l(t)|^2\rangle]$ in the scattering factor (14) averaged over all configurations is defined as follows,

$[\eqalign{ F({\bf{h}},t)&= \textstyle\sum\limits_{l\,\leq\,l_m,\,j}Q_l(t)F_{l,\,j}({\bf{h}})\exp\left(i{\bf{h}}{\bf{R}}_j\right)\exp\left[-W(h)\right],\cr Q_l(t)&=\textstyle\sum\limits_\lambda P_\lambda(t)g_{l,\lambda}|a_l(t)|^2,}\eqno(18)]$

where $[g_{l,\lambda}]$ is the degeneracy of this level in the configuration λ.

(ii) Electrons of the continuous spectrum appear due to photoionization, Auger recombination and electron-impact ionization. This subsystem can be described in terms of a classical one-particle distribution function $[f({\bf{r}},{\bf{p}},t)]$ normalized as follows,

$[\textstyle\int f({\bf{r}},{\bf{p}},t)\,{\rm{d}}{\bf{p}}\,{\rm{d}}{\bf{r}}= n_{\rm{e}}(t),\eqno(19)]$

where n_e (t) is the total number of free electrons per unit cell.

This subsystem includes all excited electrons as well because for any excitation they occupy the conduction bands following the free-electron approximation for overlapping electron shells of atoms in a crystal.

At the initial moment of time there are no free electrons, which corresponds to the following condition,

$[f({\bf{r}},{\bf{p}},0)=0.\eqno(20)]$

One should also note that although the photoionization cross section is not isotropic over the ejected electron direction (Landau & Lifshitz, 1989), the multiple electron–electron collisions lead to the loss of information about the initial velocity directions, so that the distribution function $[f({\bf{p}})]$ can be assumed to be isotropic over the momentum variable (Landau & Lifshitz, 2001 ).

(iii) The electromagnetic field is described by the wave packet

$[\eqalign{ {\bf{A}}({\bf{r}},t)&= {\bf{e}}_{s}\Phi({\bf{r}},t) \exp\left[ i({\bf{k}}{\bf{r}}-\omega t)\right],\cr I({\bf{r}},t)&= |\Phi({\bf{r}},t)|^2,}\eqno(21)]$

where $[I({\bf{r}},t)]$ is the intensity distribution function. Using the kinematical approximation of X-ray diffraction, the evolution of the electromagnetic field is not taken into account.

Let us consider the general form of the rate equations describing the atomic population dynamics (Son et al., 2011),

$[{{{\rm{d}}P_\lambda}\over{{\rm{d}}t}}= \textstyle\sum\limits_{\mu\,\neq\,\lambda} \left(W_{\mu\lambda}P_{\mu}-W_{\lambda\mu}P_{\lambda}\right),\eqno(22)]$

where $[P_{\lambda}]$ is the probability of the system occupying a configuration with index λ and $[W_{\lambda\mu}]$ is the probability of a transition between the configuration λ to μ in unit time.

Transitions between various atomic configurations during the XFEL pulse propagation are mainly caused by photoionization, Auger decay, electron-impact ionization and three-body recombination. The photoionization rate is given by

$[W^{\,({\rm{Ph}})}_{\lambda\mu}(t)= \sigma^{({\rm{Ph}})}_{\lambda\mu}\,J(t),\eqno(23)]$

where $[\sigma^{({\rm{Ph}})}_{\lambda\mu}]$ is the cross section of the photoionization process that corresponds to the transition from configuration λ to μ and J(t) is the photon flux function.

For the time-independent Auger process rate $[W^{\,({\rm{Ag}})}_{\lambda\mu}]$ we use the expressions given by Son et al. (2011) and Santra (2009) and modify them with all ionization potentials calculated in the framework of ECM (Triguk & Feranchuk, 2011).

The electron-impact ionization rate can be deduced on the basis of a collision integral calculation and has the following explicit form,

$[W^{\,({\rm{eii}})}_{\lambda\mu}= {{n_a}\over{2}} \int v^3 f(v)\,{\rm{d}}v \int{{{\rm{d}}\sigma_{\lambda\mu}^{({\rm{eii}})}(v'|v)} \over {{\rm{d}}v'}}\,{\rm{d}}v',\eqno(24)]$

where n_a is the number of atoms per unit cell and the parameter dependence of the cross section is organized in the way $[(v_{\rm{fin}}|v_{\rm{ini}})]$ .

Using the principle of detailed balance (Hau-Riege, 2011; Landau & Lifshitz, 2001), the rate of the three-body recombination process can be deduced on the basis of the electron-impact ionization rate,

$[\eqalignno{ W^{\,({\rm{tbr}})}_{\lambda\mu}=\,\,&{} \left({{2\pi\hbar}\over{m_{\rm{e}}}}\right)^3\,\, {{n^2_a}\over{2}} \int f(v)\,{\rm{d}}v &(25)\cr&\times \int v^{\prime\,3} {{{\rm{d}}\sigma_{\mu\lambda}^{\,({\rm{eii}})}(v|v')}\over{{\rm{d}}v}}\,\, f[{v'}^2-v^2-({{2}/{m_{\rm{e}}}}) E_{\mu\lambda}]^{1/2}\,{\rm{d}}v',}]$

where $[E_{\mu\lambda}]$ is the ionization potential that corresponds to the transition from configuration μ to λ.

It is important to stress that as long as the rates (24)–(25) depend on the electron density function (see below), the subsystems of free and bound electrons are coupled.

The dynamics of the free-electron gas density function is described by the Boltzmann kinetic equation and has the form (Landau & Lifshitz, 2001)

$[\eqalignno{ {{{{\rm{d}}{f}}({\bf{r}},{\bf{p}},t)}\over{{\rm{d}}{t}}}&= {{{\partial{f}}({\bf{r}},{\bf{p}},t)}\over{\partial{t}}} +{\bf{v}}{\boldnabla}_{{\bf{r}}}\,f({\bf{r}},{\bf{p}},t) +{\bf{F}}{\boldnabla}_{{\bf{p}}}\,f({\bf{r}},{\bf{p}},t)\cr& = I_{\rm{B}}[\,f({\bf{r}},{\bf{p}},t)].&(26)}]$

For simplicity and insight into the ongoing processes, let us make a number of additional assumptions. First of all, let us suppose that the system remains homogeneous in the lateral direction during the field–matter interaction due to the fact that the beam size in this direction is much larger than the size of a crystal cell. This means that all functions depend only on z (the axis parallel to the wavevector) and t; the wavefront itself depends on the variable z-ct.

Furthermore, the only vector that could cause an anisotropy in momentum space is the photon momentum, so that the anisotropy parameter

$[\xi_a\approx{{k_{\rm{ph}}}\over{p_{\rm{e}}}} \sim \left( {{\hbar\omega}\over{m_{\rm{e}}c^2}}\right)^{1/2} \sim 0.14 \ll 1,\eqno(27)]$

and due to thermalization the density function can be considered approximately isotropic over the momentum directions.

In the non-relativistic case, the net force $[{\bf{F}}]$ acting on an electron is defined by the uncompensated Coulomb field created by the other electrons of the continuous spectrum and the ionized atoms. This force becomes essential if the photon pulse has left the crystal but can be neglected during the passage of the pulse through the crystal. Moreover, in the non-relativistic case with the assumptions mentioned above, the diffusion term yields

$[{\bf{v}}{\boldnabla}_{\bf{r}}\,f({\bf{r}},{\bf{p}},t)\sim {{v}\over{c}}{{\partial{f}}\over{\partial{t}}} \ll {{\partial{f}}\over{\partial{t}}}\eqno(28)]$

and can be neglected as well.

As a result of these approximations, one can reduce the initial Boltzmann equation (26) to the form

$[{{\partial{f}(v,t)}\over{\partial{t}}}= I_{\rm{B}}[\,f].\eqno(29)]$

In the collision integral I_B, the following transitions should be taken into account: (i) electron-impact ionization of atoms (ions); (ii) three-body recombination; (iii) electron–electron elastic scattering. The corresponding collision integrals can be written as

$[\eqalignno{I^{\,({\rm{eii}})}_{\rm{B}}={}& n_a \Bigg[ {{1}\over{v^2}}\sum_{\lambda,\mu}P_{\lambda}\int{v'}^3 {{{\rm{d}}\sigma_{\lambda\mu}^{\,({\rm{eii}})}(v|v')}\over{{\rm{d}}v}} \,f(v')\,{\rm{d}}v'\cr& -{{1}\over{2}}v\,f(v) \sum_{\lambda,\mu}P_{\lambda}\int{{{\rm{d}}\sigma_{\lambda \mu}^{\,({\rm{eii}})}(v'|v)}\over{{\rm{d}}v'}}\,{\rm{d}}v'\Bigg].&(30)}]$

In order to derive the three-body recombination collision integral one can use the principle of detailed balance (Landau & Lifshitz, 2001), so that the corresponding cross section can be obtained on the basis of the electron-impact ionization cross section,

$[\eqalignno{ I^{\,({\rm{tbr}})}_{\rm{B}}={}& \left({{2\pi\hbar}\over{m_{\rm{e}}}}\right)^3 n^2_a \Bigg\{{{1} \over {2}} v \sum_{\lambda, \mu} P_{\lambda} \int {{{\rm{d}}\sigma_{\mu\lambda}^{\,({\rm{eii}})}(v'|v)}\over{{\rm{d}}v'}} \,f(v') \cr&\times f[v^2-{v'}^2-({{2}/{m_{\rm{e}}}}) E_{\mu\lambda}]^{1/2}\,{\rm{d}}v' \cr& -{{f(v)}\over{v^2}} \sum_{\lambda,\mu} P_\lambda \int {v'}^3 {{{\rm{d}}\sigma_{\mu\lambda}^{\,({\rm{eii}})}(v|v')}\over{{\rm{d}}v}}\cr&\times f[{v'}^2-v^2-({{2}/{m_{\rm{e}}}}) E_{\mu\lambda}]^{1/2}\,{\rm{d}}v' \Bigg\}. &(31)}]$

In order to describe the elastic electron–electron scattering we implement the scheme of relaxation dynamics for particle systems with Coulomb interaction as introduced by MacDonald et al. (1957 ).

It is important to stress again that, due to the dependence of the collision integrals (30)–(31) on the atomic configuration the probabilities $[P_\lambda]$ shown in equations (22) and (29) are coupled and must be solved simultaneously. However, as long as we use ECM (Triguk & Feranchuk, 2011), all cross sections introduced in the system of master equations can be calculated analytically (see Appendix A) with the necessary accuracy. The latter makes numerical simulations less expensive in time and resources.

5. Numerical results for atomic populations

In order to simulate the population dynamics we implemented the algorithm of Morgan & Penetrante (1990 ) to solve the Boltzmann equation and the system of rate equations. In order to check both the validity of the results predicted by ECM and the stability of the numerical algorithm for solving the system of master equations, we simulated the atomic dynamics of carbon gas without taking into account the crystal structure and contribution of unbound electrons. This system has been calculated by Son et al. (2011) using a full numerical treatment in terms of the Hartree–Fock–Slater (HFS) model (LANL Atomic Physics Codes, https://aphysics2.lanl.gov ) (see Appendix B for details).

The XFEL pulse used for calculations was specified to have a photon energy of 8 keV, a photon number of 10¹² per pulse, beam size of 1 µm × 1 µm (thus the fluence being 1.6 mJ µm⁻²), Gaussian shape with full duration of 40 fs (13 fs FWHM). All calculations have been performed for the example of a silicon crystal.

The energy of the Si K-line (1.8 keV) is more than four times smaller than the photon energy of 8 keV, resulting in a non-resonant photon-to-atom interaction. In order to estimate the electron density evolution for a photon energy closer to the silicon K-edge, where non-resonant effects could become non-negligible, we performed additional simulation for a 4 keV pulse with the same characteristics as defined above (the fluence being 0.8 mJ µm⁻² in this case).

Fig. 3 shows the probability of finding differently ionized ions in the silicon crystal as a function of time. It shows that the number of neutral atoms decreases during the time of interaction between the photon pulse and the crystal. At 8 keV photon energy the population probability decreases almost to zero by the end of the pulse and for the 4 keV case it decreases completely to zero already at half of the pulse length. The latter is remarkable considering the fact that the pulse energy is about 2.2 keV above the threshold of atomic K-shell ionization. Here most of the populated states are +6 and +7 at the end of the pulse. In the non-resonant case at 8 keV the interaction between the XFEL pulse and the electron subsystem of the atom is weak, so that the atoms are not so deeply ionized and the mostly populated states are ions with +1, +2 and +3 ionization charges.

Figure 3
Atomic population probabilities for Si crystal as a function of time: (a) 8 keV pulse, (b) 4 keV pulse.

Fig. 4 shows the distribution of kinetic energy of the free electrons as a function of time. At 8 keV, i.e. in the non-resonant case, one can see three vivid energy bands varying in time: the top band (at about 6 keV) describes the energy of the photoelectrons, the middle band (about 1.3 keV) corresponds to the energy of the Auger electrons, and the range close to zero energy describes the secondary electrons that appear due to the electron-impact ionization process. In contrast to this, the 4 keV result shows a broad spectrum corresponding to the photo (both spikes at about 2.0 keV and 3.8 keV) and Auger (middle spike at about 1.3 keV) electrons. Additionally the bands are broadened due to the fact that every step of ionization is accompanied by a certain decrease of the ionization potential and subsequently a reduction of the energy of every successive photoelectron. Moreover, the free electrons undergo elastic and inelastic scattering, which also results in a broadening of the energy distribution.

Figure 4
Free electron density as a function of time and energy: (a) 8 keV pulse, (b) 4 keV pulse.

Fig. 5 shows the total number of free electrons per atom in the crystal unit cell and the contribution of the different ionization channels in time. One can conclude that, in both cases, the near-resonant and the non-resonant one, the electron-impact ionization channel plays the dominant role for the creation of free electrons. The respective yield of free electrons via this process for the 8 keV case is almost seven times higher than those of photoionization and Auger processes. For the 4 keV case the relative contribution of the electron-impact ionization channel is about four times larger than that of the photoionization and Auger recombination but in absolute numbers two times larger than for 8 keV photons.

Figure 5
Contribution of different channels and total yield of free electrons per atom: (a) 8 keV pulse, (b) 4 keV pulse.

6. Evolution of the atomic scattering factor

The most relevant quantity for the formation of the diffraction peak is the average value ASF $[\bar{F}({\bf{q}},t)]$ describing the number of scattering electrons as a function of the momentum transfer q = $[{{\sin\theta}/{\lambda}}]$ , where θ is the scattering angle and λ is the photon wavelength. The statistical character of the ionization processes means that the ASF at a moment of time t is a random value which depends on the probabilities of finding a certain electron configuration of the atom $[P_{\lambda}(t)]$ . Let us define the amount of the average ASF $[\bar{F}({\bf{q}},t)]$ and its standard deviation $[\Delta F({\bf{q}},t)]$ as follows,

$[\eqalign{ \bar{F}({\bf{q}},t)&=\textstyle\sum\limits_{\lambda}F_{\lambda}({\bf{q}})\,P_{\lambda}(t),\cr \Delta F({\bf{q}},t)&= \left[\,\overline{F^{\,2}}({\bf{q}},t)-\bar{F}^{\,2}({\bf{q}},t)\right]_{\vphantom{\big|}}^{1/2},\cr \overline{F^{\,2}}({\bf{q}},t)&=\textstyle\sum\limits_{\lambda} F^{\,2}_{\lambda}({\bf{q}})\,P_{\lambda}(t),}\eqno(32)]$

where $[F_{\lambda}({\bf{q}})]$ is the stationary ASF value for the atomic configuration $[{\lambda}]$ at the momentum transfer $[{\bf{q}}]$ . Since the anomalous dispersion term is omitted we do not consider the energy range close to the exact resonance energy.

The calculation of the ASF value with probabilities $[P_{\lambda}(t)]$ related to one cell is performed by use of the ergodic hypothesis (Landau & Lifshitz, 2001) for the statistical ensemble of the atoms in the whole crystal. It is also supposed that the fluctuations of the ASF for atoms in different cells are not correlated. In this case the ASF dispersion contributes only to the X-ray diffuse scattering background and does not change the intensity of the coherent diffraction peak (Lorenz et al., 2012 ).

Because the free-electron distribution is broad in real space the value of the ASF mainly depends on the number of bound electrons in the atoms/ions. Fig. 6 shows the alteration of this number during the pulse length. It becomes evident that atoms lose about seven bound electrons in the near-resonant case, whereas in the non-resonant case the drop is less than three electrons per atom.

Figure 6
Average number of bound electrons per atom as the function of time: (a) 8 keV pulse, (b) 4 keV pulse.

Fig. 7 shows the ASF values as a function of momentum transfer, q, for three cases: the neutral free atom and the time-averaged ASF after the passage of the 4 keV and 8 keV XFEL pulses. We find a significant reduction in the ASF in the near-resonant (4 keV) case over all values of q. On the other hand, the drop in the ASF is small in the non-resonant case (8 keV) and is substantial at low q values only.

Figure 7
ASF as a function of q for the conventional case and after the passing of 8 keV and 4 keV pulses.

The properties of the ASF are studied in more detail for two different diffraction peaks: the 220 Bragg peak at q = 0.26 Å⁻¹ and the 222 Bragg peak at q = 0.31 Å⁻¹ that are affected by changes in both the valence and core shells.

Fig. 8 shows the time evolution of the ASF at the q position of the 220 and 222 Bragg reflections as a function of photon energy between 4 keV and 12 keV. Without interaction with the XFEL pulse the ASF is about 9.0 for both 220 and 222 Bragg peaks for all photon energies (these cases do not differ qualitatively but we have considered them in order to analyze stability of the algorithm). This value drops during the time of interaction of the XFEL pulse with the crystal. The amount of this drop increases with decreasing energy difference to the K-absorption edge. At 4 keV the total ASF decreases by 50% during the XFEL pulse of 40 fs. The inset of Fig. 8 shows the drop for the 4 keV pulse case with and without the contribution of the free electrons. It becomes evident that the free electrons contribute by about 20% to the time-dependent drop of the ASF. At the same time, the 8 keV pulse causes less photoionization damage, so that the ASF drop is less than 10% for the mentioned reflections.

Figure 8
Evolution of the average ASF $[\bar{F}({\bf{h}},t)]$ (a, b) and its standard deviation $[\Delta F({\bf{h}},t)]$ (c, d) as functions of the photon energy. Inset: evolution of the average ASF $[\bar{F}({\bf{h}},t)]$ and its standard deviation $[\Delta F({\bf{h}},t)]$ without (red line) and including (black line) the contribution of free electrons for the 4 keV case. (220) and (222) reflections are considered.

The critical point of X-ray diffraction with XFEL pulses is to find the photon intensity that initiates complete ionization of the atom during a time faster than that necessary for the formation of the diffraction peak, i.e. faster than the pulse time. This threshold intensity can be determined using the numerical results shown in Fig. 9. It shows the flux dependence of both reflections and their standard deviation. It becomes evident that the form factor drop is dramatic if the fluence exceeds 1.6 mJ µm⁻².

Figure 9
Dependence of the average ASF (a, b) and its standard deviation (c, d) on the fluence of the XFEL pulse at 8 keV photon energy.

In the framework of the kinematical theory, the diffraction peak intensity is defined by the square of the ASF from all atoms (Landau & Lifshitz, 2001). In the case of XFEL pulse diffraction, it is the fluctuating value that should be averaged for all configurations,

$[\bar{R}({\bf{q}},t)\sim \left\langle \textstyle\sum\limits_a\textstyle\sum\limits_b F_a({\bf{q}})\,F_b^{\,*}({\bf{q}}) \exp\left[i{\bf{q}}({\bf{R}}_a-{\bf{R}}_b)\right]\right\rangle,\eqno(33)]$

where the symbol $[\langle\ldots\rangle]$ means the average over all configurations defined in the formula (30) and the summation is performed over the coordinates $[{\bf{R}}_a,{\bf{R}}_b]$ of the same atoms with ASF $[F_a({\bf{q}})]$ in all unit cells of the crystal.

It was mentioned above that the average ASF is supposed to be the same for all unit cells and its fluctuations are not correlated. This allows one to use the formula

$[\eqalign{&\left\langle F_a ({\bf{q}})F_b^{\,*} ({\bf{q}}) \right\rangle= \bar{F} ({\bf{q}},t)\bar{F}^{\,*} ({\bf{q}},t) + \Delta^2 F({\bf{q}},t) \delta_{ab}, \cr& \bar{R}({\bf{q}},t)\sim |\bar{F}({\bf{q}},t)|^2 \Big| \textstyle\sum\limits_a \exp\left(i{\bf{q}}{\bf{R}}_a\right)\Big|^2 \,+\, N\Delta^2F({\bf{q}},t),}\eqno(34)]$

where N is the total number of unit cells in the crystal.

The first term in (34) defines the coherent diffraction intensity in accordance with the identity (Ziman, 1972)

$[\bar{R}({\bf{q}},t)\equiv \textstyle\sum\limits_{{\bf{h}}}\bar{R}({\bf{h}},t)\delta_{{\bf{q}},{\bf{h}}} \sim N^{\,2} \textstyle\sum\limits_{{\bf{h}}}|\bar{F}({\bf{h}},t)|^2 \delta_{{\bf{q}},{\bf{h}}}.\eqno(35)]$

This value is proportional to N² and is significantly larger than the diffuse scattering background defined by the ASF fluctuations in the second term in (34).

Compared with the intensity of conventional diffraction, the change of diffraction intensity induced by an XFEL pulse $[R_0({\bf{h}})]$ can be characterized by the value

$[{{R}\over{R_0}}(N_{\rm{ph}},\omega)= {{1}\over{R_0({\bf{h}})}} \int\limits_{-\infty}^{\infty} I(t)\bar{R}({\bf{h}},t)\,{\rm{d}}t,\eqno(36)]$

that is a function of the number of photons in the pulse N_ph (photon flux) and their frequency ω; the intensity slope function I(t) is defined in (21); $[R_0({\bf{h}})]$ $[\equiv]$ $[R({\bf{h}},-\infty)]$ .

As was shown above, the standard deviation of a Bragg reflection is significantly less than the average ASF, so the expression (36) can be written as follows,

$[{{R}\over{R_0}}(N_{\rm{ph}},\omega) \approx {{1}\over{|F_0({\bf{h}})|^2}} \int\limits_{-\infty}^{\infty} I(t)|\bar{F}({\bf{h}},t)|^2\,{\rm{d}}t.\eqno(37)]$

Figs. 10 and 11 show the dependence of this ratio as functions of photon energy and fluence. Fig. 10 demonstrates that the diffraction intensity decreases in a non-linear manner if the photon energy approaches the Si K-edge. The deviation from unity is about 5% for 8 keV and reaches about 20% at 4 keV and will decrease further for energies closer to the K-edge energy. The calculation of the fluence dependence of R/R₀ at 8 keV demonstrates a dramatic drop if the fluence exceeds 1.6 mJ µm⁻² (see above).

Figure 10
Integral intensity of the XFEL pulse diffraction compared with the conventional low energy diffraction as a function of photon energy with 10¹² photons per pulse.

Figure 11
Integral intensity of the XFEL pulse diffraction compared with the conventional low energy diffraction as a function of fluence with 8 keV photon energy.

7. Estimation of the role of electron density evolution in the formation of the PXBI effect

It is also important to estimate the influence of electron density evolution on the conditions of the PXBI effect mentioned in §1. The characteristic value of the gain G for PXBI was calculated previously (Baryshevsky & Feranchuk, 1984; Leonov et al., 2013; Baryshevsky et al., 2005) for a crystal with the static X-ray polarizability and was defined by the following formula,

$[G\,\,[{\rm{cm}}^{-1}]= \left[{{Q\omega_{\rm{B}}|\chi_{{\bf{h}}}(\omega_{\rm{B}})| }\over{ |\cos2\theta_{\rm{B}}|^{1/2}}}\right]^{1/4},\eqno(38)]$

with

$[Q\approx {{2\pi e_0^2}\over{\hbar{c}}}\, {{|\chi_0|^2\hbar\omega_{\rm{B}}}\over{cE}}n_{\rm{e}},]$

where $[\omega_{\rm{B}},\theta_{\rm{B}}]$ is the Bragg frequency and angle correspondingly connected with the reciprocal lattice vector $[{\bf{h}}]$ ; n_e is the particle density in the bunch of the electrons with energy E. Threshold electron current for PXBI was estimated as j_th ∼ 10⁸ A cm⁻² if the crystal length L ∼ 1 cm.

Integral loss $[\Delta E]$ of the total energy of the bunch due to the electromagnetic interaction between relativistic electrons of the beam and the crystal with atomic charge Z_a can be estimated on the basis of the Bethe formula (Landau & Lifshitz, 1982),

$[\Delta E \approx {{4 \pi Z_a^2 e_0^4 } \over {m c^2 \Omega}}\, N_{\rm{e}} L \ln({{E}/{mc^2}}),\eqno(39)]$

where N_e is the number of electrons in the bunch.

It corresponds to the effective fluence ϒ,

$[\Upsilon \approx {{4 \pi Z_a^2 e_0^4 } \over {m c^2 \Omega S_{\rm{e}}}}\, N_{\rm{e}} L \ln ({{E}/{mc^2}}),\eqno(40)]$

where S_e is the transversal section square of the bunch.

The condition for the PXBI threshold current will be fulfilled if we choose N_e = 10⁹, S_e = 25.0 µm², E = 200 MeV, that corresponds to the electron bunches produced by the laser-driven accelerators (Nakajima, 2008; Corde et al., 2013). Then in the case of a Si crystal one can estimate from (40)

$[\Upsilon\approx 2.8\times10^{-2}\,\,{\rm{mJ}}\,\,\micro{\rm{m}}^{-2}.]$

This value is less than the threshold fluence for the essential drop of the polarizability (Fig. 11). However the secondary electrons produced by ultrarelativistic particles in the crystal may lead to substantial increase of the effective fluence in comparison with the estimation (40). It means that dynamics of the PXBI effect from the bunch should be considered taking into account the evolution of the electron density in the crystal. We suppose to analyze this case in a separate paper.

8. Discussion and conclusions

A numerical algorithm and software were developed for calculation of the X-ray polarizability of a crystal and diffraction intensity during the propagation of an intense XFEL femtosecond-pulse through a crystal. Together with photoionization and Auger processes we considered additional processes related to the free electrons generated in the conduction band of the solid state.

The results of the present paper lead to the following general conclusions:

According to Fig. 5, the role of the free electrons is dominant via the process of electron impact ionization.

According to Fig. 10, our approach remains valid for photon energies about 2 keV above the K-edge. However, in order to make the simulation more precise and avoid additional errors, the accuracy of the photoionization cross sections within the ECM should be improved for this energy region. For silicon this may happen at a photon energy of 2.5 keV. Further decrease of the photon energy will result in a decrease of the ionization potential up to the value where the ionization potential becomes deeper than the photon energy itself, where single-photon transitions from the K-shell become forbidden. Generally, the approach to the solution of the rate equations becomes invalid in close vicinity to the exact resonance. Here, one should use the density matrix method in order to take into account both diagonal and non-diagonal elements for the solution of the evolution problem. A more exact treatment in terms of quantum mechanics is needed in order to consider quantum coherence effects (Rabi oscillation) that are expected if the photon energy exactly matches the energy of transition. The coherence effects become significant only if equation (9) is non-local in time, i.e. if [Y(t-t')] has a significant time spread in comparison with the pulse duration, or, turning to the frequency domain, if $[Y(\omega)]$ is sharp in comparison with the spectral width of the pulse slope function. In this case the resonance can take place and the system of equations (10) and (12) should be used to calculate the amplitudes. However, if the frequency of the X-ray pulse corresponds to the transition to the continuous part of the spectrum, $[Y(\omega)]$ covers a wide range of the X-ray frequency that is broader than the spectral width of the pulse slope function. Then non-Markovian effects can be neglected, and we come to the rate equations in the form (11) for the occupation probabilities.

The general result of our numerical investigation consists of the predicted time dependence of the atomic form factor. As shown in Fig. 8, the ASF decreases during the propagation of the intense XFEL pulse through the crystal. This results in a drop in the diffraction intensity and in the decrease of the crystal polarizability components $[\chi_{{\bf{h}}}]$ during the pulse propagation. Due to photoionization and other processes, the amount of this drop at the end of the pulse increases if the photon energy approaches the K-resonance and can reach 50% already at 4 keV. Therefore an X-ray scattering experiment using intense XFEL femtosecond-pulses cannot probe the ground-state electron density of a crystal. Using XFEL pulses the measured ASF will always be smaller than the form factor measured with conventional synchrotron radiation. The deviation of the measured electron density from the ground-state electron density increases for photon energies closer to the K-resonance. However, major changes of diffraction intensity are expected above a certain threshold of pulse fluence. This threshold can be extracted from Fig. 11 and is supposed to be close to 1.6 mJ µm⁻² using a focus spot of ∼1 µm × 1 µm. This is remarkable because diffraction is still possible in spite of the fact that this value is much greater than those found in experiment (Hau-Riege et al., 2007, 2010 ; Chalupsky et al., 2009 ). As seen in Fig. 10, this threshold decreases with decreasing photon energy, and has to be considered using a photon energy close to the K-edge. In this case the possibility of the PXBI effect from the electron bunches should be investigated additionally.

APPENDIX A

Cross-sections calculation

We have calculated all necessary cross sections by means of an analytical model of an atom with effective charges (Triguk & Feranchuk, 2011) for each shell. The simplest one is the photo-ionization cross section that can be written in atomic units as (Santra, 2009; Landau & Lifshitz, 1989)

$[\eqalignno{\sigma^{({\rm{Ph}})}_{nl}(\omega,p)={}& {{2}\over{3}} \, {{\pi\alpha\omega}\over{p}} \, g_{nl} \sum\limits_{l_j\,=\,l\,\pm\,1} {{l_\gt}\over{2l+1}} \cr&\times \left|\textstyle\int\limits_0^{+\infty}r^3R_{nl}(r)R_{pl_j}(r)\,{\rm{d}}r\right|^2,&(41)}]$

where α is the fine-structure constant, ω is the photon frequency, p = $[[2m(\hbar\omega-\varepsilon_{nl})]^{1/2}]$ is the momentum of the photoelectron, and $[\varepsilon_{nl}]$ and g_nl are the ionization potential and the occupation number of the (nl)-subshell, respectively.

Within the framework of ECM, the hydrogen-like wavefunctions for discrete [R_nl(r)] and continuous [R_pl(r)] spectra were used,

$[\eqalignno{ R_{nl}(r)&= N_{nl} \left({{2Z_{nl}r}\over{n}}\right)^l \exp\left[-\left({{Z_{nl}r}/{n}}\right)\right] \cr&\quad \times F\left(-n+l+1,2l+2,{{2 Z_{nl}r}/{n}}\right)_{\vphantom{\Big|}},\cr N_{nl}&= {{1}\over{(2l+1)!}} \left[{{(n+l)!}\over{2n(n-l-1)!}}\right]^{1/2} \left({{2Z_{nl}}\over{n}}\right)_{\vphantom{\Big|}}^{{3}/{2}},&(42)}]$

$[\eqalignno{ R_{pl}(r)&= N_{pl}(2pr)^l \exp(-ipr) \, F(i\nu+l+1,2l+2,2ipr),\cr N_{pl}&= {{Z}\over{(2l+1)!}} \left\{{{8\pi}\over{\nu[1-\exp(-2\pi\nu)]}}\right\}^{1/2} \prod_{s\,=\,1}^l (s^2+\nu^2)^{1/2},\cr&\qquad\qquad \nu\equiv{{Z}/{p}},&(43)}]$

with Z and Z_nl being the total charge and effective charge of the (nl)-subshell, respectively. These values are calculated by a universal formula derived by Triguk & Feranchuk (2011). One should note that for numerical simulation the integrals in (41) with the functions (42)–(43), and hence the photoionization cross section, can be calculated analytically.

To calculate the electron-impact ionization cross section and consider the three-body recombination process we use the binary-encounter dipole model (Kim & Rudd, 1994 ; Kai, 2010 ) with all ionization potentials being calculated within the framework of ECM.

APPENDIX B

Justification of the validity of ECM

In order to justify the validity of ECM we performed a set of numerical simulations for the atomic carbon system. As a first test, we compared the photoionization cross sections predicted by ECM with the results obtained by the HFS model. The comparison is shown in Fig. 12; both have the same functional behavior. The slight shift between the ECM and HFS values occurs due to the fact that for reasons of simplicity we made a rough estimation for the continuous spectrum radial wavefunction.

Figure 12
Photoionization cross section for the 1s shell of neutral carbon calculated using the HFS approach (black line) and ECM (red line).

Calculations of the normalized ASF for certain atomic configurations (neutral, single core-hole and double core-hole states of carbon) by ECM are shown in Fig. 13. One can conclude that these quantities are in very good agreement with the results shown in Fig. 1 of Son et al. (2011) over the whole range of the momentum transfer values.

Figure 13
Normalized [in accordance with Son et al. (2011

[Son, S.-K., Young, L. & Santra, R. (2011). Phys. Rev. A, 83, 033402.]

)] atomic scattering factor for neutral (black line), single core-hole (red line) and double core-hole (blue line) states of carbon [the agreement with Son et al. (2011

) is very good].

Figs. 14 and 15 show the results of the simulation of the electron dynamics in carbon for the cases of pulses of 8 keV and 12 keV. Our data are in good agreement with the results shown in Figs. 2 and 3 of Son et al. (2011). However, one should note that there is a small shift between the corresponding extrema of atomic population probabilities [see Figs. 15 and 3 of Son et al. (2011)] caused mainly by the fact that we implemented smaller values for the photoionization cross sections for the reason mentioned above.

Figure 14
Time-averaged charge as a function of pulse duration: (a) 8 keV pulse; (b) 12 keV pulse. Solid lines correspond to CEIX calculation, dashed lines to data from Son et al. (2011

Figure 15
Time-averaged atomic population probabilities of the single core-hole and double core-hole states of carbon: (a) 8 keV pulse; (b) 12 keV pulse. Solid lines correspond to CEIX calculation, dashed lines to data from Son et al. (2011

The good agreement with the results of HFS for carbon justify the validity of our code for the case of silicon that will be considered below.

Acknowledgements

This work was supported by the BMBF under grant 05K10PSA. The authors would like to thank Professor Dr R. Santra and members of his group as well as Dr I. Vartaniants and Dr O. Yefanov for fruitful discussions. AL and IF are grateful to Professor Dr J. Kunz for the interest to the present research.

References

Akhiezer, A. I. & Berestetzkii, V. B. (1969). Quantum Electrodynamics. Moscow: Nauka (in Russian). Google Scholar
Altarelli, M. et al. (2007). EuroXFEL Technical Design Report, DESY XFEL Project Group, Hamburg, Germany (https://xfel.desy.de/technical_information/ ). Google Scholar
Authier, A. (2003). Dynamical Theory of X-ray Diffraction, International Union of Crystallography Monographs on Crystallography. Oxford University Press. Google Scholar
Baryshevsky, V. G., Feranchuk, I. D. & Ulyanenkov, A. P. (2005). Parametric X-ray Radiation in Crystals. Berlin/Heidelberg/New York: Springer. Google Scholar
Baryshevsky, V. G. & Feranchuk, I. D. (1984). Phys. Lett. A 102, 141–144. CrossRef Google Scholar
Batterman, W. & Cole, H. (1964). Rev. Mod. Phys. 36, 681–717. CrossRef CAS Web of Science Google Scholar
Benediktovich, A., Feranchuk, I. & Ulyanenkov, A. (2014). Theoretical Concepts of X-ray Nanoscale Analysis. Berlin/Heidelberg: Springer. Google Scholar
Bostedt, C., Thomas, H., Hoener, M., Möller, T., Saalmann, U., Georgescu, I., Gnodtke, C. & Rost, J.-M. (2010). New J. Phys. 12, 083004. CrossRef Google Scholar
Bushuev, V., Samoylova, L., Sinn, H. & Tschentscher, T. (2011). Proc. SPIE, 8141, 81410T. CrossRef Google Scholar
Chalupsky, J. et al. (2009). Appl. Phys. Lett. 95, 031111. Google Scholar
Chapman, H. N. (2009). Nat. Mater. 8, 299–301. Web of Science CrossRef PubMed CAS Google Scholar
Corde, S., Ta Phuoc, K., Lambert, G., Fitour, R., Malka, V., Rousse, A., Beck, A. & Lefebvre, E. (2013). Rev. Mod. Phys. 85, 1. Web of Science CrossRef Google Scholar
Emma, P. et al. (2010). Nat. Photon. 4, 641. Web of Science CrossRef Google Scholar
Feldhaus, J., Arthur, J. & Hastings, J. B. (2005). J. Phys. B, 38, S799. CrossRef Google Scholar
Feranchuk, I. D., Gurskii, L. I., Komarov, L. I., Lugovskaya, O. M., Burgäzy, F. & Ulyanenkov, A. (2002). Acta Cryst. A58, 370–384. Web of Science CrossRef CAS IUCr Journals Google Scholar
Gnodtke, C., Saalmann, U. & Rost, J.-M. (2012). Phys. Rev. Lett. 108, 175003. Web of Science CrossRef PubMed Google Scholar
Hau-Riege, S. P. (2011). High-Intensity X-rays – Interaction with Matter. Weinheim: Wiley-VCH. Google Scholar
Hau-Riege, S. P. (2013). Phys. Rev. E, 87, 053102. Google Scholar
Hau-Riege, S. P. et al. (2007). Phys. Rev. Lett. 98, 145502. Web of Science CrossRef PubMed Google Scholar
Hau-Riege, S. P. et al. (2010). Opt. Express 18, 23933. PubMed Google Scholar
Ishikawa, T. et al. (2012). Nat. Photon. 6, 540. Web of Science CrossRef Google Scholar
Iwayama H. et al. (2009). J. Phys. B, 42, 134019. Google Scholar
Kai, T. (2010). Phys. Rev. A, 81, 023201. CrossRef Google Scholar
Kim, Y.-K. & Rudd, M. E. (1994). Phys. Rev. A, 50, 3954–3967. CrossRef CAS PubMed Google Scholar
Kissel, L., Zhou, B., Roy, S. C., Sen Gupta, S. K. & Pratt, R. H. (1995). Acta Cryst. A51, 271–288. CrossRef CAS Web of Science IUCr Journals Google Scholar
Landau, L. D. & Lifshitz, E. M. (1982). Electrodynamics of Condensed Matter. Moscow: Nauka. Google Scholar
Landau, L. D. & Lifshitz, E. M. (2001). Physical Kinetics. Moscow: Fizmatlit. Google Scholar
Landau, L. D. & Lifshitz, E. M. (1989). Quantum Electrodynamics. Moscow: Nauka. Google Scholar
Leonov, A., Benediktovitch, A., Ksenzov, D., Feranchuk, I. & Pietsch, U. (2013). Proc. SPIE, 8778, 877816. CrossRef Google Scholar
Lorenz, U., Kabachnik, N. M., Weckert, E. & Vartanyants, I. A. (2012). Phys. Rev. E, 86, 051911. Web of Science CrossRef Google Scholar
MacDonald, W. M., Rosenbluth, M. N. & Chuck, W. (1957). Phys. Rev. 107, 350. CrossRef Google Scholar
Malgrange, C. & Graeff, W. (2003). J. Synchrotron Rad. 10, 248–254. Web of Science CrossRef CAS IUCr Journals Google Scholar
Morgan, W. L. & Penetrante, B. M. (1990). Comput. Phys. Commum. 58, 127–152. CrossRef CAS Google Scholar
Nakajima, K. (2008). Nat. Phys. 4, 92–93. CrossRef CAS Google Scholar
Neutze R., Wouts, R., van der Spoel, D., Weckert, E. & Hajdu, J. (2000). Nature (London), 406, 752–757. Google Scholar
Pellegrini, C. & Reiche, S. (2004). J. Sel. Top. Quantum Electron. 10, 1393–1404. CrossRef CAS Google Scholar
Popov, V. S. (2004). Uspekhi Fiz. Nauk. 174, 921 (in Russian). Google Scholar
PSI (2009). Ultrafast Phenomena at the Nanoscale: Science opportunities at the SwissFEL X-ray Laser, September 2009 issue, p. 88. Paul Scherrer Institute, Villigen, Switzerland. Google Scholar
Roseker, W., Franz, H., Schulte-Schrepping, H., Ehnes, A., Leupold, O., Zontone, F., Robert, A., & Grübel, G. (2009). Opt. Lett. 34, 1768–1770. CrossRef PubMed CAS Google Scholar
Santra, R. (2009). J. Phys. B, 42, 023001. CrossRef Google Scholar
Schorb, S., Rupp, D., Swiggers, M. L., Coffee, R. N., Messerschmidt, M., Williams, G., Bozek, J. D., Wada, S.-I., Kornilov, O., Möller, T. & Bostedt, C. (2012). Phys. Rev. Lett. 108, 233401. Web of Science CrossRef PubMed Google Scholar
Shastri, S. D., Zambianchi, P. & Mills, D. M. (2001). J. Synchrotron Rad. 8, 1131–1135. Web of Science CrossRef CAS IUCr Journals Google Scholar
Shvyd'ko, Y. & Lindberg, R. (2012). Phys. Rev. ST Accel. Beams, 15, 100702. Google Scholar
Son, S.-K., Young, L. & Santra, R. (2011). Phys. Rev. A, 83, 033402. Web of Science CrossRef Google Scholar
Stepanov, S. A. (undated). X-ray dynamical diffraction web server, https://sergey.gmca.aps.anl.gov . Google Scholar
Ter-Mikaelian, M. L. (1972). High Energy Electromagnetic Processes in Condensed Media. New York: Wiley. Google Scholar
Triguk, V. V. & Feranchuk, I. D. (2011). J. Appl. Spectrosc. 77, 749–757. Google Scholar
Ziaja, B., Szöke, A., van der Spoel, D. & Hajdu, J. (2002). Phys. Rev. B, 66, 024116. CrossRef Google Scholar
Ziaja, B., Chapman, H. N., Fäustlin, R., Hau-Riege, S., Jurek, Z., Martin, A. V., Toleikis, S., Wang, F., Weckert, E. & Santra, R. (2012). New J. Phys. 14, 115015. Web of Science CrossRef Google Scholar
Ziman, J. M. (1972). Principles of The Theory of Solids. Cambridge University Press. Google Scholar

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.