3D multi-scale modelling of plasma-based seeded soft X-ray lasers

Modelling plasma-based seeded soft X-ray lasers from the creation of the plasma to the propagation of a femtosecond high-order harmonic (HOH) seed throughout several millimetres of inhomogeneous plasma is a complex challenge. Different spatio-temporal scales from the hydrodynamic evolution of the plasma (hundreds of micrometres and nanoseconds) to the propagation of pulses through the plasma (millimetres and tens of picoseconds), electron collisions (picoseconds or even shorter) and the evolution of the envelope of the seeded HOH (tens of femtoseconds) must be tackled in order to fully understand these systems. In this paper, we will present the multi-scale computational paradigm that we have used to perform a full ab initio simulation of a dense, Ni-like Krypton plasma amplifier of soft X-rays. Results of the modelling and expected future applications will also be shown.


Introduction
Plasma-based seeded soft X-ray lasers [1,2] are a promising, compact source of coherent XUV and soft X-ray beams. Potential applications range from coherent imaging to warm dense matter [3] and plasma diagnosis [4]. This source has demonstrated sub-picosecond amplification of high-order harmonics (HOH) [5] and the delivery of circularly polarized light [6]. However, both experimental and modelling efforts are still required to understand and control this source and use it in a day-to-day basis.
Modelling this kind of lasers is a multi-scale problem, ranging from femtoseconds to nanoseconds and micrometres to centimetres. The lasing transition takes place in closed-shell ions (Ne-like, Ni-like, Pd-like, etc.). The amplifying medium is thus a laser-created plasma, millimetre-centimetre sized that evolves hydrodynamically in a time scale of nanoseconds. The plasma electrons are energetic enough to create the population inversion by collisional excitation in picoseconds. They also ionize the lasing ion, in scales ranging from hundreds of femtoseconds to tens or even hundreds of picoseconds, depending on the density and temperature of the plasma. Finally, a HOH with a duration of tens of femtoseconds is seeded inside the amplifier. The HOH propagates throughout the amplifier. During the propagation and amplification (tens of picoseconds), the spatio-temporal complex envelope of the beam evolves in a tens to hundreds of femtoseconds timescale. The full modelling of all these processes and related scales a e-mail: eduardo.oliva@upm.es (corresponding author) requires the use and coupling of different specialized codes optimized for each of the different processes aforementioned.
In this paper, we will present the suit of codes used to fully model plasma-based soft X-ray lasers. These codes allow us to study the different physical processes (plasma hydrodynamics, atomic physics, interaction of lasers with plasmas, absorption/amplification of XUV and soft X-rays in plasmas) and their related spatio-temporal scales (from micrometres to centimetres and from femtoseconds to nanoseconds). This computational paradigm can be used to model different systems, like plasma amplifiers from gas or solid targets or even atmospheric lasing from plasma filaments. As an example of its capabilities, we will present the full modelling of a Krypton dense plasma amplifier seeded with high-order harmonics that demonstrated the delivery of sub-picosecond pulses [5]. The amplified pulse presents some structures in intensity and phase resulting from the interaction with plasma inhomogeneities [4].
The layout of the paper will be as follows. In Sect. 2, we will describe the amplifier and its particularities, relating them to the required computational models to study it. In Sect. 3, we will describe briefly the codes used to fully model this plasma-based laser. In Sect. 4, we will show the results of the full modelling of this amplifier, from the creation and hydrodynamic evolution of the plasma to the 3D spatio-temporal evolution of the complex femtosecond envelope of the amplified HOH beam. Finally, we conclude with a résumé of the current status of this multi-scale paradigm and future applications in plasma amplifiers and plasma diagnosis. In this section, we will describe the creation of the plasma amplifier that we model with our suite of codes. The amplifier is composed of eightfold ionized Krypton (Ni-like Krypton, Kr 8+ ). This ionization stage is created through optical field ionization (OFI). The field of an intense pump laser ionizes the atoms and heats the free electrons, so they are energetic enough to pump the laser transition. The Ni-like Krypton OFI plasma amplifier was the first one to demonstrate saturated amplification of a seeded HOH [1] thanks to, among other reasons, its relative simplicity.
The plasma created with this method has a relatively low electron density (n e ≈ 10 18 cm −3 ) compared to plasmas created from solid targets. Consequently, the resulting amplified pulse had a duration of several picoseconds [5]. One strategy to reduce the duration of the amplified pulse towards attaining the subpicosecond regime consists on increasing the density of the amplifier. In denser plasmas, the electron-ion collision frequency is higher, reducing the collisional depolarization time [7,8] and thus reducing the duration of the amplified beam [9]. This collisional broadening is not enough by itself to attain the sub-picosecond regime. In this amplifier, the increase in electron density results also in an increased ionization rate. When the electron density is high enough (n e ≈ 10 20 cm −3 ), both the creation of the Ni-like ion and its disappearance due to further ionization are sub-picosecond processes. This overionization quenches the amplification in hundreds of femtoseconds, and thus, it is expected that the amplified beam will have a duration of the order of magnitude of that of the gain.
However, increasing the Krypton density to the 10 19 − 10 20 cm −3 range has some deleterious effects. In OFI amplifiers, the pump pulse is usually an intense infrared (IR) pulse that must propagate throughout several millimetres of Krypton. The electron density attained when the gas is ionized is within one order of magnitude of the critical density (which is n c ≈ 10 21 cm −3 ). The result is a strong defocusing of the pump beam that hinders its propagation. This defocusing can be overcome by tailoring the electron density gradients of the plasma. When the radial electron density profile is parabolic, the resulting refractive index gradient focuses the beam, compensating the former effect and ensuring the efficient propagation of the IR pump pulse. This is a so-called plasma waveguide, also known as plasma channel [10]. The plasma channel can be created by using a sequence of two infrared laser pulses, a short one (the ignitor) and a long one (the heater) [11,12] focused using an axicon lens. These two pulses generate a plasma and heat it. The subsequent radial hydrodynamic expansion creates the desired parabolic electron density profile in some nanoseconds [13].
Once the plasma channel is created, the intense IR pump pulse is injected inside it, propagating through-out the plasma. During its propagation, the pulse creates the Ni-like ion and heats the free electrons to ensure collisional pumping and further overionization. This propagation is, nonetheless, far from trivial. The high peak intensity of the IR laser creates an overionized region in the central part of the plasma channel that defocuses the beam. The channel subsequently focuses the beam, and the IR pump pulse propagates through a series of focusing and defocusing regions, akin to the filamentation process [14]. The resulting plasma efficiently amplifies the HOH seed in a subpicosecond scale [5], but the inhomogeneities and overionized regions induced by the oscillatory propagation process aforementioned impact the intensity and phase profile of the amplified beam [4].

