Magnetic-field generation and its effect on ablative Rayleigh–Taylor instability in diffusive ablation fronts

In direct-drive inertial confinement fusion (ICF),^1,2 a high-power laser illuminates a spherical shell filled with deuterium–tritium (DT) fuel. The laser energy is convected toward the cold dense shell by heat conduction. Mass is ablated off the shell and expands, forming a structure commonly known as ablation front. During the acceleration phase, the Rayleigh–Taylor (RT) instability^3–7 grows at the ablation front degrading the integrity of the imploding shell.

In the linear regime, mass ablation stabilizes perturbations with wavelengths shorter than a certain cutoff, $kco−1$⁠, which depends on the Froude number Fr. This number is defined as the ratio between ablative convection and acceleration of the capsule; the larger the Froude number, the larger the cutoff wavelength. There are different mechanisms through which ablation damps or even suppresses the RT instability. When the Froude number is large, the main stabilizing mechanism corresponds to the “rocket effect,” by which a restoring pressure is generated on the crest of the corrugated ablation front (spikes).^8–11 In the opposite limit, the “fire-polishing effect,” by which the crest of the corrugated ablation front evaporates faster than the valleys, and thermal smoothing become the predominant mechanisms.^12,13 The rocket effect corresponds to a conservative stabilizing force, while the fire-polishing effect has a nonconservative nature (damping).

During the development of the RT instability, magnetic (B) fields are generated due to the misalignment of gradients of density and pressure, known as the baroclinic or Biermann battery effect.¹⁴ Self-generated B fields have been investigated both in numerical studies^15,16 and in experiments,^17–19 with the latter reporting measurements of the order of several Megagauss. The presence of magnetic fields modifies the transport coefficients in a plasma. Thermal conduction perpendicular to the B field is reduced, and a heat-flux component along isotherms is generated. The latter effect is known as Righi–Leduc and is important in the linear regime since it is of the first order in B.²⁰ In essence, the Righi–Leduc effect deflects the heat-flux lines, which in turn has a direct effect on the stabilization of the ablative RT instability. Depending on the orientation of the self-generated B field, mass ablation (and, therefore, stabilization) can be either enhanced or reduced, as shown in Fig. 1.

In Ref. 21, the RT instability with self-generated magnetic fields was simulated in a configuration relevant to laser fusion. It was reported that significant magnetization of the plasma could be attained even in the linear regime. In addition, a destabilizing feedback between induction and hydrodynamics was observed, which resulted in B-field generation enhancement. This behavior was not obtained for every simulation setup studied. In Ref. 22, the effect of the self-generated magnetic fields in the stagnation dynamics of ICF capsules was assessed. In the particular simulation setup of that paper, the Righi–Leduc term deflected the heat flux around the spikes, pumping heat toward the bubble, corresponding to the destabilizing configuration in Fig. 1. However, in the authors' words, the relative importance of each magnetized transport effect was expected to vary with the drive profile and perturbation combination.

More recently, a self-consistent analytical study of the effects of self-generated B fields on the Rayleigh–Taylor and Darrieus–Landau (DL) instabilities has been published in Ref. 23. The analysis assumed large Fr and focused on perturbation wavelengths longer than the ablation-front length scale. A sharp boundary model (SBM) of the ablation front and the adjacent coronal plasma led to a dispersion relation encompassing both instabilities. It was found that the B field plays a stabilizing role for the RT-unstable perturbations, contrary to the results reported in Ref. 21, while the DL-unstable ones are unaffected by them.

In this paper, we perform a linear stability analysis of the Rayleigh–Taylor instability with self-generated magnetic fields in an ablation-front length scale. We retain its diffusive structure and assume that the laser energy absorption takes place sufficiently far (no DL-unstable waves are considered). We first derive the structure of the self-generated B field and then analyze its feedback on the hydrodynamics. We find that the B field enhances the RT instability for moderate Fr and contributes to its stabilization for large Fr, consistent with Ref. 23.

