New self-similar radiation-hydrodynamics solutions in the high-energy density, equilibrium diffusion limit

This work presents semi-analytic solutions to a radiation-hydrodynamics problem of a radiation source driving an initially cold medium. Our solutions are in the equilibrium diffusion limit, include material motion and allow for radiation-dominated situations where the radiation energy is comparable to (or greater than) the material internal energy density. As such, this work is a generalization of the classical Marshak wave problem that assumes no material motion and that the radiation energy is negligible. Including radiation energy density in the model serves to slow down the wave propagation. The solutions provide insight into the impact of radiation energy and material motion, as well as present a novel verification test for radiation transport packages. As a verification test, the solution exercises the radiation–matter coupling terms and their v/c treatment without needing a hydrodynamics solve. An example comparison between the self-similar solution and a numerical code is given. Tables of the self-similar solutions are also provided.


Introduction
This study deals with a new, self-similar solution to the equations of radiation hydrodynamics (RHD) in the high-energy density regime. In particular it is an extension of the classic Marshak wave problem [1], as solved in detail by Petschek [2]. This problem has been reprised in several monographs that cover RHD [3][4][5]. The numerical solution of Marshak wave problems was covered in detail by Nelson and Reynolds [6]. It is also worth pointing out that the theory of admissible self-similar solutions is covered in some detail in Coggeshall and Axford [7].
The Marshak wave is a special, soluble problem in RHD, and can arise in high-energydensity experiments where radiation strikes a cold material (e.g. when radiation emitted from a hohlraum strikes a fusion target). As we will demonstrate below, the classical version of this problem is valid in the regime where the radiation energy flux is high, but the radiation energy density is negligible in regard to the material internal energy. Also, the material is considered to be stationary, that is, radiation energy impinges on a quiescent material and drives a radiation wave that travels faster than the speed of sound in the material. Recently, there has been a renewal of interest in RHD solutions. Lowrie and co-workers [8,9] have published several studies of radiating shock waves using different radiation models, and there has also been developments in the theory of such shocks in different regimes [10][11][12].
Part of the motivation for these recent papers on RHD behavior has been the necessity of verifying numerical simulation codes for RHD. Verification in this sense means demonstrating that the code is solving the intended equations and that numerical errors behave as expected (e.g. going to zero at the correct rate as the mesh is refined). Our solutions below, besides their inherent interest as novel solutions, can be used to this end in the verification of the radiation-matter coupling part of a simulation code. Our solutions can be used in code verification, regardless of the radiation transport model employed in the code (e.g. flux-limited diffusion, discrete ordinates, Monte Carlo, etc) because our solution is given in the equilibrium diffusion limit, an asymptotic limit of the RHD equations that is possessed by nearly all approximate transport models [13][14][15][16].
In the following section we review our RHD model and then derive the equilibrium diffusion model. In section 4 we pose the problem of radiation impinging on a cold slab and present an approach to estimating the location of the wavefront. Solution profiles are evaluated in section 5, and we demonstrate the usefulness of these solutions for code verification in section 6.

Radiation hydrodynamics model
We begin with the equations for RHD as given by Lowrie and Morel previously [17]. The hydrodynamics equations, in non-dimensional form under the condition of local thermodynamic equilibrium [5] are ∂ρ ∂t where ρ is the material density, v is the relative velocity and p is the pressure. E, the total specific energy is given by where v = | v| and e is the internal specific energy with a simple equation of state e = c v T . The specific heat, c v , is a uniform and constant parameter. The non-dimensionalization we use aŝ where hatted variables are dimensional, l is a characteristic length scale, ρ ∞ is a characteristic density and T ∞ is a reference temperature. The reference sound speed, a ∞ , is given using the constant, reference specific heat c v∞ as Also, in our model we have included the non-dimensional parameters where c is the speed of light and a r = 4σ SB /c is the radiation constant with σ SB the Stefan-Boltzmann constant. The interpretation of these parameters is as follows: C is a measure of how relativistic the flow is and P is a measure of how much energy is in the radiation field compared with the material internal energy. In particular, equations (1a)-(1c) represent the conservation of mass, momentum and energy of the fluid, respectively. We have yet to define the momentum and energy sources in these equations. Before doing so we will introduce the radiation model.
For the radiation transport we will employ a P 1 model [17,18] that uses a P r = E r /3 closure model: where The non-dimensionalization used here haŝ In these equations E r is known as the radiation energy density and it is proportional to the zeroth angular moment of the radiation specific intensity. F r is known as the radiation flux and it is proportional to the first angular moment of the radiation specific intensity. The material absorption opacity (with units of inverse scale length) is given by σ . For convenience we have ignored scattering in our radiation transport model.

