Macroscopic laser–plasma interaction under strong non-local transport conditions for coupled matter and radiation

High-intensity laser (⁠$≳1018W/cm2$⁠) interaction with a solid target has been studied extensively in the last decades, experimentally and numerically. Even though the parameters of the main pulse can be finely adjusted, less attention is paid to accurate predictive analysis of the pre-pulses appearing in these experiments. In general, pre-pulses are inevitable precursors of high-intensity main laser pulses, where especially amplified spontaneous emission (ASE) contributes significantly. As a consequence, a pre-plasma is created in front of the target, which may affect the interaction of the main laser pulse by itself or by means of the caused disruption or deformation of the target. These effects have been observed and proved their significance in proton acceleration by irradiation of thin foils [1–3], fast-ignition schemes [4,5] and also interaction of fs-pulses with high-Z targets [6,7]. Moreover, a recent paper investigates the pre-pulse effect on a $1023W/cm2$ peak intensity main pulse interaction [8].

This article focuses solely on the pre-pulse interaction with a solid target. A comparison is made as far as the material is concerned, where moderate-Z materials are represented by Aluminium and high-Z ones by Copper. Simulations for laser intensities of $1013W/cm2$⁠, $1014W/cm2$ and $1015W/cm2$ were performed. The hydrodynamic code PETE [9,10] is used for this purpose, which includes a state-of-art non-local radiation and electron heat transport model [11,12]. Even though the classical diffusion approximation is invalid for such high laser intensities, this fact is normally disregarded in simulation codes. Here, a direct comparison of simulations with the non-local radiation transport and radiation diffusion is presented, where 12 different simulations were made in total. The effects on formation of the pre-plasma then becomes evident especially in the regime of tightly coupled radiation and matter.

The rest of the paper is organized as follows. The following section, Sec. 2, briefly discusses the physics issues of non-local transport in the context of laser–target interaction. Sec. 3 describes applied numerical methods and newly developed schemes. Sec. 4 is dedicated to the virtual experiments of laser–target interaction and pre-plasma creation. Finally, discussion of the results and conclusions are given in Sec. 5 and broad outlook for upcoming research in this field is made in Sec. 6.

The diffusion approximation of energy transport has been prevalent numerical method due to its simplicity and numerical stability for many years. Even nowadays, it is still not uncommon to employ flux limiter models which limit the conduction flux $qe=−κe∇Te$ to a fraction $fe$ of the free-streaming flux $qf=nemevTe3$⁠, where $Te$ and $ne$ stand for the electron temperature and density, $me$ electron mass and $vTe$ is the electron thermal velocity $vTe=kBTe/me$ with $kB$ being the Boltzmann constant. Typically, the conductivity is modified by using the formula:

Here $κsh$ is the Spitzer-Härm conductivity and $κe$ the supposedly corrected one. This approach to handle non-locality makes no sense as it depends on the local temperature scale length $Te/|∇Te|$ whereas it is the global temperature scale length $k−1$ (with k being the spatial frequency) which determines the heat flux. A reduction of the heat-flow for large $kλe$⁠, where $λe$ is the electron mean free path, can therefore never be modelled by a flux-limiter approach.

A detailed kinetic analysis of the phenomenon of non-local transport [13] shows that the bulk collisional electrons are responsible for absorption and increase in electron temperature but contribute little to the heat flux whereas the hot collisionless electrons hardly affect the electron temperature but dominate the heat flux. If no distinction is made between cold and hot then non-locality and flux inhibition arise if the real heat flux is represented on a macroscopic level as a conductivity times a temperature gradient. Non-local transport is a consequence of macroscopic models not taking account of the specific form of the distribution function.

Non-locality in the transport coefficient can compensate for the loss of detailed information of the underlying kinetics as single-temperature fluid-type models do not distinguish between cold bulk electrons and hot tail electrons. The original approach to non-locality was to present the heat flux as a convolution over space in the form:

with a certain appropriate kernel for the conductivity [14–16]. This accounts for the strong variation of the conductivity as function of space. However, this integral approach is in general an ill-posed and intrinsically unstable problem [17], although these problems can be overcome with a corresponding overhead as far as the numerics is concerned [18]. The impossibility to extend the convolution models to more dimensions was overcome later in [19], where a linear transport equation was found for the problem. However, a direct physical reasoning behind the models remained still elusive. An answer to this principal question was later given by renaissance of Krook-type operators [20].

A more detailed analysis can be found in our recent paper [12], which was dedicated to the problem of non-local electron heat transport, formulation of a model based on a Krook-type operator and development of an efficient numerical method for this purpose. However, the problem of non-locality is not inherently connected only with electrons, but with any kind of particle in principle, provided that it can be delocalized sufficiently to escape local potential values. In such case, classical perturbation theory is not applicable and the Chapman–Enskog method [21] underlying the diffusion approximation cannot be used for expansion of the distribution function.

This situation also occurs in radiation transport, which is the main topic of this paper. Diffusive approximation is broadly applied in this area too [22,23], where radiation flux limiters are introduced into the diffusion operator in a fashion similar to (1), i.e. modifying the diffusion operator:

where $f=f(KnR)$ is the radiation flux limiter and c is the speed of light in vacuum, $κR$ Rosseland opacity and $KnR=|∇ϵR|/(κRεR)$ the radiation Knudsen number for the densities of radiation energy $ϵR$⁠. In principle, all flux limiters satisfy the diffusion limit $f(KnR)≈1$ for $KnR≪1$ and the non-local limit $f(KnR)≈3/KnR$ for $KnR≫1$⁠, where the resulting flux approaches the free streaming value $|qR|=cϵR$ with $ϵR$ being the radiation energy density. The similarity with the heat conduction is apparent here, since the flux ratio $|qR|/(cϵR)∼KnR$⁠. However, the flux limiters are still dependent on the local profiles of temperature, which can be irrelevant in many cases as manifested in the test from Sec. 3.3.2.

From the physical point of view, it can be proved that flux limiters are related to the variable Eddington factor (VEF) [24]. This connects the flux limited diffusion (FLD) with the Eddington approximation that poses another popular approach in the field. An expansion of the radiation transport equation into angular moments is made, where the first two moments are typically considered to be sufficient [25]. Termination of the expansion at the energy equation is achieved by construction of an appropriate closure independent of the higher moments. Using only the isotropic moments for construction of the closure, the VEF method is obtained, as originally proposed in [26]. Later, the method was improved by taking the radiation fluxes into account as suggested in [27]. Several examples of flux limiters directly corresponding to the closures are given in Sec. 3.3.1.

Another widely applied approach is the discrete ordinates (⁠$SN$⁠) method, where the quantities are discretized in the angular space. Following then the technique of short-characteristics, the radiation transport is solved on the level of energy exchange between adjacent computational cells. This method was originally developed by Carlson [28] and become widely known thanks to Pomraning [29]. It brings in essentially the true notion of non-locality of the energy transport.