This paper is organized as follows. In Sec. II, the governing equations are presented and linearized. In Sec. III, we describe the method to obtain the dispersion relation. In Sec. IV, we discuss the generation and structure of the magnetic field in a purely hydrodynamic RT instability. In Sec. V, the dispersion relation with self-generated magnetic fields is computed. Finally, in Sec. VI, conclusions are drawn.

We consider a semi-infinite plasma slab expanding due to the deposition of laser energy in planar geometry. We choose a hydrodynamic description of a fully ionized, single-species plasma and we assume quasi-neutrality and the same ion and electron temperatures. With these assumptions, the evolution of the electron density n, ion velocity $v→$⁠, ion plus electron (total) pressure p, plasma temperature T, and magnetic-field intensity $B→$ is given by the continuity, total momentum, and total energy conservation equations together with Faraday's induction law,

where $m¯$ is the ion mass m_i divided by the plasma atomic number Z, $E→$ is the electric field, $g→=ge→x$ is the acceleration in the frame of reference of the ablation front, and $j→=−enu→$ is the current, with $u→$ being the difference between electron and ion velocities. We have neglected the ion heat flux in Eq. (3) and electron inertia, and we use the ideal gas equation of state

The electric field and the current are given by the electron momentum equation, leading to a generalized Ohm's law,

with the current $j→$ given by low-frequency Ampère's law

We use the Braginskii²⁰ expressions and notations for the ion–electron friction force $R→$ and electron heat flux $q→e$⁠. We assume small electron Hall parameter $ωeτe=(eB/mec)τe$⁠, with $τe∝T3/2/n$ being the electron collision time,

In the linear stability analysis performed in this paper, only terms up to the first order in $|B→|$ are retained. The Hall and Ettingshausen terms are neglected since they are proportional to $|B→|2$⁠, and all the coefficients are taken in the unmagnetized limit. In addition, the thermoelectric effects cancel out when added up in Eq. (3), and the ion–electron friction becomes of the second order in the same equation. We choose to rewrite $γ0nTτe/me=K¯T5/2$⁠, with $K¯$ being the Spitzer conduction constant.

By means of Eqs. (4), (6), and (8), we can build an induction equation for the magnetic field

Equations (1)–(3) and (10) become a system of four equations for the unknowns n, $v→$⁠, T, and $B→$⁠. We decompose every variable into base plus perturbation, and we assume small perturbations. The base variables are denoted with the subindex “0” and the perturbations with “1.”

We consider a one-dimensional steady base flow. The only component of the velocity is the streamwise component u₀, and we assume $B→0=0→$⁠. The governing equations for the base flow become

with $p0=(1+1/Z)n0T0$⁠. As boundary conditions, we require that at $x→−∞$⁠; we recover $n0=na$⁠, $T0=Ta, u0=ua$⁠, where the subindex a denotes the state of the plasma upstream of the ablation front (cold shell).

The variables are normalized with the values of reference in the cold shell. A scale length can be formed based on thermal conduction at the ablation front, $La=2K¯Ta7/2/5paua$⁠. Two dimensionless numbers arise: the (isothermal) Mach number at the ablation front $Ma=m¯naua2/pa$ and the Froude number $Fr=ua2/gLa$⁠. The normalized continuity, momentum, and energy equations read

This system of equations is completed with the normalized boundary conditions $T0=u0=1$ at $x→−∞$⁠. Typically, the ablation front is quasi-stationary; hence, we assume both $Ma≪1$ and $Ma2/Fr≪1$⁠. Equation (15) can then be integrated yielding $p0=1+O(Ma2,Ma2/Fr)≈1$⁠, which stands for isobaricity and implies $u0=T0$⁠. Energy Eq. (16) can be integrated once yielding