Asymptotic analysis of radiation model
In typical applications the RHD model (i.e. the coupled Euler and radiation equations) are solved in an operator split fashion where the P 1 equations (or some other transport model) is solved coupled with a material internal energy equation that contains the radiation-matter coupling terms only: The other terms in the material energy equation are updated during the hydrodynamics solve, along with a correction to take into account momentum exchange. As part of the operating splitting procedure, the radiation solve is undertaken with the density and velocity terms evaluated at a particular time level. It is these radiation equations with a quasi-static material velocity and density that we will perform an asymptotic analysis on. Consider a small, positive parameter . The P 1 equations above, coupled to equation (7), are to be scaled under the conditions that the absorption cross-section is very large: and where the ratio of the speed of light to the speed of sound is also very large: These substitutions indicate that we are considering asymptotic solutions to the transport equation for a system where the absorption opacity is large compared to our reference length scale and the material reacts much quicker to changes in radiation than any material effects such as sound waves.
Substituting these relationships into the P 1 equations and equation (7) gives where the sources are now given by All dependent variables, (E r , F r and T ), in equations (8a)-(8c) are now expanded with a formal power series of the scaling constant . For example We will now look at the coefficients for each power of . The three orders of S E are Similarly for S F : After inspecting equations (8a)-(8c), it can be seen that no terms involve −1 . This shows that the O( −1 ) equations are equal to zero, therefore S (−1) and equation (11a) gives Gathering only the terms involving no power of in the P 1 equations, also known as O(1) equations, is as follows: The first-order equation for F r0 is Substituting this relationship and equation (11b) into equation (14b) gives after some algebra. Solving for F r (1) provides a version Fick's Law, relating the radiation flux to the gradient of the energy density, with an additional drift term Next, we look at the O( ) equations arising from equations (8a)-(8c): Substituting equation (18) into equation (19a) gives Subtracting equation (19b) from equation (20) and rearranging gives 1 P Finally, this can be given in terms of material temperature due to equation (12) and our equation of state: For a known density and material velocity (as in the case of an operator split solution), the only dependent variable in this equation is the material temperature. This equation is the equilibrium drift-diffusion limit. The drift term would be absent without material motion. Though we derived these equations starting from a P 1 model for radiation transport, previous asymptotic analyses show that equilibrium diffusion is the asymptotic limit of the most general transport model [19].

The problem of radiation impinging on a cold slab
We will now consider equation (22) in one-dimensional slab geometry with constant c v and an opacity dependent on temperature to some power. This makes equation (22), dropping the (0)'s, We note that our model equation, equation (23), can be simplified to the equation solved in several previous studies of Marshak waves [2,4,6,20]. Those studies consider problems where P is negligibly small and the quantity CP = 1. It is interesting that one need not assume that v = 0 to arrive at the standard Marshak wave solution, only that P is negligible.
The problem we will solve has an initially cold, semi-infinite material located at x 0. At time t = 0 a radiation source is turned on at x = 0, heating the material there to a temperature of T = 1. We obtain a self-similar solution to this problem with the similarity transformations where A and U are constants to be defined later and θ is a parameter that gives the magnitude of the material velocity, ξ is the scaled independent variable and u is the scaled material velocity. Note that we have prescribed the relationship between the independent variables x and t, and prescribed a form for the material velocity. We do not consider how such a form for the velocity might be formed, rather we assert that if the velocity has such a 1/ √ t dependence our selfsimilar solutions are possible. Also, the material velocity is constant in space and is only a function of time.
Implementing these transforms results in Formally introducing the temperature-dependent cross section and substituting into equation (23), and setting ρc v = 1 gives At this point the following constants will be defined to simplify the arithmetic: Using these relationships, equation (26) simplifies to Once again we point out that this equation in the limit P → 0 is equivalent to the classical Marshak wave result.
We will solve equation (29) with the boundary condition that T (0) = 1 and the condition that the temperature in front of the wave is cold (i.e. T = 0). We refer to the value of ξ beyond which T = 0 as ξ max . For a particular value of ξ max , there are an infinite number of solutions that tend to zero. However, only one of these solutions maintains a zero flux in the limit ξ → ξ max . Castor [4] outlines an iterative numerical procedure beginning with an initial guess for ξ max to determine the wave temperature profile. This value for ξ max must be adjusted through trial and error until an integration from ξ max back to 0 gives T (0) = 1.