Suite of codes
The complete modelling of the Ni-like Krypton dense plasma seeded amplifier requires to study different physical processes that act on different timescales, from nanoseconds to femtoseconds.
The creation of the plasma channel must be modelled with, at least, a 2D radially symmetric hydrodynamic code. The line focus created with the axicon lens requires a ray-tracing subroutine to properly model the energy deposition. Once the nanosecond evolution of the plasma is modelled, the resulting parabolic profile is fed to a Particle-In-Cell (PIC) code to model the propagation of the intense IR pump pulse throughout the channel and the creation of the lasing ion and overionized regions. Since a detailed description of collisional and radiative processes in the plasma is required, the results of the PIC simulation (mainly the propagated peak intensity) are given to a collisional-radiative code. Finally, all the data computed are fed to a 3D Maxwell-Bloch code that models the propagation and amplification of the HOH. In this section, we will describe the codes used for this task.

Plasma evolution: radiative hydrodynamics
The creation of a plasma when an intense enough IR laser interacts with a gas and its subsequent evolution during several nanoseconds is modelled using plasma hydrodynamic codes. These codes treat the plasma as a compressible fluid, enhancing these equations with other processes of interest like radiative transfer, electron thermal conduction, laser energy deposition and ray tracing.
We have used our code ARWEN [15,16] to model the creation of plasma amplifiers of soft X-rays. ARWEN is a 2D radiative hydrodynamics code which supports Cartesian and cylindrical coordinates. It makes extensive use of the adaptive mesh refinement (AMR) paradigm [17,18]. AMR allows simulations to have a uniform level of accuracy throughout the problem domain by refining the mesh only in regions where it is required (i.e. where sharp gradients appear).
In ARWEN, fluid equations are solved using a second-order unsplit Godunov method with multimaterial capability. Electron heat conduction is treated with a flux-limited diffusion solver. Radiative transfer is solved with a multi-group discrete-ordinate (S n ) synthetically accelerated radiation transport module [19].
The equation of state (EOS) that closes the system of compressible fluid dynamics is supplied to ARWEN in tabular form. Usually, it is based on the quotidian EOS [20] and fitted to experimental data [21]. Opacities are computed using our code BIGBART [22,23]. It will be briefly described in Sect. 3.3.

