Modeling hydrodynamic instabilities of double ablation fronts in inertial confinement fusion

A linear Rayleigh-Taylor instability theory of double ablation (DA) fronts is developed for directdrive inertial confinement fusion. Two approaches are discussed: an analytical discontinuity model for the radiation dominated regime of very steep DA front structure, and a numerical self-consistent model that covers more general hydrodynamic profiles behaviours. Dispersion relation results are compared to 2D simulations.


INTRODUCTION
The use of moderate-Z ablators in inertial confinement fusion (ICF) targets is motivated by the reduction in target preheat that was observed in the glass (SiO 2 ) compared to plastic ablators in recent work [1].Fuel preheat, caused by hot-electrons generated by two-plasmon-decay (TPD) instability, increases the shell entropy preventing the thermonuclear pellet from achieving high compression.Therefore, target preheat, as one source of compression degradation, is a limiting factor to ICF goal that should be mitigated.
Hydrodynamic profiles of the ablation region in the acceleration phase of implosion considerably change when moderate-Z ablators are used.In particular, a second ablation front appears forming a double ablation (DA) front structure [2][3][4].This DA front structure (see figure 1) is driven by two energy transport mechanisms, the electronic heat flux (q e ) and the radiative energy flux (S r ).The flow is dominated by radiation at the peak density.Therefore, inner front is called radiative ablation (RA) front.A plateau in density/ temperature develops between the two fronts.In this region, radiation and matter are almost in equilibrium, radiation energy is diffusively transported and the radiation energy flux is dominant.Around the outer ablation front, a transition layer (TL) develops where S r changes its sign.Outer front is always driven by electronic heat flux.However, the developed TL is as well a very strong emitter of radiation.Thus, this second ablation front is called electronic-radiative ablation (ERA) front [3,5].Beyond the ERA front, almost all the incoming electron heat flux is outwardly radiated.The ablation region is known to be Rayleigh-Taylor (RT) unstable during the pellet implosion.This hydrodynamic instability is critical for obtaining appropriate implosions, and it is indeed one limiting factor that prevents the achievement of ignition.
Physics involved in the stability of DA fronts are not well known yet, and only recently has been object of study [4,5].Current linear models [6,7] do not take into account both energy transport mechanisms that lead to the emergence of the DA front, but only a general thermal conductivity transport.Therefore, they cannot reproduce a DA front structure as a background flow, and are expected a e-mail: yanez@celia.u-bordeaux1.frThis is an Open Access article distributed under the terms of the Creative Commons Attribution License 2.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.to fail in the determination of the linear RT instability growth rate of imploding shells with moderate-Z ablators.The goal of this research is, then, to develop a linear stability theory for double ablation fronts in the ICF context.

ANALYTICAL DISCONTINUITY MODEL FOR STABILITY ANALYSIS
Discontinuity models have been extensively used in the linear stability study of ablation fronts.The main characteristic of this type of models is to consider the ablation front as a zero-thickness interface with steady flow at both sides of it.Solutions are reached by imposing conservation laws across the surface discontinuity.Nevertheless, the problem is not mathematically closed.A possible option is to keep one free parameter that can be taken a posteriori as an input from 1D simulations [8].
In the present model, the existence of two ablation fronts leads to keep two free parameters: r D , the ratio between the plateau and the peak density and d p , the plateau length.Two other parameters are needed to complete the picture: the Froude number, Fr t ≡ v 2 t /gL st and the material dependent parameter t ≡ 16( T 4 t K P t )( t T t )/(5P t v t ) 2 .Here, g is the slab acceleration, P , v, and T are the plasma pressure, velocity and temperature, is the Stefan-Boltzmann constant, K P , the Planck opacity, t , the Spitzer thermal conductivity and L st ≡ 2 t T t /(5P t v t ) is the Spitzer length.Subscripts t refer to the conditions at the plateau density/temperature.Dimensionless parameter t comes naturally in the dimensionless ERA front energy equation [5].Actually, the 1D theory of double ablation fronts [3] states that a value t > 1 is required for the development of the ERA front.This parameter strongly depends on the ablator material, and simulations [3,4] provide an average value of 20 for the SiO 2 .
We consider two steady ablation fronts of zero-thickness that are displaced by a small wavelike perturbation of i exp( t + ikx) (with i = 1 or 2 for the RA or ERA front, respectively).Then, the isobaric fluid equations are linearized under a subsonic flow assumption.For each homogeneous region we reach a characteristic equation that gives the possible perturbation modes (∝ e y ) developing inside.These modes must be bounded at y = ±∞ since boundary conditions state that perturbations are localized in the ablation fronts.Ahead of the RA front, there is only one bounded mode (incompressible mode = k).Inside the plateau region (a finite region), five modes can develop: two incompressible modes ( i = ±k), one vorticity mode ( v = − ) and two thermal conduction modes ) is a characteristic length of the radiative transport.The study of the region behind the ERA front (the optically thin region) is more complex and a full detailed description can be found in [5].Additionally to conservation laws, we need two other equations to close the problem.They come from assuming that both ablation front surfaces behave as an isotherm (T 1 + T 0 = 0).Thus, we arrive to a linear system of equations that leads to a lengthy analytical dispersion relation formula (expression can be found in the Appendix F of [5]).

