The effects of self-generated magnetic fields on the ablative Rayleigh–Taylor (RT) instability are investigated in the linear regime. The main governing parameters are the Froude number (Fr), which stands for the ratio between ablative convection and acceleration of the target, and the Mach number at the ablation front (Ma), assumed to be small (isobaricity). During the development of the RT instability, magnetic fields are generated due to misalignment between pressure and density gradients (Biermann-battery effect). They accumulate at the section of the ablation front where the Nernst and the plasma velocities cancel each other. The magnetic field modifies the dynamics of the instability through the Righi–Leduc term, which acts as a heat source in the energy equation. It is found that the B fields affect perturbations with short wavelengths up to the most unstable wave in the spectrum. The B field plays a destabilizing role for moderate Froude numbers and becomes stabilizing for large Froude numbers. For plastic ablators, the Fr threshold is found to be $Fr=5$.

## I. INTRODUCTION

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\u22121$, 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.^{14} 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.^{20} 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.

## II. GOVERNING EQUATIONS

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\u2192$, ion plus electron (total) pressure *p*, plasma temperature *T*, and magnetic-field intensity $B\u2192$ is given by the continuity, total momentum, and total energy conservation equations together with Faraday's induction law,

where $m\xaf$ is the ion mass *m _{i}* divided by the plasma atomic number

*Z*, $E\u2192$ is the electric field, $g\u2192=ge\u2192x$ is the acceleration in the frame of reference of the ablation front, and $j\u2192=\u2212enu\u2192$ is the current, with $u\u2192$ 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\u2192$ given by low-frequency Ampère's law

We use the Braginskii^{20} expressions and notations for the ion–electron friction force $R\u2192$ and electron heat flux $q\u2192e$. We assume small electron Hall parameter $\omega e\tau e=(eB/mec)\tau e$, with $\tau e\u221dT3/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\u2192|$ are retained. The Hall and Ettingshausen terms are neglected since they are proportional to $|B\u2192|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 $\gamma 0nT\tau e/me=K\xafT5/2$, with $K\xaf$ being the Spitzer conduction constant.

### A. Base flow

We consider a one-dimensional steady base flow. The only component of the velocity is the streamwise component *u*_{0}, and we assume $B\u21920=0\u2192$. The governing equations for the base flow become

with $p0=(1+1/Z)n0T0$. As boundary conditions, we require that at $x\u2192\u2212\u221e$; we recover $n0=na$, $T0=Ta,\u2009u0=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\xafTa7/2/5paua$. Two dimensionless numbers arise: the (isothermal) Mach number at the ablation front $Ma=m\xafnaua2/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\u2192\u2212\u221e$. Typically, the ablation front is quasi-stationary; hence, we assume both $Ma\u226a1$ and $Ma2/Fr\u226a1$. Equation (15) can then be integrated yielding $p0=1+O(Ma2,Ma2/Fr)\u22481$, 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\u2192\u221e,\u2009T0\u223cx2/5$, $dT0/dx\u223cx\u22123/5$, and therefore, this model fails when $x\u223cMa\u22125$ and the non-isobaric effects need to be accounted for.

### B. Perturbed flow

We now linearize the equations to study the linear stability. We propose a modal decomposition in time and transversal coordinate as $exp\u2009(\gamma t+iky)$. Particularly, it is convenient to expand $p=p0(x)+Ma2p1(x)exp\u2009(\gamma t+iky)$. We denote as *v*_{1}, *b*_{1} 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=\u2212T1/T02$.

The magnetic field *b*_{1} has been normalized with a reference $Bref$ derived from the Biermann battery term in the induction equation (22), $Bref=cm\xafua/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=\gamma 0\alpha 0c2Ta2/10\pi e2pauaLa\u221dTa\u22123/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*is always greater than unity. The high atomic number case $Z\u2192\u221e$ 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.

_{N}. | Z = 1
. | Z = 2
. | Z = 3
. | Z = 4
. | $Z\u2192\u221e$ . |
---|---|---|---|---|---|