Propagation: Particle-In-Cell
The propagation of an intense IR beam through a millimetre or centimetre sized plasma is modelled using Particle-In-Cell (PIC) codes. These codes solve Maxwell and Vlasov equations in a self-consistent way. In these codes, the charged particles that compose the plasma are grouped in the so-called macroparticles. The trajectories of these macroparticles are computed by solving the corresponding relativistic motion equations and the Lorentz force. The charge and current are computed from this solution and fed as sources to a Maxwell's equation solver. The resulting fields are interpolated and used to compute once again the trajectories of the macroparticles.
To model the Ni-like Krypton dense plasma amplifier, we have used PIC codes, like CalderCirc [24] and FBPIC [25]. Both codes take advantage of the cylindrical symmetry of the problem and use 2D grids, with the consequent saving on computational resources and time while obtaining a good agreement with full 3D simulations. Of course, depending on the problem to model with this computational paradigm, other PIC codes with the required capabilities (for example, Kerr nonlinearity to study filamentation) might be easily used and coupled to the other codes.
Since we are interested in the relative abundance of the different ionic species of Krypton, particular attention must be given to the ionization model. For OFI plasmas, the main ionization mechanism is tunnelling ionization [26]. The alternating electromagnetic field of the laser reduces the depth of the atomic potential well in which the electron is confined. Since the tunnelling probability increases, the electron might tunnel through the barrier, ionizing the atom. The probability can be computed using the ADK model [27,28], which is the one implemented in CarderCirc and FBPIC.

Atomic processes: Collisional-radiative codes
Many plasma properties required to understand and model these systems, such as opacity, emissivity and ionization, depend on the populations of the energy levels of atoms and ions. For rapidly evolving plasmas, in a sub-picosecond scale, solving the local distributions requires the use of time-dependent models. Atomic-level populations and spectra are thus computed using codes based on the collisional-radiative model (CRM) [29,30]. These codes take into account the most relevant atomic processes that govern the species abundance. To obtain the local distributions, a system of differential equations must be solved. This system describes the rates of population and depopulation of the different charge states. The inclusion of electron-ion collisional dynamics is of utmost importance for a proper description of plasma transport properties and non-equilibrium plasma evolution.
The CRM is the preferred method to study ionization balance in plasmas out of equilibrium. These models are based on detailed atomic structure and tend to be in good agreement with highly resolved experimental data. Their key ingredients are the atomic states to be included and the rates coupling these states. In the case of complex plasmas, in order to deal with the large number of levels involved, codes based on CRM must employ an average procedure over these levels and the transition rates connecting them [31,32]. Furthermore, to speed up processing, most of the existing codes assume a thermalized distribution of free electrons when obtaining collisional rates.
The emission and absorption coefficients needed for the hydrodynamic code are computed with the atomic physics code BIGBART [33]. This code can perform local thermodynamic equilibrium (LTE) and non-LTE time-dependent atomic calculations. In this case, it is used to provide ARWEN with detailed opacity and emissivity data to compute radiative transport simulations.
For the calculation of the absorption coefficient and opacity, the contributions from all the absorbing processes with non-negligible rates are considered: resonant photoabsorption (bound-bound), photoionization (bound-free) and inverse bremsstrahlung (free-free). We solve the Saha equation, that relates the populations of charge states in LTE, and obtain the distribution of ionic states within the plasma. To simplify the calculations, the Flexible Atomic Code (FAC) package [34] is coupled to BIGBART, in order to generate atomic potentials to compute configuration wave functions for the states and radiative transitions among them.
The optical field ionization of the plasma by the intense IR pump pulse and consequent evolution is modelled with the collisional-radiative code OfiKin-Rad [35]. This code takes into account non-equilibrium effects when computing electron collision cross sections. The evolution of atomic populations is calculated along the electron energy distribution function by solving a Fokker-Planck equation. The resulting rates and atomic populations are fed to our Maxwell-Bloch model in order to compute accurately the populations of the levels involved in the lasing transition with its crude three-level model.