Sometimes, a combination of diffusion with a $SN$ method can be encountered in astrophysical codes [30], where the latter model covers only the free streaming limit of the transport. However, in the context of laser–target interaction, radiation transport occupies the space between the two limits parametrically as was observed in Sec. 4.2 and a single numerically effective method with correct limits is necessary. We presented such a method in [11] and its brief description is given in Sec. 3.3.

The main aim of this text is then to show the discrepancies between the diffusive and non-local transport models in real scenarios involving interaction of high-intensity lasers with a solid target. An intensity scan is made as well as a scan in the material of the target, because all of these parameters strongly affect the results. Moreover, the role of non-local energy transport and its interplay with the hydrodynamics is also discussed. In order to clarify that the principally correct description including non-locality of the transport and self-consistency of the laser absorption must be applied, some trade-off in overall realism of the simulations was made to keep the simulations illustratively simple, but still covering important phenomena.

Detailed radiation-hydrodynamic simulations of the ablation process is of fundamental importance for inertial confinement fusion (ICF) as well as pre-pulse characterization of high-power laser–matter interaction. In particular, the following issues require more detailed understanding:

• extend and profile of the sub-critical plasma corona in front of the target,

• location of the critical density layer $nc≈1021λ02[μm]cm−3$⁠, with $λ0$ the laser wavelength in vacuum,

• velocity and position of the created shock wave.

The macroscopic, single-fluid, two-temperature equations of radiation hydrodynamics governing laser–plasma interaction in the Lagrangian frame of reference take the following form [9,10,12]:

where ρ is the mass density, $u$ single fluid velocity, $pe,pi$ are electron and ion pressures respectively, $εe,εi$ represent the specific internal energies, $Te,Ti$ temperatures of the species, $εR$ stands for the density of radiation energy and G is the energy exchange coefficient. Finally, $qe,qi,qR,qL$ are energy/heat fluxes of electrons, ions, radiation and laser field (i.e. the Poynting vector). Note that the electron and radiation heat fluxes $qe,qR$ represent both, the non-linear heat/radiation diffusion fluxes or the non-local heat/radiation fluxes respectively, depending on the configuration. The radiation transport models are detailed in Sec. 3.3. The ion heat flux $qi$ is obtained only from the diffusion approximation (⁠$qi=−κi∇Ti$ with $κi$ being the ion heat conduction coefficient), since the ion Knudsen numbers $Kni=λi|∇Ti|/Ti$ with the ion mean free path $λi$ are typically very small (⁠$Kni≲10−3$ in our simulations) and diffusive approximation is then sufficient. However, the formulation of the energy equations (6) and (7) is only formal, because both closure models (non-local radiation/heat transport and heat diffusion) implicitly solve the isochoric heating of the electrons through these equations. The compression parts of (6) and (7) are explicitly calculated together with the rest of the hydrodynamic equations (4) and (5) by the compatible hydrodynamic scheme on a staggered computational mesh [31].

Full details of the non-local electron models are out of scope of this text and can be found in [9,12], but some of the key features must be highlighted. The model is based on the first principles, where no arbitrarily chosen function enters the model as it is common for the convolution models mentioned in Sec. 2.1. It is essentially built on the Bhatnagar–Gross–Krook type collision operator [32]. Its applicability in the context of the non-local transport for macroscopic models was widely analysed in [33]. Moreover, angular distribution of transported species is retained without any limiting assumption on anisotropy of the distributions. Furthermore, the formulation itself allows construction of discretization fully satisfying energy conservation, since it naturally uses energy intensities in a very similar way to the radiation transport model described in Sec. 3.3. Finally, the energy exchange term $G(Ti−Te)$ is implicitly included to provide proper coupling of the electron and ion equations (6) and (7).

The laser absorption model determines the laser energy deposition, which appears as the explicit source term $∇·qL$ on the right-hand-side of the equation of energy conservation (7). In order to solve the laser propagation in an inhomogeneous medium, the algorithm based on stationary Maxwell's equations is employed. It was originally introduced in [34] and later extended and successfully applied in [35,36], but never detailed thoroughly. In this text, a complete description of the method is presented, including a brief inference of the numerically solved equations. The main advantage of the model can be seen in its full self-consistency, where only an appropriate closure model for the collision frequency must be supplied as described in Sec. 3.2.2, but no empirical or ad hoc reflectivity coefficients of the critical plane are needed.

The method originates from the Helmholtz equation for transverse component of the electric field $EL$ in 1D (or more precisely the phasor of the harmonic field in time with angular frequency ω) that takes the form:

where c is the speed of light in vacuum. Note that polarizations of the waves are not distuinguished here, since only normal incidence of a lineraly polarised laser wave is assumed. The relative complex permittivity ε reads [37]:

where $ωpe=e2ne/(ε0me)$ stands for the electron plasma frequency with e the elementary charge and $ε0$ the permittivity of vacuum. The electron collision frequency $νe$ determines the efficiency of inverse-bremsstrahlung laser absorption and the closure model for it is given in Sec. 3.2.2.

In order to proceed, it must be taken into account that the numerical method uses staggered spatial discretization, where ε is approximated by a piece-wise constant function on the computational cells. The solution in the interior of the cells then consists of a superposition of the inward $P=P0exp(−iωcεz)$ and the outward $R=R0exp(iωcεz)$ propagating plane waves (⁠$EL=R+P$⁠). The reflection coefficient of the waves is defined as $V=R/P$⁠. This allows to cast (8) into an equivalent system of first-order equations for V and P [34]:

where $k=ωcε=k0ε$ is the wave number.

However, a full solution of (10) in terms of the electric field $EL=P(V+1)$ is not needed, since the point of interest is only the time-averaged Poynting vector $qL$⁠. Consequently, information about the phase of P can be removed completely, avoiding ad hoc determination of the phase in the boundary condition for P completely. Following this idea, the coefficient $A=14|P|2$ is introduced. These two coefficients, V and A, are then chosen as the primary variables and the system (10) is transformed into:

where $χ=Imε$ is the imaginary part of the complex refractive index. It must be stressed at this point that (11) is no longer equivalent to (10) and to the original Helmholtz equation (8) consequently. The irreversibility of the transformation is implied by the definition of the second order quantity A. Hence, the resulting equation for A is real and (11) is composed of only two linear first-order ordinary differential equations (one complex and one real).

The differential equations (11) are integrated inside the cells and interface relations are applied at the computational nodes, which represent optical step interfaces in the given discretization. Hence, the left and right values of the quantities relative to the interface are distinguished and denoted with a $−$ and $+$ subscript respectively. After some algebra, one may derive from the interface conditions for electric field $P−+R−P++R+$ and magnetic field $ε−(R−−P−)=ε+(R+−P+)$ the expressions for V and A in the form:

In order to solve the resulting global system of equations, two boundary conditions are set. The first one is based on the assumption that V vanishes behind the critical density plane $zc$ (i.e. the point, where $ω=ωpe$⁠) and sets $V(zmin)=0$⁠. The position of the point $zmin$ is given by the requirement $∫zminzcχdz>Czminλ0$⁠, where the constant $Czmin$ is typically 2–4. The other boundary condition originates from the fact that A coincides with the definition of Poynting vector in vacuum and it states $A(zmax)=IL$⁠, where $IL$ is the intensity of the incoming laser beam.

Finally, the solution is numerically obtained for V using (12) at the interfaces and the analytic solution of (11) in the form $V∼exp(2ikz)$ in the interiors of the cells, starting from the boundary $zmin$⁠. Then, the solution of A can be constructed from the opposite boundary $zmax$ utilizing the Eq. (13) at the interfaces and the solution $A∼exp(2k0χz)$ within the cells. After the solution is found for A and V, the Poynting vector $qL$⁠, which is of the main interest here, is obtained in the form:

where $n=Reε$ is the real refractive index. Finally, the resulting Poynting vector $qL$ of the laser field is substituted into (7), satisfying the power conservation of the laser field and absorbed power in the plasma this way.

Viability of the method was verified by a test given in  Appendix A. The numerical solution is compared to the analytic one for a simplified problem resembling the typical configuration of absorption in vicinity of the critical plane.

The previous section described the basic model of laser absorption, but several relevant issues occurring in the resulting numerical method have not been mentioned. Their appearance can be seen as a direct consequence of the made physical approximations. First of all, the wave nature of the Helmholtz equation (8) and consequently the complex equation for V (11) requires spatial resolution of the mesh comparable to the wavelength of the laser wave. In particular, we obtain from the Nyquist–Shannon sampling theorem [38] the requirement $Δx≤122π/[2Re(k)]=14λ0/n$⁠, where $Δx$ is the spatial step and $λ0$ is the vacuum wavelength of the laser. In practice, the resolution needs to be even better to attain a good accuracy of the solution. An additional limit comes from the fact that the piece-wise constant approximation of ε is applied, which implies the condition $Δx≤|k|/|∇k|∼|ε|/|∇ε|$⁠. Therefore, it is not feasible to construct the main Lagrangian computational mesh with typically sub-micron spatial step that must be guaranteed throughout the whole simulation, despite the fact that the plasma corona can reach even centimetre scale. Hence, an adaptive mesh refinement is performed instead within each time step, where the main mesh is refined to satisfy both mentioned requirements. In each computational cell, linear reconstruction of the hydrodynamic quantities is performed and the Barth–Jespersen limiter is then applied to avoid creation of new local extremes [39]. However, it must be stressed that the method maintains conservation of the quantities and this property becomes essential for steep gradients especially. Moreover, local linear regression over one wavelength is applied instead of point-wise projection to remove the periodic component of the solution in the case the main mesh does not satisfy the given requirements any more. This procedure prevents occurrence of aliasing effects in turn.

An additional level of refinement is provided by mapping of the temperature finite elements used within the high-order DG-FEM scheme for non-local electron heat and radiation transport (see Sec. 3.3). However, a full description of the method is out of scope in this text and can be found in [10]. It is a possible topic of future research to develop a higher order laser absorption algorithm in order to increase the computational efficiency and partially ease the requirements on the spatial resolution.

For purposes of laser absorption, we use the global model of the electron collision frequency $νe$ proposed in [40], which covers a wide range of temperatures and densities during the pre-plasma formation. This point is essential in early times of the simulation, when the laser irradiates a solid target and a dense plasma is being formed in front of the surface. The effects of this abrupt evolution should be included in accurate simulations of laser–target interaction in order to properly capture timing of the shock wave launching [41].

In the case of a hot plasma the model is identical with the classical Spitzer-Härm formula [42]:

The Coulomb logarithm $Λ$ is approximated by the expression $Λ=max(2,ln1+(bmax/bmin)2)$ [43], where $bmax=kBTe/memax(ωpe,ω)$ and $bmin=max(Z¯e2kBTe,ℏmekBTe)$ with $Z¯$ being the mean ionization and $ℏ$ the reduced Planck constant. In contrast, electrons are degenerated in the solid state limit and dependency on $Te$ vanishes. Scattering of electrons by phonons must be taken into account [44] and the collision frequency is approximated in this case by [40]:

Here, $vF$ denotes the Fermi velocity $vF=ℏ3π2ne3/me$ and $kS$ is an empirical constant used for adapting the model to experimental results. Both frequencies, (15) and (16), are then interpolated by the harmonic mean $νe−1=νei−1+νel−phonon−1$⁠. This formula satisfies both limiting cases, but it strongly overestimates collision frequencies for intermediate temperatures near the Fermi temperature, where none of the previous formulas (15) and (16) hold [45]. A reasonable criterion is then added to avoid the cases, when the mean free path is shorter than a characteristic inter-ion distance [40]:

where the characteristic electron velocity $ve$ is approximated by $ve=vF2+vTe2$ to satisfy both limits of a hot plasma and a cold solid.

The presented collision frequency models are compared in Fig. 1 for Copper. The value of $kS$ was determined by the procedure described in [40] based on reflectance of the material. It has value $R≐0.41$ for the parameters used in the plot [46], which yields $kS≐396$⁠. All dependencies are artificially prolonged outside the area of validity to show their limit behaviour. The limitation of the collision frequency in the intermediate regime given by (17) is significant. However, the effects on evolution of the simulated system are minimal in real scenarios, since the system stays outside this region most of the time due to simultaneous expansion and heating of the matter.

Radiation transport has an essential role in the interaction of high intensity lasers with a solid target, especially for high-Z materials. In the context of the current work, we limit ourselves to only gray body approximation with opacity scaling laws taken from [47], where the average ion model was used to calculate the LTE electronic levels occupancies. The reason for these simplifications is that we want to clearly show the effects of non-locality without additional considerations about the spectral distribution. However, it is evident that full physical realism can be achieved only with proper treatment of the spectra independently of the followed (non-)local approach [48].

As it is customary, radiation transport is described by the mean radiation intensity $IR=IR(t,x,n)$⁠, depending on the time t, spatial coordinate $x$ and the transport direction unit vector $n$⁠. The classical gray-body radiation transport equation with linear scattering reads [29]:

where $j=j(t,x,n)$ represents the mean emissivity, $κ=κ(t,x)$ is the average opacity (see below) and $σs=σs(t,x,n)$ stands for the mean scattering coefficient. The symbol $I¯$ is the angular average of the intensity, i.e. $I¯=1/(4π)∫4πIdn$⁠. Finally, it must be noted that the convection term for intensity is meant in the ”physical” sense, so it is defined as the directional derivative $n⋅∇I=∂I/∂n$ in order to be applicable in all coordinate systems [49].