Since the problem lacks a spatial reference, we choose to solve this equation by imposing that the origin x = 0 corresponds to the point of maximum temperature gradient; hence, $T0(x=0)=7/5$⁠. It is interesting to notice that when $x→∞, T0∼x2/5$⁠, $dT0/dx∼x−3/5$⁠, and therefore, this model fails when $x∼Ma−5$ and the non-isobaric effects need to be accounted for.

We now linearize the equations to study the linear stability. We propose a modal decomposition in time and transversal coordinate as $exp (γt+iky)$⁠. Particularly, it is convenient to expand $p=p0(x)+Ma2p1(x)exp (γt+iky)$⁠. We denote as v₁, b₁ the transversal velocity and magnetic field, respectively, and we absorb the imaginary unit within them. Normalizing with the variables in the cold shell $(na,Ta,ua,La)$⁠, we find

where the density perturbation has been related to temperature through the linearized equation of state in the quasi-isobaric limit, $n1=−T1/T02$⁠.

The magnetic field b₁ has been normalized with a reference $Bref$ derived from the Biermann battery term in the induction equation (22), $Bref=cm¯ua/eLa$⁠. In the right-hand side of this equation, we have the magnetic-field diffusion inversely proportional to the magnetic Reynolds number $Rem=uaLa/Dm$⁠, with D_m being the magnetic diffusivity given by $Dm=γ0α0c2Ta2/10πe2pauaLa∝Ta−3/2$⁠. The numerical expression for the Mach and magnetic Reynolds numbers is

where m_p is the proton mass. The scale length of the ablation front is

with Λ standing for the Coulomb logarithm. The Froude number becomes then

For practical applications for experimental designs, we recall that, according to the SBM in Ref. 23, the laser intensity I can be related to the ablation pressure and velocity as $I=25paua/32Ma2$⁠.

The three numerical coefficients appearing in the Righi–Leduc term $(cR)$⁠, current enthalpy convection $(cH)$⁠, and Nernst term $(cN)$ only depend on the atomic number

The numerical coefficients as given by Braginskii are shown in Table I. It can be seen that c_R decreases with the atomic number and that c_N is always greater than unity. The high atomic number case $Z→∞$ shall be understood as Z large enough that the Braginskii values can be used while the energy flow is still dominated by electronic heat conduction.

It is important to highlight that, to the first order, the effect of the magnetic field on the hydrodynamics is restricted to the last term in parentheses in Eq. (21), which acts as a heat source. Two effects are accounted for here: the first one is the Righi–Leduc term and the second one comes from the difference between electron and ion enthalpy convection (which in the following is called current enthalpy convection). They are proportional to $Ma2$⁠; therefore, they should be neglected within the framework of isobaricity. However, we choose to retain them and we prove numerically that the effect of the B field in the RT instability turns out to be significant even for small Mach numbers.

We define the range of applications of this study, and we compare it to the analysis in Ref. 23. Isobaricity is maintained as long as the perturbation wavelength, $k−1$⁠, does not reach the sonic point; hence, we need $k≫Ma5$⁠. As shown in Ref. 8, for large Froude numbers, the cutoff takes place at $kc0∼Fr−5/3$⁠. Consequently, this analysis is valid for Froude numbers satisfying $Fr≪Ma−3$⁠. On the other hand, the analysis in Ref. 23 requires $Fr≫1$ for the SBM to be applicable. Hence, both analyses overlap in the range $1≪Fr≪Ma−3$⁠.

We use the same method as in Ref. 25 to derive the dispersion relation; that is, for every k > 0, we compute the value of γ that allows for the existence of nontrivial solutions of the homogeneous system of governing Eqs. (18)–(22). Since we are mainly interested in unstable waves, we restrict the analysis to positive values of γ. The system of equations can be rewritten in a matrix form as

with $q→={u1,p1,v1,T1,dT1/dx,b1,db1/dx}T$ being the state variable vector. The matrix M is given in  Appendix A. We require that nontrivial solutions be bounded at $x→±∞$⁠. Let ${q→jc,j=1,7}$ and ${q→jh,j=1,7}$ denote the modes in the cold shell, $x→−∞$⁠, and in the hot plasma, $x→∞$⁠, respectively. We arrange them in such a way that the bounded modes in the cold shell correspond to the first j^c modes, while the unbounded modes in the hot plasma correspond to the first j^h modes.