Amplification: Maxwell-Bloch
The propagation and amplification of the HOH beam throughout the plasma are modelled solving the socalled Maxwell-Bloch equations [36]. The particularities of the problem allow the use of several approximations that ease the algorithms and result on efficient 3D codes. The elongated geometry of the problem (the plasma is some hundreds micrometres wide and several millimetres long) defines a preferential direction (that of propagation), and the paraxial approximation can be used. This approximation decouples the propagation (an advection equation at constant velocity solved efficiently by interpolation) and the diffraction of the beam. Since the frequency of the HOH beam is much higher than the plasma frequency, the fast oscillation of the electromagnetic field can be neglected by the use of the slowly varying envelope approximation.
This equation is enhanced with a constitutive relation that defines the polarization and rate equations for the populations of the two energy levels involved in the lasing transition. These equations are obtained from Bloch equations: the polarization follows the dynamics of the non-diagonal term of the density matrix of the system, while the diagonal terms correspond to the populations [37].
These equations, along with a 3D description of the plasma and the capability to use inputs from different collisional-radiative and hydrodynamic codes, are implemented in DAGON [38], our 3D Maxwell-Bloch model. The extensive output of the code is saved in HDF format and visualized using VisIt [39].

Modelling results
In this section, we will show the results of the full multi-scale modelling of a Ni-like Krypton dense plasma amplifier seeded with high-order harmonics. As explained above, each code solves some part of the full problem and the resulting data are post-processed and fed to the other codes as initial conditions or input data.
A simple sketch of the modelled system is shown in Fig. 1. The plasma is created by two IR beams focused in a line (not depicted in the figure). The plasma expands and creates a waveguide throughout which the IR pump pulse propagates. The HOH seed arrives after some delay to be amplified in the plasma. These processes are further described in the next subsections.

Creation of the plasma waveguide
The first step consists in modelling the creation and evolution of the plasma waveguide. Experimentally [5], two infrared laser beams (λ = 800 nm) are focused by an axicon lens in a 5 mm ×500 µm nozzle with a backing pressure of 150 bar. These two pulses, the ignitor (130 mJ, 30 fs) and the heater (690 mJ, 600 ps), create and heat the plasma. The subsequent hydrodynamic expansion leads to the formation of the plasma waveguide.
In order to model this problem, we take advantage of its symmetry and use cylindrical coordinates (r − z) in ARWEN. The z direction corresponds to the length (5 mm) of the nozzle. The 5-mm line focus created by the axicon lens is modelled using the 3D ray-tracing subroutines implemented in ARWEN. Figure 2 shows the radial electron density profile at two different times: 1.55 ns (upper panel) and 2.55 ns (lower panel) after the arrival of the ignitor pulse. Figure 2 shows that a parabolic density profile, capable of guiding an intense IR pulse, is achieved. Moreover, its radius and density gradient vary with time.
(The radius increases with time, due to the expansion, while the radial gradient decreases with time.) Thus, it is possible to match the plasma waveguide with the intense IR pulse varying its injection time [14].