Numerically, the Eq. (18) is solved by a generalization of the discrete ordinates method (typically denoted as $SN$⁠) using the high-order Discontinuous Galerkin Finite Element Method (DG-FEM). The finite element method is applied in both, the spatial dimension and the angular dimension [50]. Details about the discretization and construction of the numerical method can be found in [9,11]. At this point, we would like to stress that the method includes an implicit coupling of the radiation and the matter, that means inclusion of the isochoric part of the electron energy equation (7) into the radiation transport equation (18). This feature is not usually seen in numerical codes and we find it important for establishing accurate radiation wave profiles and stability of the simulation. Finally, the intensities are integrated over all solid angles to obtain the resulting radiation energy fluxes as $qR=∫4πIndn$⁠, which are then substituted into the electron energy equation (7). We will refer to this model as the non-local radiation transport, since it inherently includes non-locality in the description. A concurrent option is to compute the radiation energy fluxes $qR$ from radiation diffusion instead, as described in the dedicated Sec. 3.3.1. However, both methods are validated to approach the very same limit in the diffusive regime, see Sec. 3.3.2. Even from the numerical point of view, the non-local transport scheme converges very well for the 3rd order polynomial approximation of the radiation intensities in space and polar angle (split in forward and backward direction) used solely in this text, as it was partially shown in [11]. Details of the numerical features will be presented in more detail in subsequent papers.

The physical scenario of the laser target interaction allows to apply additional simplification as far as the radiation transport equation (18) is concerned. The mean emissivity is set equal to $j=κσ/πTe4$⁠, where σ is the Stefan–Boltzmann constant, and the scattering is neglected completely [51]. Finally, the notion of the average opacity must be explained. When Planck spectrum is assumed, it is widely known that the mean opacities in the gray-body approximation approach the Rosseland mean opacity $κR$ in the diffusive limit (see Sec. 3.3.1). On the other hand, in the non-local limit the Planck mean opacities $κP$ are more appropriate in this case [49]. The regime in between is very complicated from the point of view of modelling and a multi-group treatment is more adequate. However, an interpolation must be made between the two limits in the gray body approximation. Such an ad hoc interpolation was proposed in [48] and its good performance was verified on problems of radiative shocks. The average opacity κ then reads:

where $fR$ is the reduced radiation flux $fR=|qR|/(cεR)$⁠. Since this factor has the limiting values $fR→0$ for the diffusive regime and $fR→1$ in the free-streaming limit, it presents a very natural choice for parametrization and is equivalent to the radiation Knudsen number$KnR=|∇εR|/(κεR)$ for radiation diffusion in fact. A test problem dedicated to validation of the model was performed in Sec. 3.3.2. However, the model is only empirical without any direct physical justification. Hence, we propose this problem as a topic of our future research.

The equation of radiation transport (18) can be simplified even further by taking into account the assumptions mentioned in the preceding paragraphs and present only a small deviation from the isotropic local distribution. Moreover, the given medium is optically thick or the gradients are small, i.e. the radiation Knudsen number $KnR$ is small. Then truncation after the first two terms of spherical harmonics expansion and stationarity of the fluxes yield the non-equilibrium diffusion approximation [29]:

where $cV=(∂εe/∂Te)ρ$ is the specific heat of electrons and $a=4σ/c$ presents the radiation constant. As customary for diffusion models, an appropriate flux limiter $f=f(KnR)$ must be applied in cases when the requirement on smallness of $KnR$ is not fully met.

The flux limiters must respect the diffusive limit $f(KnR)→1$ for $KnR→0$ and the free streaming limit $f(KnR)→3/KnR$ for $KnR→∞$⁠. Although the limits are known, the most important regime in the context of laser plasma is about unity as can be deduced from the simulations in Sec. 4. In this region, the flux limiters may differ as shown in Sec. 3.3.2. Several flux limiters based on $KnR$ have been implemented in the simulation code PETE. Simpliest one is the flux limiter denoted as ”sum”, which is defined as $fsum=1/(1+KnR/3)$⁠. Another choice presents so-called Larsen limiter [25]:

where $n=2$ was chosen for the simulations in Sec. 3.3.2 corresponding to the optically thick limit. For purposes of the virtual experiments in Sec. 4, we use the limiter developed by Levermore and Pomraning [52]:

Another popular choice is the limiter developed by Minerbo [53], that is based on anisotropy considerations, and reads:

The list of the flux limiters is certainly not complete and many different formulas can be found in the literature, especially for astrophysical problems (see e.g. [23] and references therein). However, we have not observed any significant differences for our physical scenario, but all flux limiters are compared theoretically in Sec. 3.3.2.

Regarding the system (20 and 21), $εR$ in the radiation diffusion equations can be replaced by the radiation temperature, that is defined by $εR=aTR4$⁠. It fully characterizes the diffusion description and the system is then typically denoted as 2-T (two-temperature) or 3-T (three-temperature) including the ion energy equation (6). However, $TR$ serves only as an auxiliary parameter in the case of the non-local radiation transport.

Numerically, the system of equations (20 and 21) is solved by a scheme based on the method of mimetic operators [54]. The full inference of the scheme is out of scope of this text, only a brief description is given. First of all, the quantity $θ=aTe4$ is introduced. Then temporal discretization with backward Euler method in (20) is performed with the fixed opacities $κR$⁠, $κP$ and specific heats $cV$ at the time level n giving rise to the semi-discretized implicit system:

where $μ=cκP$ is a relaxation parameter, the symbol $τ=ρcV/(4aμ)(Ten)−3$ is the characteristic time of relaxation, the diffusion coefficient is defined as $d=fc/(3κR)$ and $Δt$ is the time step size.

The linear ordinary differential equation (26) can be integrated numerically or analytically with parametric dependence on $εRn+1$⁠. In both cases, the solution can be written in the form $θn+1=βθn+(1−β)ϵRn+1$ [22]. The fully implicit time differencing gives the value of the parameter $βBE=1/(1+Δt/τ)$⁠. In contrast, the analytic solution yields $βexp=exp(−Δt/τ)$⁠. Both coefficients give exactly the same limit behaviour $β=1−Δt/τ+O(Δt/τ)$ for $Δt/τ≪1$ corresponding to the forward differencing provided that $εR$ is constant. The coefficient $βexp$ is more suitable for larger ratios of $Δt/τ$⁠, in particular when τ is nearly constant especially. However, the former coefficient $βBE$ preserves symmetry of the equations (20) and (21) discretized in time and symmetry of the relaxation consequently. In the context of this work, we chose the coefficient $βBE$⁠, because we have not observed any significant differences for our physical settings and reliability was preferred over speed of computation. In any way, the method requires that $Δt≲τ$ holds. Therefore, the hydrodynamic time step is divided into smaller time steps proportional to τ and opacities are updated in each sub-step.

Following the terminology of [54], the Eq. (25) is rewritten using abstract operators:

where the abstract operators and function of the right hand side are defined on the computational domain V as:

This formulation unifies the equations inside the domain and at its boundary and allows further construction of mimetic operators. Newton boundary conditions in the form $qR⋅nV=ψ(ϵR−ϵout)$ are applied, where $nV$ is the outer normal of V and ψ and $ϵout$ are given functions at $∂V$⁠. The inference of the mimetic operators and full discretization of the system (27 and 28) follows the reference [54]. At this point, it must be noted that the presented modifications compared to the original paper do not violate the necessary properties of the operators for construction of the scheme. Namely, the operator Ω remains positive for $ψ>0$ and its discrete analogue has a diagonal matrix. Finally, an equation for the fluxes $qR$ is obtained and solved numerically. However, it should be stressed that the scheme satisfies conservation of $ϵR+ρcVTe$⁠, when the subsequent update of $Te$ is made by the expression:

that follows from summing of the discretized equations (20) and (21). Still, one may argue that the conservation of $εe$ is slightly violated, because of the temporal variation of $cV$ due to weak temperature dependency. This problem is addressed by the method known as Symmetric Semi-Implicit (SSI) [55]. Picard iterations are then performed until convergence of the non-linear system in $εR$ is reached. Even though the scheme may be less efficient than the solvers based on Newton–Krylov method for example [56,57], we find it highly reliable and robust for the reasons given in the preceding text and positive properties of the mimetic operators mentioned in [54].

In order to verify the non-locality of the radiation transport and evaluate it quantitatively, a very simplistic problem is proposed. It is inspired by the test problem well-known in non-local electron transport [58]. The initial periodic equilibrium temperature (⁠$TR=Te$⁠) perturbation $δT⁡cos(kz)$ from the base temperature $T0=1 keV$ decays by the action of a given radiation transport model. The asymptotic ratio $|qR|/|qRd|$ (with $qRd$ being the radiation diffusion flux) is then evaluated for $δT≪T0$⁠, where we chose $δT=5⋅10−2T0$⁠. The non-locality of the transport is set by an appropriate choice of the wave number k, while the the opacities are kept constant $κP≐100.8cm−1,κR≐13.9cm−1$ (corresponding to Copper with the density chosen as $ρ=1g/cm3$⁠). A uniform computational mesh on the domain $(zmin,zmax)$ with 5000 cells was constructed. The size of the computational domain was set proportionally to k as $|zmax−zmin|k=2π⋅50$ in order to retain constant number of the waves inside. The boundary conditions for the radiation diffusion were placed at local maximum and minimum of the waves, implying zero flux boundary condition there. In the case of the non-local transport, reflective boundary conditions were set due to the symmetry of the problem.

The results of the problem are plotted in Fig. 2, where the Planck opacity $κP$ is used as well as the average opacity κ introduced in Sec. 3.3. In addition, the results obtained by the radiation diffusion are shown. The dependency on the radiation Knudsen number $KnR∼k/κP$ reveals significant disagreement of the radiation diffusion and the non-local radiation transport model for the Planck opacity in the diffusive regime (⁠$KnR≪1$⁠), where the fluxes are underestimated by a factor of $≈7.26$⁠, that is essentially equal to the ratio of $κP/κR$⁠, because both models depend linearly on the inverse of the opacity. This result shows the necessity to use Rosseland opacity in this regime and justifies the construction of the average opacity, which then recovers precisely the diffusive limit. In the free-streaming limit (⁠$KnR≫1$⁠), the flux ratio declines showing inappropriateness of the diffusion approximation. The limit behaviour can be found in the form $|qR|/|qRd|∼(k/κP)−2$ as can be deduced from the linear radiation transport equation (18) by a Fourier analysis similar to [17].

The radiation diffusion model with various flux limiters described in Sec. 3.3.1 gives closely grouped curves with minor differences for each of the limiters. However, the limit behaviour for the diffusion model is identical in all cases, that is the dependency of the flux ratio $∼KnR−1$⁠, since the flux approaches the constant value $cεR$ in the free streaming limit $KnR≫1$ as described in Sec. 3.3.1. Unfortunately, this implies that the flux is dependent on the base temperature $T0$ in this test, that is independent of $δT$ for the perturbations. Consequently, the position of the resulting curves for all flux limiters is absolutely arbitrary in the sense that it depends on the ratio $δT/T0$ independently of $k/κP$⁠. This shows the non-feasibility of the flux limited diffusion model, where any flux limiter based only on local knowledge of the profiles (like on $KnR$ most typically) cannot correct the values of radiation fluxes. In principle, only knowledge of global profiles or global Fourier space is sufficient for description of fully non-local behaviour. In turn, this fact then justifies construction of a non-local transport model, where the flux ratios are absolutely independent of the $δT/T0$ ratio as we verified for the presented model.

The equation of state (EOS) data are taken from the SESAME tables [59,60] to properly close the equations of non-local radiation hydrodynamics (4–7). Discrete values of pressure and internal energy are bilinearly interpolated on the original grid provided by the SESAME tables. This approach is robust and guarantees that the interpolated data are consistent with the original data.

The momentum loss collision frequencies taken from [61] are used for calculation of the heat diffusion coefficients $κe$⁠, $κi$ and the energy exchange coefficients $G=ρ(∂εe/∂Te)ρνeiϵ$⁠, where $νeiϵ$ stands for the energy loss collision frequency. It is worth to note that the common problem of electron and ion heat exchange coefficients symmetry arises, because a general EOS does not provide the specific heats that could satisfy $(∂εe/∂Te)ρ/(∂εi/∂Ti)ρ=νieε/νeiε=Z¯$⁠, where the last equality relies on elasticity of the collisions and the quasi-neutrality condition $ne=Z¯ni$⁠, that is assumed to hold in the hydrodynamic description. We solved this problem by using ideal gas specific heats at this place, where the symmetry relation is clearly satisfied and a single heat exchange coefficient for electrons and ions is obtained.

The closure model for mean free paths of delocalized electrons is built on the collision frequencies taken from [61]. The formula for collision frequency has a form similar to (15) and holds with a good accuracy in hot plasma. However, the situation changes in colder pre-plasma, where the classical Coulomb logarithm $Λei=ln(bmax/bmin)$ approaches $Λei≈1$ as scattering on larger angles becomes more frequent. The usage of Λ defined in Sec. 3.2 gives at least physically relevant values of the electron mean free paths in cold material of the order of magnitude $∼1Å$ [62], but the results with non-local transport in this regime must be taken with a grain of salt and the improvement of the models in this area is left for a future work.

The present comparative study considers two materials: aluminium (⁠$Al$⁠, $Z=13$⁠, $ρ0=2.7g/cm3,T0=0.03 eV$⁠) and copper (⁠$Cu$⁠, $Z=29$⁠, $ρ0=8.96g/cm3,T0=0.03 eV$⁠), where Z is the proton number, $ρ0$ initial (solid) density and $T0$ initial temperature. The ablation of planar solid targets made of these materials is considered for laser intensities $1013W/cm2$⁠, $1014W/cm2$ and $1015W/cm2$⁠. Whereas the electron heat transport is non-local for all simulations, the transport of radiation is computed using the non-local as well as the non-equilibrium diffusive approach for comparison (see Sec. 3.3 for details).