We integrate the system in Eq. (28) starting in the cold shell with a linear combination of the bounded modes: $q→(x→−∞)=∑j=1j=jcajcq→jc$⁠. When integrating in x, the state variable vector will project into each mode in the hot plasma. Let $alh$ be the component of such projection in the unbounded mode l in the hot plasma. It linearly depends on ${ajc,j=1,jc}$ through

The function f_jl represents the projection of the integration of $q→jc$ onto the unbounded hot mode $q→lh$ and must be computed numerically. The dispersion relation D is obtained by requiring the projections $alh$ be zero for a nontrivial set of ${ajc,j=1,jc}$⁠, which results in

The dispersion relation will, therefore, be discrete if the number of bounded cold modes is equal to unbounded hot modes, that is, as long as

We now proceed to derive the modes in the cold shell and hot plasma. This is a crucial step to appropriately compute the functions f_jl. However, it is not essential to understand the main results of this paper; therefore, the reader not interested may skip to Sec. IV. The only important result is that the condition Eq. (31) is satisfied with $jc=jh=3$⁠.

In the cold shell, we have $T0=1$ and $dT0/dx→0$⁠. The matrix M in Eq. (28) becomes constant and it is diagonalizable. The modes evolve exponentially as $q→jc=q→j0c exp (λjcx)$⁠. The induction equation becomes uncoupled from the hydrodynamics; therefore, the modes can be split into five hydrodynamic modes and two magnetic modes. The eigenvalues ${λjc}$ can be grouped as ${λhydroc,λmagc}$⁠, with

and

The first hydrodynamic mode corresponds to vorticity flow and the second and third correspond to incompressible flow motion, $∇2p=0$⁠, whereas the last two modes are thermal modes. Only the second and fourth modes are bounded. It can be seen that, when $k≪1$⁠, the thermal mode decays faster than the incompressible mode, $λ4c∼1$ compared to $λ2c=k$⁠. The two magnetic modes are of a diffusive nature, and only one of them is bounded. In summary, there are three bounded modes in the cold plasma.

In the hot plasma, we have $T0≈(5x/2)2/5$⁠. The point $x→∞$ is an irregular singular point of the system Eq. (28). In order to obtain the modes in this region, we must follow the procedure thoroughly described in Wasow,²⁴ for which the leading order of the matrix M does not suffice to characterize the boundedness of the modes and we must take into account higher orders. This procedure is summarized in  Appendix B. The key point is to find a “shearing transformation” that enables full diagonalization of the matrix M by transforming the Jordan blocks that may appear in higher-order terms in the asymptotic expansion of M. At this point, the procedure greatly differs when considering or not the coupling between induction and the hydrodynamics. If it is not considered, a shearing transformation for our particular problem corresponds to introducing $q→=S·Y→$⁠, with $Y→={u1,p1,v1,T05/2T1,T05/2dT1/dx}T$ and S the corresponding diagonal shearing matrix. A similar transformation was used in Ref. 25 also in the context of ablative RT instability. However, when the coupling is considered, a rather different transformation needs to be performed: $Y→={u1,p1,v1,T03/4T1,T03/4dT1/dx,T07/4b1,T07/4db1/dx}T$⁠. In the following, only the latter case is treated.

It is convenient to introduce the variable $η=2T0(x)5/2/5$ so that the system in Eq. (28) can be written as

where

The matrix $M¯$ can be expanded as $M¯=∑i=0i=∞M¯iη1−i/10$⁠, with $M¯i$ being constant matrices. According to Wasow, the seven independent solutions (modes) ${Y→j,j=1,7}$ of such a system can be sought in the form