Propagation of the intense IR pulse through the waveguide
Once the plasma waveguide is created, an intense (I = 3×10 18 W · cm −2 ) IR pulse is injected in the waveguide 1.55 ns after the arrival of the ignitor pulse (Fig. 2 upper  panel). We model the propagation of the IR pulse with PIC codes (FBPIC).
The simulation window spans 70 µm in the axial (propagation) direction and 30 µm in the radial direction. The corresponding mesh is 1600 × 200 cells in each direction, respectively, with 2 macroparticles per cell and per direction. The electron density profile of the waveguide is given by the hydrodynamic simulations aforementioned. In order to simplify the coupling between codes, the radial shape is approximated using analytical formulas. Higher-order approximations, by using Padé approximants, for example, can be easily implemented if required. In our simulations, a parabolic radial profile has shown to be an accurate enough approximation. In addition to this, we assume that Kr 3+ is the predominant species in the waveguide, as hydrodynamic simulations suggest. Since Kr density at r = 0 µm (the valley of the parabola) is n 0 = 1.5 × 10 19 cm −3 , assuming a threefold ionized plasma results in a minimum electron density of n e = 4.5 × 10 19 cm −3 and of n e = 1.2 × 10 20 cm −3 when the lasing ion (Kr 8+ ) predominates. These densities are high enough to be within one order of magnitude of the critical density of the IR pulse n c ≈ 10 21 cm −3 , and the propagation velocity of the IR pulse is reduced. Indeed, in order to properly track its propagation, the simulation window must move with a velocity slightly lower than that of the light in vacuum. The upper panel depicts the initial regions of the plasma (370 µm ≤ z ≤ 430 µm). At the centre of the plasma (r = 0 µm), there is no Kr 8+ since the beam is intense enough to create higher ionized species (Kr 9+ , Kr 10+ ). This overionization in the central region produces an excess of electron density that defocuses the laser beam. At higher radii, the beam intensity has decreased to barely sustain the creation of Kr 9+ or even only Kr 8+ , being the latter the most abundant in these regions.
After propagating some millimetres, as depicted in Fig. 3 central panel, the intensity of the beam has been reduced, so Kr 8+ is abundant now at lower radii regions, while at higher radii the intensity of the laser beam is so low that no Kr 8+ is produced. An overionized region in the central part of the plasma is still observed. This is no more the case when the beam has propagated four millimetres, as shown in Fig. 3 lower panel. The lasing ion only appears in the central region of the amplifier (r ≤ 5 µm).
This overionization induced defocusing competes with the guiding effect of the plasma. A substantial part of the initial energy is lost in the first millimetre of the plasma due to this effect. Afterwards, the energy is lost due to further absorption and ionization in the plasma and the mismatch between the beam spatial profile and the waveguide.
PIC modelling allows unveiling, on the one hand, the deleterious effect that overionized regions have in the propagation of the IR pump pulse, defocusing and depleting it, reducing the intensity of the beam and thus the plasma volume where lasing ion is produced. The resulting amplifier is strongly inhomogeneous, with central regions presenting electron density peaks and low abundance of lasing ions during the first millimetres of the amplifier. These regions impact the amplification of seeded HOH [4]. On the other hand, the data show that the lasing ion is produced during all the amplifier length, thus ensuring an efficient amplification of the HOH seeded throughout several millimetres of plasma.