The laser pulses are modelled as flat with total length of $10ns$ with a small Gaussian ramp with both, the temporal FWHM and offset of the intensity plateau, equal to $1.2 ns$⁠. The main reason for this configuration is that the pulses roughly approximate pre-pulses of high-intensity main pulses on laser facilities like ELI Beamlines (Extreme Light Infrastructure) for example [63–65]. The wavelength of the laser is taken as $λ0=0.35μm$ corresponding to the third harmonic frequency of Nd:glass laser systems, such as the L4 laser system at ELI Beamlines [66].

In order to attain maximal comparability of the numerical experiments, the identical computational mesh design was used in all cases. The simulation domain extends to the depth of $1500μm$ into the target. The mesh was composed of two parts, where the thin surface layer was only $10−2μm$ thick and contained $Nz=10$ computational cells geometrically spaced with geometrical factor $Cg=0.94$⁠:

Importance of this area is mainly for the initial phase of laser absorption, where the laser with intensity far beyond the ablation threshold of the given material interacts with the cold target and the surface material is very rapidly ablated. The evanescent wave penetrating the surface delocalizes the absorption provided that the mesh resolution is sufficiently high. The rest of the target has $Nz=1100$ computational cells spaced with geometrical factor $Cg=0.992$⁠.

The scan of intensities is presented in Fig. 3 for two materials, aluminium and copper, at the final time $10 ns$ in order to show the increasing importance of non-locality. The results for intensity $1014W/cm2$ were omitted for brevity, since they scale in a mostly predictable way between the two shown extremes. Spatial profiles of the electron density $ne$ normalized to the critical density $nc$ together with the temperatures of electrons $Te$⁠, ions $Ti$ and radiation $TR$ are presented. The Knudsen numbers of electrons $Kne=λe|∇Te|/Te$ and radiation defined as $KnR=1/κ|∇ϵR|/ϵR=4/κ|∇TR|/TR$ are also shown, where $λe$ is the electron mean free path. In addition, the Boltzmann number $Bo=52(pe+pi)|u|/(σTR4)$ is included. The corresponding energy fluxes including the hydrodynamic flux $qH=pu$ and the laser absorption rate, i.e. negative divergence of the Poynting vector (for details see Sec. 3.2), are plotted in the Fig. 4.

The plots are divided into several parts to distinguish different regimes of laser–target interaction. The leftmost area denoted by 0 corresponds to a cold solid target unperturbed by any shock waves or non-linear heat waves. The transition between the area 0 and I is strictly given by the position of the propagating shock wave into the target and I then comprises the shock wave plateau. The area II then starts at the head of the rarefaction wave and ends at the place where the laser gets absorbed. From the physical point of view, the conduction zone spans over this area and various energy transfer mechanisms are competing here, forming specific profiles of quantities in the process as discussed below. Finally, the area III extends over the laser absorption zone and is composed of mostly very rarefied coronal plasma expanding into vacuum.

The results show typical plasma profiles for laser–target interaction, that are composed of a laser absorption wave in the area III, an ablation wave in II and a shock wave in the area I. However, it can be recognized from the plot of $ne$ that the head of the laser absorption wave does not coincide with the position of $zc·(ne(zc)=nc)$⁠. This phenomenon is not normally observed in classical simulations without radiation transport (and we have also not observed this effect when the radiation transport was inactive). It can be identified as a consequence of strong radiation-matter coupling associated with the well-known structure called double ablation front (DAF) in the ICF community that forms in the plasma [67]. It basically arises from the balance of the hydrodynamic fluxes, radiation fluxes and electron heat fluxes. Following the terminology of the reference [67], the electron ablation wave forms in the vicinity of the maximal laser absorption point due to the collective effect of radiation and electron transport. However, this area is very small as can be noticed from the electron fluxes and roughly corresponds to $λe≤1μm$⁠. In contrast, photons can penetrate much deeper into the conduction zone, since their mean free path $λR=κ−1$ is of the order of millimetres. Due to high radiation fluxes in area II, the energy is effectively transported from the so called transition layer, where the radiation flux changes its sign, to the radiation ablation wave that removes material from the shock wave in zone I. The balance between the radiation, electron and hydrodynamic transport mechanisms leads to the creation of the so called density plateau between the fronts (which is not totally flat in logarithmic plot in fact, indicating continuous expansion of the matter). Matter and radiation are nearly equilibrated in this area as manifested by nearly identical electron and radiation temperatures. However, the transport is strongly dominated by radiation as can be deduced from the ratio of the fluxes or the Boltzmann number $Bo$⁠, that becomes $Bo≪1$ in this case. Steady state solutions can be derived for this problem using radiation and heat diffusion [67], giving an expression for the length of the plateau $∼Bo−1λR$⁠. Our results approximately agree with this scaling, since the Boltzmann number is about $∼10−1$ for the intensities $1015W/cm2$ and $λR∼1 cm$ at the right edge of II (as assumed in the reference), giving the scale length of a few millimetres. Hence, the ablation fronts are tightly coupled and evolution in densities affects both of them. Progressive enlarging of the shock wave plateau in I is followed by a natural relaxation of the ablative pressure and decrease of the radiation ablation head densities. Consequently, the electron ablation wave head density decreases gradually and crosses the critical density at a certain point. The larger separation that can be observed in the figures is then due to the very long intense laser pulses, where the DAF structure has evolved fully and shifted the maximum laser absorption region to the underdense plasma. A relatively sharp edge of the absorption region remains, resembling the classical scenario with absorption at the critical density, since the complex permittivity given by (9) strongly depends on the temperature and density, which exhibit a rapid change there. The laser absorption then affects these profiles again and the resulting profiles are given by a balance of the self-consistent laser absorption and the transport mechanisms. The radiation transport dominates at this place and it can be noticed that there exists a region, where the radiation flux $qR$ exceeds in absolute value the Poynting vector $qL$⁠. More importantly, the profiles of $qR$ and $qL$ are almost parallel in the vicinity of the maximal laser absorption point, indicating nearly identical divergences preventing the matter to heat up significantly. The plasma then evolves in the strong inverse bremsstrahlung regime following the terminology of [68], where the laser is absorbed before reaching the critical density. This assertion is also indicated by nearly full absorption in all simulations, that is characteristic for this regime. As qualitatively predicted by the models of self-regulating rarefaction wave, the separation increases with decreasing laser intensity [68,69]. The critical density $nc$ loses its predictive value for estimation of the maximum laser absorption point, which is then located further out in the underdense plasma. However, unchanged $nc$ is used in this paper for normative purposes in the figures, but the coronal plasma area III is distinguished correctly, independently of the value of $nc$⁠.