An approximation to ξ max
To find the ξ max we begin by integrating both sides of equation (29) over ξ < ξ < ξ max , giving Simplifying using T (ξ max ) = 0 yields For conciseness, the function g(ξ ) = T (ξ ) + PT 4 (ξ ) is defined. From this, the mean value of g(ξ ) from ξ to ξ max can also be defined as Using this relationship allows one to write the first two terms of equation (31) as Then, due to the fact that (ξ max − ξ ) ξ max near the radiative wave front, the second term can be neglected. Therefore, equation (33) reduces to Nelson and Reynolds illustrate the previous approximation in [6, equations (6), (7)] and state that the relative error is of order (ξ max − ξ )/ξ max . Modulo this error, equation (31) reduces to Notice that through the chain rule Because of this, equation (35) can be rewritten as Dividing both sides by T (ξ ) gives Now, applying a reverse application of the chain rule gives This reduces to A final integration over ξ and again taking advantage of the fact that T (ξ max ) = 0 gives Furthermore, evaluating the remaining integral using right-hand Riemann sums yields This method maintains accuracy to order (ξ max − ξ )/ξ max . Equation (41) now reduces to Solving this equation for T (ξ ) provides an initial temperature approximation which is defined as T 1 (ξ ): With this temperature approximation we can evaluate the integral in equation (41) instead of using right-hand Riemann sums. This will generate a more accurate expression in which we will explicitly solve for T (ξ ) below. Substituting T 1 (ξ ) into the integral term gives Evaluating the integral and simplifying yields Solving explicitly for T (ξ ) gives a more accurate approximation, defined as T 2 (ξ ) This expression allows for a more accurate approximation of the radiation wave front moving through the cold medium. In the next section, ξ max will be determined and tabulated for various values of P and θ under both n = 0 and 3.
Looking at the case when P → 0, i.e. radiation energy is negligible, equation (47) becomes which is identical to T 1 (ξ ) and the approximation given by Castor [4]. In comparison, Nelson and Reynolds [6] derive a more accurate expression for the approximate wave front.

Evaluated the self-similar solution profiles
To obtain solution profiles, T (ξ ) and ξ max for different problems we used the computer software Wolfram Mathematica 9. The solution procedure is as follows. We begin with an initial guess for ξ max and then solve equation (26) using the boundary conditions T (ξ max − δ) = T 2 (ξ max − δ) and T (ξ max − δ) = T 2 (ξ max − δ). Then, based on whether T (0) from this solution is greater than or less than 1, we adjust our guess for ξ max . In our case we used δ = 10 −10 , and the function NDSolve and FindRoot to integrate the ordinary differential equation and find the converged ξ max . Our calculations using NDSolve and FindRoot were performed with better-than machine precision arithmetic by setting the parameter WorkingPrecision to 32. To verify our solution procedure we computed ξ max using P = 0 for n = 0 and 3 and compared with the results from Nelson and Reynolds. Our solutions agree with theirs to six significant digits. Using our solutions we can look at the impact of including the radiation energy density terms in a stationary material, i.e. having P > 0 and θ = 0. To do this, we look at a typical Marshak wave solution that has appeared several times in the literature [16,21,22]. It has c v∞ = 0.1 GJ g −1 keV −1 , ρ ∞ = 3.0 g cm −3 , κ 0 = 300, T ∞ = 1.0 keV and l = 1 cm. Using these parameters, C = 948.027 and P = 0.045 73. Based on our solutions for this problem ξ max = 1.221 16 in the n = 0 case and ξ max = 1.106 22 in the n = 3 case. This compares to the 'low P' solution given by Nelson and Reynolds of ξ max = 1.231 17 in the n = 0 case and ξ max = 1.119 93 in the n = 3 case. For both of these the difference is about 1%. As shown in table 1, with no material motion the radiation wave travels slightly slower in the similarity variable ξ (but not necessarily in physical distance x) when P increases (the radiation energy relative to material energy density increases). Based on numerical results from P ranging from 0 to 2, we find that, for θ = 0, ξ max behaves as Note we have not included more decimals in the linear and quadratic coefficients of the model because the least-squares fit indicated that these were the only digits that were significant with respect to the standard error in the coefficients.
In table 1 we also show the effect of having a moving material in the entries having a non-zero θ . It is apparent, as would be expected, that having the material moving toward the right proportional to t −1/2 causes the wave to move farther. The effect of material motion is also related to P: the larger the value of P the farther the wave travels. This is due to the fact that the advection term in equation (26) is scaled by P. In the P = 0.045 73 case increasing θ to 10 from 0 has about a 15% effect in ξ max for n = 3. The same change at P = 1 leads to over 450% effect in ξ max . In figure 1 we show thermal wave profiles for P = 0.045 73 with n = 0 and 3. In these figures we can see the impact of increasing θ discussed above. The value of θ does not appear to affect the steepness of the solution at the wavefront. We did not show the smaller values of θ in these figures because their curves are very close to that for θ = 1.