c _{R} | 5.75 | 3.45 | 2.81 | 2.51 | 1.71 |

c _{H} | 3.08 | 1.78 | 1.39 | 1.20 | 0.68 |

c _{N} | 1.28 | 1.31 | 1.36 | 1.41 | 1.83 |

. | Z = 1
. | Z = 2
. | Z = 3
. | Z = 4
. | $Z\u2192\u221e$ . |
---|---|---|---|---|---|

c _{R} | 5.75 | 3.45 | 2.81 | 2.51 | 1.71 |

c _{H} | 3.08 | 1.78 | 1.39 | 1.20 | 0.68 |

c _{N} | 1.28 | 1.31 | 1.36 | 1.41 | 1.83 |

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\u22121$, does not reach the sonic point; hence, we need $k\u226bMa5$. As shown in Ref. 8, for large Froude numbers, the cutoff takes place at $kc0\u223cFr\u22125/3$. Consequently, this analysis is valid for Froude numbers satisfying $Fr\u226aMa\u22123$. On the other hand, the analysis in Ref. 23 requires $Fr\u226b1$ for the SBM to be applicable. Hence, both analyses overlap in the range $1\u226aFr\u226aMa\u22123$.

## III. RESOLUTION METHOD

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\u2192={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\u2192\xb1\u221e$. Let ${q\u2192jc,j=1,7}$ and ${q\u2192jh,j=1,7}$ denote the modes in the cold shell, $x\u2192\u2212\u221e$, and in the hot plasma, $x\u2192\u221e$, 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*modes.

^{h}We integrate the system in Eq. (28) starting in the cold shell with a linear combination of the bounded modes: $q\u2192(x\u2192\u2212\u221e)=\u2211j=1j=jcajcq\u2192jc$. 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\u2192jc$ onto the unbounded hot mode $q\u2192lh$ 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$.

### A. Modal analysis in the cold shell

In the cold shell, we have $T0=1$ and $dT0/dx\u21920$. The matrix *M* in Eq. (28) becomes constant and it is diagonalizable. The modes evolve exponentially as $q\u2192jc=q\u2192j0c\u2009exp\u2009(\lambda 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 ${\lambda jc}$ can be grouped as ${\lambda hydroc,\lambda magc}$, with

and

The first hydrodynamic mode corresponds to vorticity flow and the second and third correspond to incompressible flow motion, $\u22072p=0$, whereas the last two modes are thermal modes. Only the second and fourth modes are bounded. It can be seen that, when $k\u226a1$, the thermal mode decays faster than the incompressible mode, $\lambda 4c\u223c1$ compared to $\lambda 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.

### B. Modal analysis in the hot plasma

In the hot plasma, we have $T0\u2248(5x/2)2/5$. The point $x\u2192\u221e$ 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,^{24} 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\u2192=S\xb7Y\u2192$, with $Y\u2192={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\u2192={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 $\eta =2T0(x)5/2/5$ so that the system in Eq. (28) can be written as

where

The matrix $M\xaf$ can be expanded as $M\xaf=\u2211i=0i=\u221eM\xafi\eta 1\u2212i/10$, with $M\xafi$ being constant matrices. According to Wasow, the seven independent solutions (modes) ${Y\u2192j,j=1,7}$ of such a system can be sought in the form

with $\lambda j,i$ eigenvalues and $Y\u2192j,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 $\eta \u2212i/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*_{7} stands for the identity matrix in seven dimensions. It can be seen that the leading order $\lambda j,0$ and $Y\u2192j,0$ correspond to the eigenvalues and eigenvectors of $M\xaf0$. At this point, the procedure differs depending on the algebraic multiplicity of each leading order eigenvalue.

If the multiplicity of a given eigenvalue $\lambda j,0$ is one, then left multiplying Eq. (37) by $Y\u2192j,0T$ enables obtaining $\lambda j,i$ recursively. At each “*i*” step, the eigenvector $Y\u2192j,i$ can be derived from Eq. (37) once $\lambda j,i$ is known. Now $Y\u2192j,i$ is not fully determined. According to Wasow, we have the freedom to choose the projection of $Y\u2192j,i$, *i* > 0 into $Y\u2192j,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 $\lambda j,0$ has the algebraic multiplicity two. The associated eigenvectors ${Y\u2192j1,0,Y\u2192j2,0}$ span a two-dimensional space. Let ${\mu 1,\mu 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 $\lambda 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 $\lambda j,1$ is also a double root, then ${\mu 1,\mu 2}$ cannot be determined and the aforementioned procedure must be repeated successively in higher orders until we reach a $i=i*$ for which $\lambda 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\u22125)$, 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 $(\lambda 6,14)$ can be simplified based on the order of magnitude estimations. The temporal growth rate is of the order of $\gamma \u223ck/Fr$. It is well known that the cutoff wavenumber *k _{c}* decreases with the Froude number. For large Froude numbers, we have $kc\u223cFr\u22125/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\gamma \u2273Ma$ 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.

Mode . | Number (j) . | First nonzero eigenvalue . |
---|---|---|

$Vorticity$ | 1 | $\lambda 1,14=\u2212(25)2/5\gamma $ |

$Incompressible+Thermal\u2009(I)$ | 2 | $\lambda 2,10=\u2212k,\u2003\lambda 2,13=\u2212MakcR27/1053/101\u2212cN$ |

$Incompressible+Thermal\u2009(II)$ | 3 | $\lambda 3,10=\u2212k,\u2003\lambda 3,13=MakcR27/1053/101\u2212cN$ |

$Incompressible+Thermal\u2009(III)$ | 4 | $\lambda 4,10=k,\u2003\lambda 4,13=\u2212MakcR27/1053/10cN\u22121$ |

$Incompressible+Thermal\u2009(IV)$ | 5 | $\lambda 5,10=k,\u2003\lambda 5,13=MakcR27/1053/10cN\u22121$ |

$Induction$ | 6 | $\lambda 6,14=\u2212(25)2/5Fr\gamma \u2212Ma2cRFr(1\u2212cN)$ |

$B\u2009diffusion$ | 7 | $\lambda 7,0=4Rem(1\u2212cN)$ |

Mode . | Number (j) . | First nonzero eigenvalue . |
---|---|---|

$Vorticity$ | 1 | $\lambda 1,14=\u2212(25)2/5\gamma $ |

$Incompressible+Thermal\u2009(I)$ | 2 | $\lambda 2,10=\u2212k,\u2003\lambda 2,13=\u2212MakcR27/1053/101\u2212cN$ |

$Incompressible+Thermal\u2009(II)$ | 3 | $\lambda 3,10=\u2212k,\u2003\lambda 3,13=MakcR27/1053/101\u2212cN$ |

$Incompressible+Thermal\u2009(III)$ | 4 | $\lambda 4,10=k,\u2003\lambda 4,13=\u2212MakcR27/1053/10cN\u22121$ |

$Incompressible+Thermal\u2009(IV)$ | 5 | $\lambda 5,10=k,\u2003\lambda 5,13=MakcR27/1053/10cN\u22121$ |

$Induction$ | 6 | $\lambda 6,14=\u2212(25)2/5Fr\gamma \u2212Ma2cRFr(1\u2212cN)$ |

$B\u2009diffusion$ | 7 | $\lambda 7,0=4Rem(1\u2212cN)$ |

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

## IV. MAGNETIC-FIELD GENERATION AND STRUCTURE

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*_{1} at $x\u2192\xb1\u221e$. 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\u2261cN/(cN\u22121)$. 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\u226b1$, 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\u22121)+Td(x\u2212xv)$, with $Td=(cN\u22121)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*. For a highly conductive plasma, the structure of the magnetic field close to

_{B}*x*derived from induction equation yields

_{v}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\u2192\xb1\u221e$ toward $x\u2192xv\xb1,$ respectively. It can be seen that there is a critical value of the growth rate below which the magnetic field diverges, $\gamma cr=(cN\u22121)5/2/cN5/2$. For high-

*Z*ablators, we have $\gamma cr\u22480.14$, while for hydrogen plasmas, we have a lower threshold $\gamma cr\u22480.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, $\gamma <\gamma 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^{27} to derive the structure of such layer. We strain the spatial coordinate as $\eta =(x\u2212xv)/\u03f5$, with $\u03f5\u226a1$. 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 $\eta \u2192\u2212\eta $, the solution is decomposed as a linear combination of symmetric and skew-symmetric functions, for example, $\phi =Cinsym\phi sym+Cinskw\phi skw$, with $\phi sym=1,\u2009d\phi skw/d\eta =1$ at *η* = 0 and $Cinsym,\u2009Cinskw$ constants. Expanding each function for $\eta \u2192\u221e$ gives

with $\Gamma (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\alpha /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 $\gamma =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\u2192\u221e$ cases. In the latter, the magnetic field diverges at the section where the bulk and Nernst velocities cancel out, $T0(xv)=cN/(cn\u22121)$, 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.

## V. LINEAR STABILITY ANALYSIS WITH B FIELD

In this section, we consider the feedback of the self-generated magnetic field into the Rayleigh–Taylor instability. We recall that, in the linear regime, the B field only appears as a heat source proportional to the square of the Mach number in the energy Eq. (21).

The dispersion relation is computed in Fig. 3 for large atomic numbers, $Z\u2192\u221e$, and different Froude numbers. The results for the hydrodynamics coupled to the magnetic field are compared to the uncoupled case $(Ma=0)$. The case without Nernst is shown for comparison. A perfectly conductive plasma, $Rem\u2192\u221e$, is assumed, and we have chosen a small Mach number such that $cRMa2=10\u22122$. For small wave numbers, the dispersion relation is similar to the RT instability of immiscible fluids with Atwood number equal to unity, $\gamma =kg$. Ablation stabilization becomes effective for larger wave numbers until the RT instability is suppressed at a certain cutoff *k _{c}*. The cutoff wavelength decreases at larger Froude numbers.

The self-generated B field significantly modifies the dispersion relation even for small Mach numbers. Its effect becomes important when ablation stabilization takes place (near the maximum of the spectrum). It is not monotonic and can either enhance or stabilize the RT instability depending on the Froude number. For $Fr=1$, the self-generated B field plays a destabilizing role and enlarges the region of unstable modes. This effect is greatly amplified by the Nernst term, notably close to the cutoff. For $Fr=10$, however, the stabilizing character of the B field switches when the Nernst term is considered, shortening the range of unstable modes *k*. This transition in the stabilizing character is completed in the case $Fr=100$, where the self-generated B field shortens the range of unstable modes in both with and without Nernst cases. In the former case, this effect is specially pronounced since the cutoff is reduced by almost 40%.

The RT instability is more sensitive to the effect of the self-generated B fields when the Froude number is large. In this case, the cutoff takes place at a wavelength larger than the ablation-front scale length. The perturbations penetrate deeper into the hot plasma, and the Righi–Leduc term, which is proportional to $T05dT0/dx\u223cT07/2$, becomes more important in the energy equation (21). The sensitivity of the RT instability to the self-generated B field for large Fr is in agreement with the SBM analysis in Ref. 23, where the B field is stabilizing and the Righi–Leduc term rescaled with the variables in the coronal plasma becomes of order unity.

The dispersion relation for smaller atomic numbers is shown in Fig. 4. Only the case with Nernst is considered, and the Righi–Leduc coefficient is kept constant to $cRMa2=10\u22122$. The black curves $(Z\u2192\u221e)$ are the same as in Fig. 3. It can be seen that the stabilizing effect of the B field is enhanced for lower atomic numbers. According to this, a hydrogen ablation front is more sensitive to the self-generated B field than plastic ablators.

An interesting feature observed is that, when the Nernst term is considered, the stabilizing character of the B field is maintained for a given Froude number. Therefore, there exists a threshold Fr depending exclusively upon *Z* for which the stabilizing character of the B field switches. This is shown in Fig. 5, where the difference in growth ratio is shown for the most unstable mode for every Froude number. The dependence of threshold Fr on *Z* is rather weak, since we obtained $Fr=6,5,3.4$ for $Z=1,4,\u221e$, respectively. Although this figure is shown for $cRMa2=10\u22122$, the same threshold Fr is obtained when varying the Mach number.

Thus far, only a perfectly conductive plasma was considered. The effect of magnetic diffusivity is shown in Table III, where the coefficient computed measures the similarity of the growth rate to the pure hydrodynamic case (close to zero) or to the perfectly conductive case (close to one). It can be seen that, for $Fr=1$, the system is very sensitive to the magnetic Reynolds number. Ablation fronts in ICF tend to be very diffusive, and $Rem$ is usually low. Therefore, the effect of the B field on perturbations whose wavelength is comparable to the ablation front length scale is suppressed by resistivity. Contrary to this behavior, long-wavelength perturbations are insensitive to $Rem$, as can be seen for the case $Fr=10$. At this Froude number, unstable perturbations penetrate deep into the hot conductive plasma, and the effect of the B field (which, for these perturbations is stabilizing) is unaltered by diffusion.

. | $Rem=10\u22122$ . | $Rem=10\u22121$ . | $Rem=1$ . | $Rem=10$ . |
---|---|---|---|---|

$Fr=1,k=10\u22121$ | 0.076 | 0.28 | 0.62 | 0.97 |

$Fr=10,k=5\xd710\u22123$ | 1.16 | 1.19 | 1.18 | 1.13 |

. | $Rem=10\u22122$ . | $Rem=10\u22121$ . | $Rem=1$ . | $Rem=10$ . |
---|---|---|---|---|

$Fr=1,k=10\u22121$ | 0.076 | 0.28 | 0.62 | 0.97 |

$Fr=10,k=5\xd710\u22123$ | 1.16 | 1.19 | 1.18 | 1.13 |

The transition in the stabilizing character of the magnetic field lies in the change of sign of the Biermann battery term (B. B.), which is depicted in Fig. 6. Two cases are shown: one where the effect of the self-generated B field is destabilizing, which corresponds to the most unstable wave for $Fr=1$, and the other where it is stabilizing, which corresponds to the most unstable wave for $Fr=100$. In all cases, the Mach number is assumed to be zero (induction equation decoupled) although very similar profiles are obtained for different $Ma$. In the destabilizing case, the Biermann battery term remains positive for moderate values of *kx*. However, in the stabilizing case, it changes sign in the hottest region. Since the Biermann battery is the source term in the induction equation, it has a direct effect on the sign of the magnetic field, which in turn affects the direction in which the heat-flux lines are bent. This ultimately dictates the stabilizing character of the B field. In the section considered in Fig. 6 (right side of the spike), negative B. B. regions, which create positive B-field regions, converge the heat-flux lines toward the spike, enhancing ablation. Note that it is this positive part of the B-field profile the one observed in Fig. 5 in Ref. 23, where the perturbed profiles in the conduction layer are shown for $kLa=211Ma5\u226a1$, for which the B field is stabilizing. Finally, the Nernst term favors the transition toward stabilization as it convects the B field from the hotter regions of the corona where it is positive toward the ablation front.

## VI. CONCLUSION

The effects of self-generated magnetic fields on the ablative Rayleigh–Taylor instability are investigated using a semi-analytical approach. The analysis assumes isobaricity (small Mach number) and is restricted to the linear regime. The solution is sought as a composition of single modes (eigenfunctions). The main governing parameter is the Froude number, which stands for the ratio between ablative convection and acceleration.

The mathematical procedure used to numerically determine the eigenfunctions of a homogeneous system of equations with boundary conditions at infinity is reviewed. The key point of the analysis corresponds to the treatment of modes at the irregular singular points of the system, which arise in ablation-front structures.

The structure of the self-generated B field during the development of the RT instability is derived. The region where the Nernst convection and the bulk plasma velocity cancels out contains a point of B field accumulation. The B-field can develop a singularity in this section for a perfectly conductive plasma and growth rates smaller than a certain threshold, which depends on the atomic number. This singularity is smoothed out when diffusion is included.

Finally, the effect of B fields on the development of the RT instability is analyzed. In the linear regime, the coupling takes place through the Righi–Leduc term, which acts as a heat source in the energy equation. Although this term is proportional to the Mach number squared, its effect on the RT dynamics is significant even for small Ma. It can be interpreted as a deflection of the heat-flux lines with respect to the temperature gradient direction. This effect is found to be not monotonic but depends only on the Froude number. For moderate Froude numbers, the Righi–Leduc term enhances the RT instability by diverging the heat-flux lines and, therefore, reducing ablation. For large Froude numbers, however, the B field becomes stabilizing. The Righi–Leduc converges the heat-flux lines toward the spike. As a result, mass ablation is increased and the range of unstable wave numbers is reduced. The threshold between the stabilizing and destabilizing behavior of the B field occurs at a critical Froude number that depends weakly upon the atomic number of the plasma *Z*. For large *Z*, it takes place at $Fr=3.4$, while for hydrogen ablation fronts, it occurs at $Fr=6$. Magnetic diffusivity at the ablation-front scale is found to reduce the effect of B fields for moderate Froude numbers. Unstable perturbations in the large Fr regime are insensitive to magnetic diffusion since they penetrate deep into the hot and highly conductive corona.

## ACKNOWLEDGMENTS

This work was supported by the Department of Energy Office of Science, Fusion Energy Sciences program Grant Nos. DE-SC0016258 and DE-SC0014318. J. Sanz was also supported by the Spanish Ministerio de Economía y Competitividad, Project No. RTI2018-098801-B-I00. H. Aluie was also supported by U.S. DOE Grant Nos. DE-SC0020229 and DE-SC0019329, U.S. NASA Grant No. 80NSSC18K0772, and U.S. NNSA Grant Nos. DE-NA0003856 and DE-NA0003914. This material is based upon work supported by the Department of Energy National Nuclear Security Administration under Award No. DE-NA0003856, the University of Rochester, and the New York State Energy Research and Development Authority. This report was prepared as an account of work sponsored by an agency of the U.S. Government. Neither the U.S. Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the U.S. Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the U.S. Government or any agency thereof.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: MATRIX OF THE LINEARIZED EQUATIONS

The matrix *M* in the system of Eq. (28) is

with

and

The prime denotes derivative with respect to *x*.

### APPENDIX B: IRREGULAR SINGULAR POINT ANALYSIS

In this appendix, we summarize the procedure described in Wasow^{24} to derive the asymptotic expansion of the solution ordinary differential equations in the neighborhood of an irregular singular point. More precisely, we seek to derive the asymptotic expansion at $x\u2192\u221e$ of the solution of the system

with $q\u22650$ integer and $A(x)$ square matrix holomorphic at $x\u2192\u221e$. We assume that the matrix $A(x)$ admits a power-series expansion $A(x)=\u2211r=0r=\u221eArx\u2212r$. We obtain as many independent solutions (modes) as the order of $A(x)$.

Before explaining the procedure, we recall some definitions and theorems. We define the shifting matrix *H _{n}* as a null

*n*-square matrix with ones on the first upper diagonal. We introduced the direct sum ⊕ of given matrices

*A*,

*B*as

Any constant square matrix can, therefore, be reduced by similarity transformations to direct sums of diagonal and Jordan blocks. The Jordan blocks arise whenever the geometric multiplicity (g. m.) of an eigenvalue is lower than the algebraic multiplicity (a. m.). For example, any given matrix *A* with eigenvalues *λ*_{1} with a. m. = g. m. = 1 and *λ*_{2} with a. m. = 2, g. m. = 1 can be reduced by similarity to

We recall Theorem 4.1. in Wasow. Let *M* and *N* be square matrices with dimensions *m* × *m* and *n* × *n*, respectively. Let *X* be a *m* × *n* unknown matrix and *Q* a *m* × *n* matrix. Then, the system of equations $M\xb7X\u2212X\xb7N=Q$ has a unique solution *X* if *M* and *N* do not share any eigenvalue. Finally, we rewrite the Lemma 19.1 in Wasow. Let *H* and *K* be shifting matrices of orders *h* and *k*, respectively. Let *Q* be a *h* × *k* matrix whose first *h* − 1 rows are known constants, while the entries $\alpha 1,\alpha 2,\u2026,\alpha k$ of the last row are unknown. Then, the numbers $\alpha 1,\alpha 2,\u2026,\alpha k$ can be uniquely determined in such a way that the equation $H\xb7X\u2212X\xb7K=Q$ for the *h* × *k* matrix *X* is soluble.

We proceed now to derive the asymptotic expansion at $x\u2192\u221e$. We introduce the following change of variables:

such that Eq. (B1) transforms into

with

The aim of this transformation is to make $B(x)$ as simple as possible. More precisely, we want it to be composed of direct sum of matrices *B ^{ii}*, which allows Eq. (B4) to be split into problems of lower order. Ideally, the aim is to perform successive transformations such that $B(x)$ becomes a diagonal matrix. The problem would, therefore, be reduced to a set of independent equations of order one. Using the nomenclature in Sec. III B, the resolution of each problem would yield the eigenvalues $\lambda j,r$, while the eigenvectors $Y\u2192j,r$ are related to the matrices

*P*.

_{r}We choose $P0=I$, which yields $B0=A0$. The higher-order terms are given by the following recursion formula:

with

with the last term being absent for $r\u2212q\u22121<0$. Let us assume that *A*_{0} is already formed by direct sum of diagonal and Jordan blocks. Next, we group these blocks into bigger blocks with no eigenvalues in common between them. Without loss of generality, we write

We partition accordingly

and can propose

Introducing this into Eq. (B6) gives the following set of equations:

Since $A011$ and $A022$ have no eigenvalue in common, then by virtue of Theorem 4.1 in Wasow, this system of equations has a solution and we can obtain $Br11,\u2009Pr12,\u2009Pr21$, $Br22$. Applying this procedure allows one to decompose Eq. (B4) into independent subsystems of a lower order for each different eigenvalue of *A*_{0}, that is, every block $A0ii$. Notice that our particular proposition for *P _{r}* corresponds to requiring that the eigenvectors $Y\u2192j,r$,

*r*> 0 be perpendicular to $Y\u2192j,0$ in Sec. III B.

Now, depending on the algebraic and geometric multiplicity of every eigenvalue, we can have the following:

- a. m. = g. m. $=1$. The corresponding block has order one and can be solved directly; that is, we have(B9)$x\u2212qdzdx=(b0+b1t+\cdots )z1,$whose solution is(B10)$z=exp\u2009[\u222bxcx(b0+b1s+\cdots )sqds]$
with

*x*and arbitrary constant._{c} - a. m. = g. m. > 1. The corresponding block has order larger than one and is diagonal $A0ii=\lambda iI$. Performing the change of variables $Z=exp\u2009[\lambda 1xq+1/(q+1)]W$ gives$x\u2212(q\u22121)dWdx=(\u2211r=0r=\u221eBr+1iixr)\xb7W.$
The leading term is now $B1ii$ and can be further decomposed into subsystems of lower dimension by applying the change of variables in Eq. (B3) again.

- a. m. > g. m. The corresponding block $A0ii$ cannot be diagonalized. The transformation $Z=exp\u2009[\lambda 1xq+1/(q+1)]W$ transforms the matrix into direct sum of shifting matrices $A0ii=H1\u2295H2\u2295\cdots \u2295Hs$, where the order of one of them is two at least (Jordan block). In this case, we need to perform a shearing transformation to transform the Jordan block and make the matrix diagonalizable. But before this, it is necessary to precondition the system as follows. We partition the $Brii$,
*r*> 0 in the same structure as $A0ii$. Now introducing another change of variables $W=R(x)\xb7V$, and applying Lemma 19.1 in Wasow, we can simplify $Brii$,*r*> 0 in such a way that the only nonzero entries occur in the rows corresponding to the last row of the blocks*H*, $k=1,2,\u2026,s$. Once the system is preconditioned as described, we perform the shearing transformation. Let_{k}*S*be a diagonal matrix such that $S(x)=Diag(1,x\u2212g,x\u22122g,\u2026,x\u2212(n\u22121)g)$, with*n*the dimension of $A0ii$, and let $V=S(x)\xb7U$. Introducing this, we can choose*g*such that one of the off-diagonal terms becomes as important as the entries in $A0ii$. Now defining*p*as the smallest integer such that*pg*is an integer, and introducing $x=\alpha \eta p,\u2009\alpha =p1/(g\u2212q\u22121)$, we arrive to(B11)$\eta \u2212hdUd\eta =B(\eta )\xb7U,$with $h=p(q+1\u2212g)\u22121$ and $B(\eta )=\u2211r=1r=\u221eBr\eta \u2212r=S\u22121\xb7Bii\xb7S\u2212x\u2212qS\u22121\xb7dS/dx$. The aforementioned precondition ensures that either $B0$ has more than one eigenvalue or it presents shifting matrices of a higher dimension than the shifting matrices of $A0ii$. If the latter case occurs, then performing successive shearing transformations will eventually yield a $B0$ that is composed of a single shifting matrix of order

*n*. Then, another final shearing transformation will yield a $B0$ with different eigenvalues. When this is the case, the corresponding system in Eq. (B11) can be further partitioned into smaller subsystems by the methodology described, until every block gets a. m. = g. m. = 1 and every mode can be independently obtained.

The procedure described here makes it possible to obtain the asymptotic expansion of each mode in a systematic way. However, in real cases, it can get rather troublesome due to the number of transformations required. For this reason, it is important to try “home-made” transformations based on a direct observation of the system and experience that can allow one to simplify it in a more straightforward way, or at least precondition it in such a way that this method can be applied with fewer number of transformations required. This was the case in Sec. III B, where the shearing matrix introduced *ad hoc* preconditions the system in such a way that the method can be straightforwardly applied, which is tantamount to propose solutions such as Eq. (36).

Finally, we illustrate the systematic procedure with the following example:

The leading term *A*_{0} can be expressed as a direct sum of shifting matrices $A0=H1\u2009\u2295\u2009H2$. A preconditioning needs to be done since the following term in the expansion, *A*_{1}, has nonzero entries in the second row. We, therefore, propose a change of variables $Y=P(x)\xb7Y2$ with the choice

motivated by the Lemma 19.1 in Wasow. This transforms the system in Eq. (B12) into

We now perform a shearing transformation, $Y2=S\xb7Y3$ with the choice *g* = 2/3. Introducing accordingly the change of variable $\eta =3x1/3$ leads to

A direct change of basis $Y3=T\xb7Y4$, with $T={{0,0,\u22121/2},{1,0,0},{0,1,0}}$, rearranges the previous system as

The leading order of the matrix is a shifting matrix of a higher order *H*_{3}. Notice that straightforwardly applying the shearing transformation to Eq. (B12) would not increase the order of the shifting matrices of the resulting system and, hence, the benefit of the precondition. A second shearing transformation $Y4=S\xb7Y5$ applied to Eq. (B16) with *g* = 1 yields finally

This system is fully diagonalizable, and the modes are related to the different eigenvalues of the matrix ${\lambda j}={8,(7+35)/2,(7\u221235)/2}$. The most unstable mode of the system in Eq. (B12) evolves then as $\u223cx8/3\u22122/3=x6/3$.