with $λj,i$ eigenvalues and $Y→j,i$ eigenvectors. The reason for not considering higher-order terms i > 20 in the first summation lies in the fact that these terms would be exponentials of an decreasing function of η. Accordingly, they allow for a series expansion in powers of $η−i/10$⁠, being redundant with the second summation term in Eq. (36). Inserting this solution into Eq. (34) and gathering terms with the same power in η yield a set of system of equations for each mode j,

The matrix I₇ stands for the identity matrix in seven dimensions. It can be seen that the leading order $λj,0$ and $Y→j,0$ correspond to the eigenvalues and eigenvectors of $M¯0$⁠. At this point, the procedure differs depending on the algebraic multiplicity of each leading order eigenvalue.

If the multiplicity of a given eigenvalue $λj,0$ is one, then left multiplying Eq. (37) by $Y→j,0T$ enables obtaining $λj,i$ recursively. At each “i” step, the eigenvector $Y→j,i$ can be derived from Eq. (37) once $λj,i$ is known. Now $Y→j,i$ is not fully determined. According to Wasow, we have the freedom to choose the projection of $Y→j,i$⁠, i > 0 into $Y→j,0$⁠, and we choose them to be perpendicular.

The procedure is slightly more complicated if the multiplicity of the eigenvalue is larger than one. Without loss of generality, let us assume that $λj,0$ has the algebraic multiplicity two. The associated eigenvectors ${Y→j1,0,Y→j2,0}$ span a two-dimensional space. Let ${μ1,μ2}$ be two arbitrary multipliers, taking the order i = 1 in Eq. (37), we can form

This system yields a nontrivial solution if

This condition is used to obtain $λj,1$⁠. Generally, two different values are obtained; hence, this double mode splits into two different modes in an order higher than the leading order (in this case, it would split at order i = 1). From this order on, the two modes can be treated as explained for modes with a multiplicity equal to one. However, if $λj,1$ is also a double root, then ${μ1,μ2}$ cannot be determined and the aforementioned procedure must be repeated successively in higher orders until we reach a $i=i*$ for which $λj,i*$ has two different values.

The eigenvalues obtained are summarized in Table II. They correspond to vorticity $(j=1)$⁠, combinations of incompressible plus thermal modes $(j=2−5)$⁠, induction $(j=6)$⁠, and B diffusion $(j=7)$⁠. Since the highest-order eigenvalues are rather lengthy expressions, only the first nonzero eigenvalues are shown as they characterize the boundedness of the mode. We can see that vorticity, incompressible + thermal (I) and (II) are bounded modes. The boundedness of the induction mode depends on the Nernst term. The numerator of this mode $(λ6,14)$ can be simplified based on the order of magnitude estimations. The temporal growth rate is of the order of $γ∼k/Fr$⁠. It is well known that the cutoff wavenumber k_c decreases with the Froude number. For large Froude numbers, we have $kc∼Fr−5/3$⁠. On the other hand, to avoid breaking the isobaricity assumption, we restrict this study to modes that do not penetrate deep into the hot plasma, $k>Ma5$⁠. Therefore, we can estimate $Frγ≳Ma$ and neglect the squared Mach number terms. We can then simplify this eigenvalue to

which is unbounded when the Nernst term is included and bounded when it is not. The analysis is exactly the opposite for the B-diffusion mode, which is bounded with Nernst and unbounded without. Therefore, we always have three unbounded modes in the hot plasma.

Finally, we write the leading-order eigenvectors. For vorticity, we find

while with the combined incompressible plus thermal modes, we have

For induction and B-diffusion, we have

We first analyze the generation of the magnetic field assuming no coupling with the hydrodynamics (strict limit $Ma=0$⁠). We assume that the hydrodynamic profiles are smooth and decay sufficiently fast at infinity. The magnetic-field structure is given by the induction equation (22). It must be integrated imposing the boundedness of b₁ at $x→±∞$⁠. The magnetic field is generated by the Biermann battery term and undergoes diffusion and convection. The latter is given by both the bulk plasma and the Nernst velocity. In the particular configuration treated in this paper, the Nernst velocity is opposed to the plasma velocity and convects the B field toward the cold plasma.