Collisional pumping and further collisional ionization
As explained, the IR pump pulse creates the lasing ion by optical field ionization (in the regions where the  intensity is adequate). The resulting ionized electrons are energetic enough to pump the lasing transition in Ni-like Krypton (3d 9 4d J=0 → 3d 9 4p J=1 ) and to further ionize Krypton in a hundreds of femtoseconds timescale.
The collisional-radiative code we are using, OfiKin-Rad [35], requires as inputs the propagated intensity of the IR laser and the initial ionization. We obtain these values from the PIC and hydrodynamic modelling, respectively. We use as input intensity a value of I = 10 17 W · cm −2 . PIC simulations show that the propagated intensity is lower than the initial intensity of the laser beam, due to the overionization induced refraction that takes place at the centre of the waveguide during the first millimetre of propagation and subsequent loses due to the mismatch between the pulse and the waveguide. The propagated beam has thus an intensity high enough to produce Kr 8+ , maintaining the parabolic radial profile of the waveguide, but not high enough to overionize the central region. Figure 4 shows the temporal evolution of the gain (upper panel) and mean ionization (lower panel). As it can be seen on the inset, the gain has a full width at half maximum (FWHM) duration of 260 fs. This shorttimed gain is explained by the so-called collisional ionization gating mechanism [5]. As shown in the lower panel of Fig. 4, the IR pump pulse immediately creates the lasing ion (Kr 8+ ) by optical field ionization. However, the electron fluid is energetic and dense enough to collisionally ionize the lasing ion in a timescale of hundreds of femtoseconds. When the abundance of Kr 8+ is low, the gain is shut off. It is worth mentioning that, in addition to the ionization and gain, OfiKinRad computes the populations and related collisional and radiative rates of 92 levels of the lasing ion, along with other quantities of interest as the electron density and temperature. Dagon, our Maxwell-Bloch code, solves a simplified three-level atomic model. To ensure its accuracy, the time-dependent rates previously computed by OfiKin-Rad are tabulated and fed to Dagon. Using these rates in the three-level model ensures that its results (mainly, the population inversion) are close enough to those of the full collisional-radiative model while requiring much less computational time.

Amplification of high-order harmonics
The propagation of the IR pump pulse creates the lasing ion and energetic electrons. These electrons collide with Krypton ions, exciting them and creating a population inversion between the 3d 9 4d J=0 and 3d 9 4p J=1 levels. The gain is quenched several hundreds of femtoseconds after its creation due to further collisional ionization that depletes the lasing ion abundance. A HOH with wavelength tuned to match that of the transition is seeded in the plasma to be amplified.
The amplification process in this kind of dense OFI plasma amplifiers is quite complex. The plasma is seeded with a λ = 32.8 nm, 1 nJ, 100 fs FWHM Gaussian HOH. The delay between IR pump pulse and HOH is 200 fs. As expected, the HOH starts to generate a wake that is strongly amplified [40], although the strong collisional broadening, due to the high density of the amplifier, reduces the duration of this wake [9]. Since the gain duration is less than 500 fs (as shown in Fig. 4 inset), the amplified HOH should have a duration and envelope matching that of the gain. However, the velocity mismatch between the IR pump pulse (and consequently of the gain generated by this beam) and the HOH results in the head of the HOH overtaking the IR beam. The resulting HOH pulse has a complex spatial shape, as shown in Fig. 5. In addition to this effect, the inhomogeneities in ionization and electron density contribute also to distort the envelope of the HOH [4].
Still, the resulting amplified HOH has a subpicosecond duration, an energy of the order of 1.5 µJ and a good wavefront, suitable for applications. Moreover, the amplified HOH can act as a probe of the plasma amplifier since the distortions observed in the intensity profile and wavefront are related to the lasing ion abundance and electron density.

Conclusions and future prospects
In this paper, we have presented the multi-scale computational paradigm we are using to fully model plasmabased seeded soft X-ray lasers. This paradigm tackles the variety of problems and spatio-temporal scales of these lasers, from nanoseconds (hydrodynamics) to picoseconds (propagation of IR pulses, collisional processes) and femtoseconds (gain duration, HOH envelope dynamics).
As an example of the capabilities of our computational paradigm, we have modelled a dense plasma amplifier of Ni-like Krypton. The different physical processes involved have been studied ab initio, showing a good agreement with experiments [4,5] and unveiling the structure of the plasma and the amplified HOH [4].
Since this model has been validated for plasma-based soft X-ray lasers, it will be possible to apply it in the near future to model transient collisional excitation [41] and quasi-steady-state [42] X-ray lasers or even more complex problems, as plasma amplification chains [43][44][45], recombination lasers [46], the interaction of structured beams with plasmas and the diagnosis of warm dense matter with HOH [3]. In addition to this, we plan to apply this model to study lasing from atmospheric filaments.

Author contributions
All authors contributed to the paper.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Data Availability Statement
This manuscript has no associated data, or the data will not be deposited. [Authors' comment: The data used is freely available upon request.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.