04008-p.2 IFSA 2011
Some dispersion curves given by the analytical formula can be seen in figure 1 (right).If we look to the left side of that figure, one can see that the dispersion relations for small wavenumbers follow the trend toward the classical growth rate value of RT instability as expected ( ≈ cl = √ kg).Before reaching the maximum growth rate max , there is the intermediate region of the dispersion curve which is rather dependent on the plateau length.A sufficient separation between both fronts (d p = 30L St or longer in the case of the figures) makes the dispersion curve suit the growth rate of the equivalent single RA front (inner one), which is the most unstable of the two fronts.This feature may indicate a weak influence of couple modes.Nevertheless, the situation changes when the plateau length decreases, strengthening an additional stabilization due to the interaction between both fronts (couple modes).It is noticeable the apparition of a local minimum that gives a characteristic double-hump shape to the dispersion curve (see curve (a) of figure 1 (right)).Finally, it is worth noting that the cutoff wavenumber is almost equal to the equivalent single RA front one and it weakly depends on d p .

NUMERICAL SELF-CONSISTENT MODEL FOR STABILITY ANALYSIS
Discontinuity models are valid when the wavelength of the perturbation is much longer than the characteristic length of the ablation fronts (L RA and L ERA ).Otherwise, the model fails since the effect of a smoothly varying density (L is finite) is not considered.In the classical RT instability, the unstable response of long-wavelength modes (kL 1) is different from short-wavenlength modes (kL 1) that are localized within the smooth interface.Results from both limits are reproduced in the asymptotic formula = √ kg/(1 + kL).In ICF target, the ablation process complicates the evolution of the instability.Thus, for smooth density profiles, instead of using discontinuity models, we should apply self-consistent models, which take into account the hydro-profiles of the background flow.Self-consistent models [6,7] start from studying the temperature and density profiles in the ablation region, for, next, imposing over those profiles linear perturbations, which are analyzed, numerically or analytically, in order to complete the stability study.
In the present model, we use a simplification of an existing 1D hydrodynamic radiation theory [3] as the background flow.It assumes two differentiated radiative regimes of the matter occurring, respectively, in regions covering either the outer (optically thin) and the inner (optically thick) ablation fronts.Equations that provide the hydrodynamic profiles are: where P 0 is a constant pressure across the double ablation structure, R = ¯ R T ν and = ¯ T 5/2 are the Rosseland and the Spitzer thermal conductivities, and K P = KP T −11/2 is the Planck mean opacity for a fully ionized plasma.We consider then a wave-like perturbation and linearize the equations (1) in their first order perturbation.This leads to an eigenvalue problem expressed as V 1 = AV 1 , where V 1 is the vector of perturbed hydrodynamic quantities normalized by their base flow values at the plateau region, A is a matrix depending on the profile parameters, Fr t , the wavenumber k and the growth rate , and the prime denotes derivative with respect to 0 = T 0 /T t .Modal analysis of the matrix system in the limits y = ±∞ gives the boundary conditions needed to integrate it.Actually, only bounded eigenmodes are physically valid, since perturbations can only grow in the RT unstable region (ablation fronts) and must decay to zero away from it.Therefore, boundary conditions are made up as a linear combination of the bounded modes of each limit.There are two positive (bounded) eigenvalues at y = −∞: 1/2 /2.We can then write the left boundary condition as and a 2 are still two undefined parameters.Boundary condition at y = +∞ is also a linear combination of two bounded eigenmodes ( < 0) and their expression can be found in [5].The computational method of the dispersion relation can be summarized as follows: we integrate the matrix system from the boundary condition V − 0 .Normalizing the solution with the most unbounded mode, it will tend to a constant vector c when 0 1.Each component of c is a linear combination of a 1 and a 2 , that is to say c i = a 1 f i (1, 0) + a 2 f i (0, 1), and shall be null to ensure a bounded solution.In order to have a non-trivial solution, we select any two components of c and force the determinant to be zero (f i (1, 0)f j (0, 1) − f j (1, 0)f i (0, 1) = 0), that leads to the dispersion relation curve given the parameters t , ν, Fr t , and r D .
In figure 2 we show the example of two dispersion curves computed with this numerical self-consistent approach (L finite) compared to the analytical discontinuity model results (L 1).Background flow in both examples only differs in the power index of the radiative conductivity ν.This results in a steeper RA front in the case of larger ν.In the limit ν 1, we would have L RA 1 and, consequently, the discontinuity interface assumption for the RA front would be valid.It is noticeable that the discontinuity model approximates rather well the results of the self-consistent model for the steeper profile, specially in points as the cutoff wavenumber, but, on the other hand, we can see a considerable disagreement with the smoother profile, where, for instance, the cutoff wavenumber is overstimated (roughly the double).Therefore, in those ICF targets that present smooth hydro-profiles it is convenient to estimate the RT growth rate with a self-consistent model.
In order to validate the theory, 2D simulations have been carried out with the code CHIC [9].These simulations take a 20 m thick SiO 2 layer irradiated by a laser pulse with a maximum intensity of about 04008-p.4 IFSA 2011 200 TW/cm 2 .Results are taken in the time interval 1.5-2.5 ns at the RA front, and compared to the dispersion relations given by both the discontinuity and the self consistent models.This comparison is shown in figure 3, and a remarkable agreement with the self-consistent model can be seen.

Figure 1 .
Figure 1.Left: diagram of the perturbed DA front structure.Right: dispersion relation given by the discontinuity model for t = 20, Fr t = 1, r D = 0.25 and three different plateau lengths.(a) d p = 2, (b) d p = 5 and (c) d p = 30.Thin dotdashed line corresponds to the classical growth rate = √ kg.

1 04008-p. 3 EPJFigure 2 .Figure 3 .
Figure 2. Left: density profiles for t = 20, r D = 0.35 and d p ≈ 30 and (a) ν = 10 and = 182 and (b) ν = 5 and = 83.Right: dispersion relation for these density profiles and Fr t = 1.Solid line corresponds to the selfconsistent model, dashed line to the discontinuity model and dotdashed line to the classical growth rate = √ kg.