Important differences can be found between the radiation diffusion and the non-local radiation transport (see Sec. 3.3). The main effect can be seen in different length, slope and magnitude of the density and temperature profiles in the conduction zone (area II). This comes from the fact that the $KnR$ easily approaches values of the order of magnitude about unity. The diffusion approximation becomes absolutely insufficient in such cases and considerable differences with respect to the non-local transport were observed in Sec. 3.3.2. However, the true nature of non-locality resides in transport of energy on scales comparable with $λR$ in presence of steep plasma profile gradients. In this case, most of the radiation is produced in the vicinity of the critical plane [70]. The radiation then propagates deep into the conduction zone, where it is absorbed and partially re-radiated. The radiation mean-free-path $λR∼1mm$ agrees with the spatial extent of the area II. The diffusion approximation can barely capture this effect of non-locality, which increases continuously the ablative pressure and densities while lowering the temperatures in II. This leads to a slight increase of the Boltzmann number in case of diffusion, which is more than compensated by a decrease of the opacities giving longer double ablation front structure. Moreover, non-locality is also manifested by a slightly stronger separation of the radiation and electron temperatures there, because the strong radiation originating from the high temperature gradient is not fully equilibrated there. On the other side, the coronal plasma is totally transparent for the thermal radiation, that escapes then the plasma freely, because $λR$ is normally several times larger than the dimensions of the whole plasma corona. Consequently, $TR$ stays almost constant and independent of $Te$ and $Ti$⁠, but this effect is reproduced by the diffusive model too, since a non-equilibrium description is used (see Sec. 3.3.1). The electron temperatures are significantly higher in the simulations with the non-local radiation transport. This effect can be attributed to the higher radiation fluxes in II and III (i.e. radiative cooling) for the radiation diffusion, which can be considered overestimated for the given values of $KnR$⁠, since the medium is far from the diffusive limit (see Sec. 3.3.2).

Non-local transport of electron heat is less effective than radiation transport for both, aluminium and copper, in terms of magnitude of the energy fluxes. Its significance increases with increasing intensity and considerable fluxes can be seen only for aluminium target and the intensity of $1015W/cm2$ especially. However, the effect of non-local electron heat transport is significantly differs from the radiation transport. The electron Knudsen numbers are at least two orders of magnitude lower in area III and four orders in II. This indicates more localized transport that acts mainly in the vicinity of the critical plane, where it disperses the laser absorption area and forms the electron ablation wave. However, it must be stressed that most of the electron heat is transported by super-thermal electrons with the velocity $v≈3.7vTe$ in reality [71], which have significantly longer (⁠$≈200$⁠) mean free paths due to $v4$ dependency in (15). Since only single-group approximation is used, these spectral effects may be underestimated here, since these high-velocity electrons do not travel along whole trajectory during the transport. However, this non-linearity of the transport is taken into account by a local approximation [12]. The one-group electron mean free path is considerably shorter in II than in III and only a small pre-heating appears there compared to the dimensions of the area. However, electrons penetrate much deeper into the coronal plasma participating in the formation of the spatial electron temperature profile. The corona then undergoes non-isothermal expansion into the vacuum as indicated by the decreasing temperature in the right part of III. Unlike the heat diffusion approximation, where temperature saturation may occur, the flux inhibition effect of non-local models does not make this behaviour possible for the given conditions [72,73].

The intensity scaling has a mostly predictable behaviour, as the electron temperatures and the shock wave heights and velocities increase with the intensity. Gradually, the ion temperatures depart from the electron ones in the corona. It can be noticed from the ion energy equation (6) that there are no other global energy sources than the heat exchange term (except the mechanical work term $−pi∇⋅u$⁠, but it causes only energy losses in an expanding corona). Hence, the ratio of the ion to electron temperatures is lower for higher intensities, because the electron–ion collision frequency (15) decreases with electron temperature. The radiation temperatures are very low in all cases and only slightly exceed the temperature in II, because of high transparency of the corona due to high temperatures and low densities of the plasma.

Compared to Aluminium, Copper is a significantly better radiator and the simulations are completely dominated by radiation transport. Therefore, the coronal temperatures are decreased by a strong radiation cooling effect. This phenomenon partially explains considerably lower coronal temperatures in the simulations with copper target compared to the simulations with aluminium ones. The maximal electron temperatures are lower by a factor 2–3. Considering also higher ionization, the resulting plasma is more collisional (see Eq. (15)), heat transport is less effective in the corona and the electron Knudsen numbers are lower even for comparable spatial profiles. In contrast, Aluminium reaches an approximately twice as high shock wave velocity compared to Copper. In particular, they are $≈12 km/s$ and $≈59 km/s$ for Al and $≈6 km/s$ and $≈28 km/s$ for Cu for intensities $1013W/cm2$ and $1015W/cm2$⁠, respectively. On the other hand, the density plateau is longer in the case of Copper due to slightly lower Boltzmann numbers, which are proportional to the ratio between the hydrodynamic and radiation fluxes.

The importance of non-local radiation and electron transport was verified and several key aspects were identified in a series of simulations with different materials and laser intensities. In many cases, the radiation and electron heat diffusion approximation is not sufficient for realistic simulation of the considered effects.

One of the prominent features of the non-local transport models is self-consistent heat flux inhibition. Heat diffusion codes must rely on various heat flux limiters like the one defined by (1) or simply by redefinition of the heat fluxes $qe=min(−κe∇Te,feqf)$⁠. Similarly, the radiation diffusion models use radiation flux limiters as described in Sec. 3.3.1. This treatment is usually necessary at places, where the usage of diffusion approximation is not well justified, because radiation or electron Knudsen number is large (in particular $Kne≫0.06/Z$ [71]). The transport regime of radiation and electrons then cannot be considered diffusive and the application of diffusion models leads to overestimation of the calculated energy fluxes (see Sec. 3.3.2). However, proper energy flux limitation cannot be self-consistently determined from local values of $Te$ or $TR$⁠, since large gradients may lie in a range of a few mean free paths like the laser absorption wave for example (see Sec. 4.2). Therefore, the radiation fluxes differed considerably. As a result, formation of the double ablation front structure was affected, where significantly shorter and denser structures were observed in the case of the non-local treatment.

As expected for moderate-Z materials, radiation transport was weaker for the aluminium target and other energy transport mechanisms like electron heat transfer partially contributed as well. However, the situation was completely different for high-Z copper target, where radiation transport absolutely dominated over the other transport mechanisms. Considerable amount of absorbed laser power was converted to radiation and radiated away.

We have observed that radiation transport also evoked formation of a double ablation front structure for both materials, that led to the separation of the maximum laser absorption point from the critical plane due to moderate intense and long laser pulse. This interesting phenomenon requires deeper analysis of all the complex mechanisms involved and construction of a theoretical model for further explanation may be beneficial. This topic will be addressed in a future paper.

A limitation of the simulations is given by usage of scaling laws for opacities, that may not be sufficient for the description of all regimes and even the gray body approximation is highly simplifying. Also the single-group electron approximation may not be sufficient and the dynamics of the double ablation front can be significantly affected by these applied simplifications [74]. However, it is evident that proper non-local radiation transport plays a major role for hydrodynamic simulations with high intensities, especially with high-Z materials. In addition, non-local electron heat transport increases its importance for moderate-Z materials and non-local treatment of both transport mechanisms is necessary.