Self-similar solutions for code verification
Above we have explored the impact of including both radiation energy and material motion on the classical Marshak wave problem. We also believe these solutions are useful for test problems for RHD codes. This can be accomplished by making the code fix the material velocity everywhere to have u(t) = θU/ √ t and performing the radiation solve with this prescribed u(t). In principle there could be an issue with evaluating u(0) because of the singularity, but in practice most radiation solvers expect to receive the velocity evaluated at either the end of the time step or some intermediate time level.
To demonstrate this is indeed possible we took an existing radiation solver based on the spherical harmonics approximation for slab geometry problems [16,23] and ran a problem with P = 0.045 73 and θ = 10 and n = 3. Because the code works in dimensional units we picked c v∞ = 0.1 GJ g −1 keV −1 , ρ ∞ = 3. For illustration of the difference between our solution and that without material motion we have also plotted the solution when θ = 0; this demonstrates that the material motion is noticeably affecting the solution, and if the code did not have any treatment of the relativistic terms in the coupling, the solution would be noticeably different. The numerical solution, obtained using 30 spatial cells and a time step size of 1.6667 × 10 −4 ns follows the self-similar profiles very closely.
We also did a simple convergence study for this problem by looking at the relative error at 50 ns between the numerical and semi-analytic solutions at x = 0.2 cm as a function of the mesh spacing x; see figure 3. Convergence results for a transport code can be misleading because the  x ) of 0.3. We observe first-order convergence before the mesh resolution pushes the numerical solution outside the asymptotic drift-diffusion limit. See the text for further discussion.  Table 3. Solutions T (ξ ) for a problem with P = 0.045 73 and n = 3. diffusion limit is only an asymptotic limit of the full transport system. In this problem the scaling parameter is on the order of 10 −1 and once the mesh and time step start to resolve a meanfree path or the mean-free time, the numerical method no longer is equivalent to the asymptotic drift-diffusion equation. Therefore, we expect that beyond some level of resolution, the solution will no longer converge to the self-similar profile we obtained. This is evident in our results: below a certain mesh resolution the error begins to increase. Before this point we do observe first-order convergence as we would expect for a problem with a non-smooth solution. We know of no prior work on the correct behavior of transport numerical methods in the transition region between the drift-diffusion limit. It is well known, however, that if a numerical method did not possess the proper diffusion limit, the solution would not show convergence until a mean-free path was resolved (see, for instance, [23]).
To allow other to use these solutions for code verification we have included T (ξ ) for various values θ at P = 0.045 73 in tables 2 and 3.

Conclusion
We have generalized the classic Marshak wave problem to include both radiation energy density terms and material motion. The material motion in our solutions is uniform in space and proportional to t −1/2 . In problems without any material motion, we observe that the greater the radiation energy density, as measured by the parameter P, the slower the wave moves with respect to the similarity variable ξ . We phenomenologically quantified this effect with a quadratic in P model. Besides providing insight into the effects of radiation energy and material motion on Marshak waves, our solution can also be used to verify the radiation-material coupling treatment in a simulation code.