A crucial point is the position where the plasma velocity equals the Nernst velocity, which takes place at $x=xv$ such that $T0(xv)=Tv≡cN/(cN−1)$⁠. Both to the left and to the right of this surface, the B field is convected toward it and accumulates. For a highly conductive plasma, $Rem≫1$⁠, the B field becomes singular under certain conditions.

To study this point more in detail, we expand the base temperature around $x=xv$ as $T0=cN/(cN−1)+Td(x−xv)$⁠, with $Td=(cN−1)3/2/cN5/2$⁠. The Biermann battery term is approximately constant in the neighborhood of x_v. It is given by the hydrodynamic problem and we set it equal to S_B. For a highly conductive plasma, the structure of the magnetic field close to x_v derived from induction equation yields

with

the super-index “±” referring to the left “–” and right “+” of x_v and $Cout$ being constants. They are given by the integration of the induction equation from $x→±∞$ toward $x→xv±,$ respectively. It can be seen that there is a critical value of the growth rate below which the magnetic field diverges, $γcr=(cN−1)5/2/cN5/2$⁠. For high-Z ablators, we have $γcr≈0.14$⁠, while for hydrogen plasmas, we have a lower threshold $γcr≈0.022$⁠. The physical reason behind the existence of a growth rate threshold lies in the fact that the temporal growth rate is given by the Biermann battery term (hydrodynamics). When it is slow, $γ<γcr$⁠, the B-field accumulation due to convection is stronger than the temporal rate at which the B field can increase, leading to a singularity. This is a consequence of the asymptotic nature (in time) of the solutions we are seeking (modes). In this case, strong gradients take place at $x=xv$ and a diffusive boundary layer arises.

We apply the singular perturbation theory²⁷ to derive the structure of such layer. We strain the spatial coordinate as $η=(x−xv)/ϵ$⁠, with $ϵ≪1$⁠. Introducing

the induction equation, retaining diffusion, becomes

It can be inferred that the width of the boundary layer is

Equation (42) then becomes

Since this equation is invariant to the change $η→−η$⁠, the solution is decomposed as a linear combination of symmetric and skew-symmetric functions, for example, $φ=Cinsymφsym+Cinskwφskw$⁠, with $φsym=1, dφskw/dη=1$ at η = 0 and $Cinsym, Cinskw$ constants. Expanding each function for $η→∞$ gives

with $Γ(z)$ being the Euler Gamma function. Matching with the outer solution (39) gives

The magnetic field b_1, therefore, exhibits a peak of the order of $O(Remα/2)$ for large magnetic Reynolds number.

Self-generated magnetic-field profiles have been computed in Fig. 2. The Biermann battery term is obtained from the hydrodynamic profiles in the most unstable mode for $Fr=1$⁠, which corresponds to k = 0.037 and $γ=0.085$ (see Fig. 3). The profiles have been normalized with the maximum value of the Biermann battery term, which takes place close to the ablation-front position, x = 0. In this mode, the magnetic-field profile greatly differs between the Z = 1 and $Z→∞$ cases. In the latter, the magnetic field diverges at the section where the bulk and Nernst velocities cancel out, $T0(xv)=cN/(cn−1)$⁠, assuming a perfectly conductive plasma. This singularity may suggest the development of a magnetic shield screening the ablation front. However, this is not the case in typical ICF implosions, as the ablation front is very diffusive (low magnetic Reynolds number). When diffusivity is taken into account, Fig. 2(b), the singularity is smoothed out and the magnetic-field peak decreases for lower $Rem$⁠. The profile without Nernst convection is also depicted for comparison. In this case, the maximum of the B field is substantially lower since there is no accumulation, but it decays more slowly and reaches hotter regions of the corona. The features described are in agreement with the analysis of Nernst thermomagnetic waves in high-beta flows performed in Ref. 26.