Finally, an accurate laser absorption method is absolutely substantial for the non-local transports, since it defines the gradient of temperature and density in the vicinity of the critical plane, where most of the non-local photons and electrons is created. A classical WKB is not sufficient for this purpose, since it is not able to describe this area realistically and produces a discontinuity in the Poynting vector at the critical plane essentially. Hence, a self-consistent laser absorption algorithm based on wave-optics was presented in Sec. 3.2 and several numerical methods were applied to increase its accuracy. However, a higher order laser absorption algorithm attaining a better numerical efficiency is a prospect of future work.

In addition, the topics of future work include improvements of the physical closure model for the heat and radiation transport models, that must provide the necessary coefficients for a wide range of input parameters. Division of energy spectrum of electrons and photons followed by a separate transport calculation (so called multi-group transport) is expected to also increase accuracy of the results, but at the cost of proportionally higher computational demands. In this text, only the non-local electron heat transport model was considered, but a fruitful insight may be also obtained from a comparison with the classical heat diffusion model.

From the global point of view, realistic and predictive simulations of laser–plasma interaction require at least the following extension:

• 2D axisymmetric geometry to handle basic interaction configurations,

• use of a more convenient collision operator rather than an approximative BGK approach in the non-local electron heat transport model,

• a better laser absorption model which goes beyond inverse Bremsstrahlung at the critical surface, i.e. inclusion of collective absorption processes driven by parametric instabilities in the plasma corona,

• accurate laser absorption method taking into account an oblique incidence angle and resonance absorption consequently,

• generation of hot electrons by the two preceding mechanisms and their transport in the plasma,

• self-consistent magnetic field effects which affect the non-local heat transport and which are not considered in the present analysis at all,

• the one-energy-group or gray body approximation should be replaced by a multi-group approach supplied by NLTE atomic physics models for closure coefficients.

It is clear that accurate radiation hydrodynamic simulations are mandatory for ICF to fully describe creation of the strong shock waves [75] and preheating caused by non-local electron transport [76], but also for future extreme high-field physics studies using multi-PW laser systems in order to characterize the plasma originating from the pre-pulse [77].

The authors declare that there is no conflicts of interest.

The authors would like to thank Vladimir Tikhonchuk for fruitful discussions.

The authors acknowledge support from the project High Field Initiative (HiFI) (CZ.02.1.01/0.0/0.0/15_003/0000449) and ELI Tools for Advanced Simulation (ELITAS) (CZ.02.1.01/0.0/0.0/16_013/0001793), both from European Regional Development Fund, the Czech Science Foundation project 18-20962S and Czech Technical University grant SGS16/247/OHK4/3T/14. This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement number 633053 (EUROfusion project CfP-AWP17-IFE-CEA-01).

Verification of the presented laser absorption algorithm (see Sec. 3.2) was performed by comparison of numerical results with an analytic solution for a simplified problem. In particular, the problem was designed as laser absorption on the density step function, where the left part was overcritical with the electron density $nel=2nc$ and the right part was underdense with the density $ner=(1−10−5)nc$ (⁠$nc$ is the critical density). The electron temperature has been chosen $Tel=1 eV$ on the left and $Ter=1 keV$ on the right. The material under test was fully ionized hydrogen and the computational domain spans from 0 to 40 microns with the interface placed at 10 microns. The computational mesh is uniform with the number of computational cells $Nz$ The laser wavelength was chosen as $λ0=1μm$ and only $νei$ defined in (15) was taken for simplicity. The reasons for this test scenario are two-fold, because the problem resembles a typical absorption in vicinity of the critical plane, where large gradients of the quantities appear and performance of the laser absorption solver must be validated. Moreover, the piece-wise homogeneous medium allows application of the formulas (10) and (11) for the analytic solution. Finally, it must be noted that the problem includes also the vacuum–plasma interface at the right boundary as it is typical for realistic configurations.

Analytic and numerical solutions were compared as shown in Fig. 5. In total, three procedures were used to obtain the analytic solution in order to provide a better overview of their inner working. The first one is the direct solution of the Helmholtz equation (8) decomposed into the one-sided wave equations for P and R. The left boundary condition is set to $R=0$⁠, because electric field of the reflected wave diverges on the left. On the opposite side, the incident wave is set $P=1+0i$ at the outer side of the vacuum interface. Application of the proper interface conditions for electric and magnetic field in addition to the boundary conditions and the exponential solution for P and R (see Sec. 3.2) yields a $8×8$ complex matrix, that must be inverted. On the other hand, the second analytic solution of the system (10) can be solved directly for V using the interface condition (12) and then for P in a similar fashion without involving any matrices at all. However, it is also worth noting that the solution in terms of the Poynting vector is independent of the phase provided at the right boundary condition. Lastly, the solution of the system (11) can be constructed analytically, but the electric field cannot be recovered anymore as discussed in Sec. 3.2. This system is then solved numerically. It has been practically verified that all the analytic solutions give an identical solution in terms of the resulting Poynting vector and the electric field eventually.

The results show that the discrete points of the numerical solution with $Nz=40$ lie almost exactly at the analytic curve unlike the rest of the configurations. This effect originates from the fact that the discontinuity precisely coincides with a computational node in this case. The error is then given only by the point projection of the solution. However, the discrete $L1$ error defined as $∑j|(qLa)j−(qLn)j|/∑j|(qLa)j|$ for discrete points of the analytic solution $(qLa)j$ and numerical $(qLn)j$ gives values on the level of machine precision and do not change significantly with increasing $Nz$ provided it still lies precisely at a computational node, i.e. $Nz$ is divisible by 4 (the ratio between size of the domain and position of the discontinuity). This effect is not repeated in the other configurations, where the $L1$ error is 0.31 for 38 computational cells, but decreases to 0.18 for 78 and 0.12 for 118 cells (⁠$Nz$ is chosen in the way that the remainder after division by 4 is constant). Expanding the series even further the first order convergence in $Nz$ is attained. This is an inevitable consequence of the introduced discontinuity and no algorithm can perform better in this test strictly speaking. This supports the effort to refine the computational mesh in vicinity of the critical plane as proposed in Sec. 3.2.1. However, the density profiles cross the value $nc$ smoothly in realistic configurations and this error is minimal. The piece-wise analytic solution then gives very accurate and smooth profiles of $∇⋅qL$ as observed in Sec. 4.2.

To summarize, the numerical solution agrees with all of the analytic solution reasonably. The piece-wise analytic scheme gives stable numerical solution almost independently of $Nz$ provided that the permittivity ε varies slowly in space. Under extreme conditions, where ε is discontinuous, the first order convergence is achieved at least. Besides that, the numerical method is able to recover the evanescent wave behind the critical plane, that is crucial for the initial phase of laser–target interaction.