Self-similar solutions for the non-equilibrium nonlinear supersonic Marshak wave problem

Similarity solutions to the nonlinear non-equilibrium Marshak wave problem with a time dependent radiation driving source are presented. The radiation transfer model used is the gray, non-equilibrium diffusion approximation in the supersonic regime. These solutions constitute an extension of existing non-equilibrium supersonic Marshak wave solutions which are linear, to the nonlinear regime, which prevails in realistic high energy density systems. The generalized solutions assume a material model with power law temperature dependent opacities and a material energy density which is proportional to the radiation energy density, as well as a surface radiation temperature drive which obeys a temporal power-law. The solutions are analyzed in detail and it is shown that they take various qualitatively different forms according to the values of the opacity exponents. The solutions are used to construct a set of standardized benchmarks for supersonic non-equilibrium radiative heat transfer, which are nontrivial but straightforward to implement. These solutions are compared in detail to implicit Monte-Carlo and discrete-ordinate transport simulations as well gray diffusion simulations, showing a good agreement, which demonstrates the usefulness of these solutions as a code verification test problem.

The Marshak wave, developed in the seminal work [27], which was further generalized in Refs.[28][29][30][31][32][33][34][35][36][37][38], is a fundamental concept that arises when intense energy deposition occurs in a material, leading to a steep temperature gradient.In such circumstances, radiative energy transport plays a pivotal role, and its description is far from straightforward.The Marshak wave problem poses complex questions related to the sudden energy deposition, the subsequent thermalization of the material, and the rapid emission and transport of radiative energy.At high enough temperatures, the radiative heat wave propagates faster than the speed of sound, which results in a supersonic Marshak wave [10,33,36], for which the material motion can be neglected.The original Marshak wave problem assumes local thermodynamic equilibrium (LTE) between the radiation field and the heated material, which is usually only valid at long times or for optically thick systems.Pomraning [39] and subsequently Su and Olson [40] developed a widely used solution for the non-equilibrium Marshak wave problem, by defining a material model for which the heat transfer is linear.To that end, they assumed a temperature independent * menahem.krief@mail.huji.ac.il opacity, a material energy density which is proportional to the radiation energy density and a time independent radiation temperature drive.Bennett and McClarren have recently developed a non-linear benchmark whose results are obtained numerically from detailed simulations [41].There are currently no known exact solutions to the non-equilibrium Marshak wave problem in the nonlinear regime, which prevails in most high-energy-density systems.
In this work we develop new solutions to the nonequilibrium supersonic Marshak wave problem, which are nonlinear and self-similar.The radiative transfer model used is the gray, non-equilibrium diffusion approximation.We assume a material model with power law temperature dependent opacities, which is a good approximation for many real world materials in a wide range of thermodynamic conditions.We also assume a material energy density which is proportional to the radiation energy density.It is shown that nonlinear self-similar solutions exist for a surface radiation temperature drive which obeys a temporal power-law, which is related to the absorption opacity exponent.We use the generalized solutions to define a family of standardized benchmarks for supersonic non-equilibrium radiative heat transfer.These benchmarks are compared in detail to numerical transport and gray diffusion computer simulations.

II. STATEMENT OF THE PROBLEM
In supersonic high energy density flows, where radiation heat conduction dominates and hydrodynamic motion is negligible, the material density is constant in time, and the heat flow is supersonic.Under these conditions, neglecting hydrodynamic motion, the non-equilibrium 1group (gray) radiative transfer problem in planar slab symmetry for the radiation energy density E (x, t) and the material energy density u (x, t), is formulated by the following coupled equations [42][43][44][45][46]: where c is the speed of light, k a the radiation absorption macroscopic cross section (which is also referred to as the absorption coefficient or absorption opacity), and U = aT 4 where T is the material temperature and a = 8π 5 k 4 B 15h 3 c 3 = 4σ c is the radiation constant.The effective radiation temperature T r is related to the radiation energy density by E = aT 4  r .In the diffusion approximation of radiative transfer, which is applicable for optically thick media, the radiation energy flux obeys Fick's law: where the radiation diffusion coefficient is given by: where k t = ρκ R is the total (absorption+scattering) macroscopic transport cross section, which we also refer to as the total opacity coefficient, where κ R the Rosseland mean opacity and ρ is the (time independent) material mass density.The initial conditions are of a cold material and no radiation field: E (x, t = 0) = U (x, t = 0) = 0. (5) As for the boundary conditions, we consider an imposed surface radiation temperature at x = 0 which obeys a temporal power law: T r (x = 0, t) ≡ T s (t) = T 0 t τ , (6) so that the radiation energy density at the system left boundary is: We note that this boundary condition of an imposed surface temperature is different than the common Marshak (or Milne) boundary condition which is employed in the non-equilibrium Marshak wave problem [39,40], which represents the flux incoming from a heat bath with a prescribed temperature.Nevertheless, as will be shown in section III C, a closed form relation exists between the surface radiation temperature and the heat bath temperature, so that the problem can be defined alternatively by a Marshak boundary condition with a prescribed bath temperature as a function of time.This approach was previously employed for LTE waves in Refs.[6,7,38,47].
In this work, we assume a power law temperature dependence of the total opacity, which we write as: and define a similar form for the absorption opacity: We note that the form 8 is equivalent to the common power law representation [7,10,23,33,36,37,[48][49][50] of the Rosseland opacity κ R (T, ρ) It will be noted in Appendix A that in order to obtain self-similar solutions to the problem at hand, we need to assume that u is proportional to U , that is, there should be a quartic temperature dependence for the material energy density.We write this in the common form: where ϵ is a dimensionless constant which represents the ratio between the radiation and material energies, at equilibrium.The form 10 is a special case of the general power law u (T, ρ) = FT β ρ 1−µ (see Refs. [7,10,23,33,36,37,[48][49][50]), with F = a ϵ , β = 4 and µ = 1.The temperature dependence in equation 10 is the same as employed in the theory of linear Marshak waves [39,40].
Using the material model defined in equations 8-10, equations 1-2 are written in closed form as a set of nonlinear coupled partial differential equations for E and U : where we have defined the dimensional constants:

III. A SELF-SIMILAR SOLUTION
As shown in detail in Appendix A, the problem defined by the nonlinear gray diffusion model in equations [11][12] with the initial and boundary conditions 5,7, has a selfsimilar solution, only for the specific surface temperature exponent give by: for which the following quantity is a dimensionless constant: Under the constraint 15, the problem can be solved using the method of dimensional analysis [23,36,42,51,52], resulting in a self-similar solution which is expressed in terms of self-similar profiles: with the dimensionless similarity coordinate: where the similarity exponent is: The radiation and material temperature profiles are: By plugging the self-similar form 17-18 into the nonlinear gray diffusion system 11-12 and using the relations ∂ξ ∂t = − δξ t and ∂ξ ∂x = ξ x , all dimensional quantities are factored out, and the following (dimensionless) second order ordinary differential equations (ODE) system for the similarity profiles is obtained: The surface radiation temperature boundary condition (equation 7), is written in terms of the radiation energy similarity profile as: It is evident that the dimensionless problem defined by equations 23-25 depends only on the dimensionless parameters α, α ′ , ϵ and A.
We assume that the solution has a steep nonlinear heat front, that is, f (ξ) = g (ξ) = 0 for ξ ≥ ξ 0 , where ξ 0 is finite and represents the similarity coordinate at the heat front.According to equation 19, the heat front propagates in time according to: By using equations 13,20, and 4 the heat front can be written as: which is essentially a generalization of the familiar diffusion law x F ∝ √ Dt, with a time dependent diffusion coefficient D s (t) = c/3k t (T s (t)) evaluated at the (time dependent) surface radiation temperature (see also Refs.[36,48]).Since opacities of plasmas usually decrease in magnitude with temperature, that is, α > 0, α ′ > 0, we always have δ > 1  2 and the heat propagates faster than classical diffusion: a result which is due to the temporal increase of the surface temperature (as τ = 1 α ′ > 0), which in turn increases the characteristic diffusion coefficient.According to equation 20, we have an accelerating heat front (δ > 1) if α ′ < α, a constant speed front (δ = 1) if α = α ′ and a decelerating heat front (δ < 1) if α < α ′ .In the limit α ≪ α ′ , the change in surface temperature does not increase the characteristic diffusion coefficient, and we recover the classical diffusion x F ∝ √ t propagation law, which is the familiar behavior of LTE Marshak waves with a constant driving temperature [27,28,33,35,37,45,46].
The value of ξ 0 is obtained by iterations of a "shooting method", applied on the numerical solution of the ODE system 23-24, which is integrated inwards from a trial ξ 0 to ξ = 0.The iterations adjust the trial ξ 0 until the result obeys the boundary condition from equation 25 at ξ = 0.This is essentially the same procedure which is employed in the integration of LTE Marshak waves [27,28,35,36,45,46,53].In Fig. 1, the similarity radiation and matter temperature profiles f 1/4 (ξ), g 1/4 (ξ) are presented for various cases.
We note that the self-similar solution 17-18 is a special solution of the radiative transfer problem defined by equations 1-10, for which the ratio between the material and radiation temperatures is time independent at any point ξ ∝ x/x F (t), since T (x, t) /T r (x, t) = g 1/4 (ξ) /f 1/4 (ξ) depends on x, t only through the dimensionless coordinate ξ ∝ x/x F (t).This means that the LTE limit, for which T (x, t) → T r (x, t), is never reached, even at long times.This behavior occurs for the unique value of the temporal exponent τ = 1/α ′ (equation 15), which enables a self-similar solution, by setting the temporal rate in which radiation energy enters the system via the boundary condition 7 (or equivalently, as will be discussed below, equation 37), such that the ratio between the radiation and matter energies are constant in time.We note that numerical gray diffusion simulations of the radiative transfer problem 1-10 with τ ̸ = 1/α ′ , for which the solutions are not self-similar, indeed show that the ratio T (x, t) /T r (x, t) is a function of time and space, and not a simple function of ξ ∝ x/x F (t).Moreover, it will be shown below in section III D, that for the special case α = α ′ the ratio g 1/4 (ξ) /f 1/4 (ξ) does not depend on ξ as well, so that the ratio between radiation and matter temperatures is the same at all times and at all points ξ ∝ x/x F (t) across the heat wave.

A. The LTE limit
The dimensionless constant A quantifies the strength of the emission-absorption process.If t is a typical timescale in the problem, the typical absorption coefficient is ka = k ′ 0 T 0 tτ −α ′ , and from equation 16 we find A = c ka t, which is the typical absorption-emission energy exchange rate, multiplied by the typical timescale.In general, the matter equation 2 can be written as: which shows that the quantity ϵA determines the equilibration rate.This is to be expected, since even for strong emission-absorption (A ≫ 1), a material with a large heat capacity (ϵ ≪ 1) can remain out of equilibrium, as it takes a long time to heat the material.On the other hand, even for a low emission-absorption rate (A ≪ 1), a material with a small heat capacity (ϵ ≫ 1) can reach equilibrium, as it takes a short time to be heated.Therefore, the LTE limit which prevails in the case of strong radiation-matter coupling, should be reached when ϵA ≫ 1.The fact that the quantity ϵA determines how close the problem is to LTE, is also evident in figures 1,2 and 3.In addition, equation 24 shows that when ϵA ≫ 1 we have f (ξ) ≈ g (ξ), that is, the radiation and material temperatures are approximately equal, and local equilibrium is reached.We now consider the case of strong emission-absorption (A ≫ 1) and finite ϵ, such that ϵA ≫ 1.The total energy density is with the total energy density similarity profile h (ξ) = f (ξ) + 1 ϵ g (ξ), which at LTE is simply h (ξ) = 1 + 1 ϵ f (ξ), and therefore obeys boundary condition h (0) = 1 + 1 ϵ .By summing the ODEs 23-24 in the limit f (ξ) ≈ g (ξ), a single second order ODE for h (ξ) is obtained.This equation, along with equation 19 is equivalent to the LTE Marshak wave problem (see equations 3 and 10 in Ref. [48] or equations 6 and 9 in Ref. [36]), in the specific case of a temperature surface temperature T 0 t τ with τ = 1 α ′ and a power law equation of state which includes the material and radiation energies, that is u (T, ρ) = FT β ρ 1−µ with F = a 1 + 1 ϵ , β = 4 and µ = 1.

B. The solution near the origin
As we now show, the behavior of the solution near the system boundary, can be analyzed without having to solve the full coupled ODE system 23-24.To that end, we write the first order expansion of the solution near ξ → 0: Since f (0) = 1, the value of g 0 ≡ g (0) is a measure of the deviation from equilibrium at the system's boundary: By substituting the expansion 29-30 into the matter equation 24 and equating the zero order terms, we find: Similarly, equating the first order terms gives a direct relation between the ratio of derivatives and g 0 : equation 32 is a nonlinear equation for g 0 that can be solved numerically by standard root finding algorithms, such as the Newton-Rapshon method.First we note that in the LTE limit, since ϵA ≫ 1, it is evident from equation 32 that g 0 ≈ 1 and from Eq 33 that f ′ (0) ≈ g ′ (0), as expected, since in the LTE limit we have f (ξ) ≈ g (ξ).In the general non-LTE case it is evident from equation 32 that g 0 depends on α ′ while it does not depend on α at all, and that it depends on A and ϵ only through their product ϵA.This is to be expected, since the matter-radiation equilibration process is dictated by the absorption opacity rather than the total opacity (which determines the spatial heat propagation), and since the radiation-matter coupling rate scales as the inverse material heat capacity, which is proportional to ϵ (as discussed in section III A).We also note that it is evident from equation 32, that g 0 is strictly monotonic with respect to α ′ , which is to be expected since a stronger absorption temperature dependence leads to a faster equilibration.These trends can be seen in Fig. 1, as well as in Fig. 2 where we show the resulting solutions for g 1/4 0 as a function of ϵ and A in a wide range for selected values of α ′ , and in Fig. 3 as a function of α ′ and ϵA.
Various solutions for the radiation f 1/4 (ξ) (in red) and matter g 1/4 (ξ) (in blue) temperature similarity profiles.The dimensionless defining parameters, namely the total and absorption opacity temperature powers α and α ′ , the heat capacity ratio ϵ and the coupling parameter A, are listed in the figure titles together with the resulting heat front coordinate ξ0 and temperatures ratio at the origin, g 1/4 0 .As discussed in the text, it is evident that the radiation and matter temperature become close only if ϵA is large.

C. Marshak boundary condition
As opposed to the surface radiation temperature boundary condition 7, the classical non-equilibrium Marshak wave problem [39,40], as well as many other nonequilibrium benchmarks [54][55][56][57][58][59][60][61][62], are specified in terms of a given incoming radiative flux.The former boundary condition is more natural to apply in the diffusion approximation, while the latter is more natural to use in the solution of the radiation transport equation, which has the angular surface flux as a boundary condition (see below in section IV A).Nevertheless, these two different boundary conditions are closely related.
The incoming flux boundary condition, also known as the Marshak boundary condition [39,40,47,62] where F inc (t) is a given time-dependent incoming flux at x = 0.For a medium coupled to a heat bath at temperature T bath (t), the incoming flux is F inc (t) = ac 4 T 4 bath (t).The Marshak condition 34 is an approximation of the Milne boundary condition, which is valid in the diffusion limit of radiation transport (see section IV A below).Since the surface radiation temperature is E (x = 0, t) = aT 4 s (t), the Marshak boundary condition 34 can be written as: This is a statement of the Marshak boundary condition 34, written in terms of the bath temperature, the surface radiation temperature and the net surface flux.For the self-similar problem considered in this work, the radiation surface temperature is given by equation 6, and by using equations 3-4, the radiation flux can be written in terms of the similarity profiles: so that the time dependent bath temperature, according to equation 35 is given by: where we have defined the bath constant: It is evident that only for α = α ′ the bath temperature is given by a temporal power law, which has the same temporal power τ of the surface temperature.This special case will be discussed in detail in section III D.
It is also evident that B decreases as k , so that T bath (t) ≈ T s (t) for optically opaque problems.

D. An exact analytic solution for the case α = α ′the non-LTE Henyey wave
It is shown in Appendix B, that for the special case where the absorption coefficient and the Rosseland opacity have the same temperature dependence, that is, when α = α ′ , equations 23-24 have a simple exact analytic solution of the form: where g 0 = g (0) is obtained by solving the nonlinear equation 32, and the heat front coordinate is given by: It is evident that the solution 39-40 agrees with the general relation 33, which for α = α ′ gives g ′ (0) /f ′ (0) = g 0 .We also note that the solution 39-40 is special as it results in a ξ independent ratio between the material and radiation temperatures, which does not hold in the more general case α ̸ = α ′ , for which this ratio depends on ξ.Moreover, it is evident from equation 37 that only for α = α ′ the Marshak boundary condition with a prescribed bath temperature is given by a temporal power law, which is proportional to the surface radiation temperature: where the bath constant B defined by equation 38, is and the derivative f ′ (0) can be calculated from the solution 39, to give the exact expression for the bath constant: As for the LTE limit of the exact solution with ϵA ≫ 1 and finite ϵ, we have g 0 = 1 and: In addition, in the LTE limit B → 0, so that T bath (t) → T s (t).Finally, we note that it is not surprising that an exact analytical solution exists for α = α ′ for which τ = 1 α .It is known that an exact analytic solution for the LTE radiation diffusion equation, called the Henyey Marshak wave, exists for the specific temporal exponent τ H = 1 4+α−β (see Sec. II-B of Ref. [33], Appendix A of Ref. [10] and Refs.[38,47]).Since in this work we assumed the material energy temperature power β = 4, we have τ H = 1 α = τ , so our solution for α = α ′ must reproduce the exact Henyey solution in the LTE limit.In the more general non-LTE case (a finite ϵA), the solution 39-40 in fact represents a generalization of the LTE Henyey heat wave to non-LTE gray diffusion.This generalized solution attains some characteristics of the LTE Henyey solution: it has a constant speed heat front x F (t) ∝ t and a material temperature profile of the form

E. The total and absorption optical depths
The total optical depth T , which is defined as the number of mean free paths within the heat wave is given by [47]: Taking the material temperature from the self-similar solution 22, we find: It is evident that for a Henyey wave with α = α ′ the optical depth does not depend on time.However, in this case the solution is given analytically by equation 40, for which the integral in equation 47 diverges logarithmically due to the steep temperature decrease near the front.This result is in agreement with the analysis in section IV of Ref. [47], in which LTE Henyey Marshak waves are considered, assuming a material energy density in a general temperature power law form.
A simpler estimate for the optical depth can be obtained using the mean-free-path at the system's boundary k t (T (x, t)) ≈ k t (T (x = 0, t)).Since the actual temperature profile is decreasing, this results in a useful lower bound for the optical depth which can be written in the following equivalent forms: Similarly, we define the absorption optical depth, which is the number of absorption mean free paths within the heat wave: which for our self-similar solution is given by: As for the total optical depth, if α = α ′ the absorption optical depth is time independent and the dimensionless integral diverges logarithmically.Making the same approximation by taking the absorption mean-free-path at the system's boundary k a (T (x, t)) ≈ k a (T (x = 0, t)), we obtain the following lower bound: As expected, T a increases with A, and in the LTE limit, as A → ∞, we get T a → ∞.The ratio between the typical total and absorption optical depths is given by: Finally, we write the effective optical depth (see equation 1.98 in Ref. [63]): which is evidently time-independent.The effective optical depth sets the overall thermalization rate, that is, LTE is reached if and only if T eff ≫ 1.This means that in general, LTE can be reached even when T a is not large, in a highly scattering medium with T ≫ 1 such that T eff = √ T T a ≫ 1.This can be understood in the framework of random walk, where a large number of scattering events increase the time between absorption events.However, equation 53 shows that for the specific problem defined in section II, T eff is large only if A is large and therefore, LTE will be reached only if T a is large.
In summary, the absorption and total optical thicknesses are essentially independent.Specifically, it is possible to define a heat wave which is optically thick but thin with respect to absorption, that is, the diffusion approximation holds, but the matter and radiation temperatures are significantly different.This will be demonstrated in Sec.IV.

A. Transport setup
In this section we construct a setup for a transport calculation of the diffusion problem defined in Sec.II.The general one dimensional, one group (gray) radiation transport equation in slab symmetry for the radiation intensity field I (x, µ, t) is given by [41, 44-46, 57, 62, 64]: where µ is the directional angle cosine, k a (T ) and k s (T ), are, respectively, the absorption and elastic scattering macroscopic cross sections.This transport equation for the radiation field is coupled to the material via the following material energy equation: The radiation energy density is given by the zeroth angular moment of the intensity via: and the effective radiation temperature is T r (x, t) = (E (x, t) /a) 1/4 .For optically thick problems (when the optical depth T ≫1, see Sec.III E), the diffusion limit holds, and the transport problem 54-55 can be approximated by the gray diffusion problem defined by equations 1-3, with the total opacity k t (T ) = k s (T ) + k a (T ).Hence, a transport setup of the diffusion problem defined in Sec.II should have the following effective elastic scattering opacity: We note that unless α = α ′ , this temperature dependence of the scattering opacity does model well real materials, but is used here to construct a transport problem which is equivalent to a gray diffusion problem with power law total and absorption opacities.It is also important to note that since the scattering opacity must be positive, the transport problem is well defined only if k t (T ) ≥ k a (T ) for the relevant temperatures in the problem.This constraint does not have to hold for the analogous diffusion problem, which is well defined for any k t , k a .We note that since opacity spectra of mid or high-Z hot dense materials can be extremely detailed, the total (one group "Rosseland") opacity which is dominated by spectral dips near photon energies close to 3.8k B T , is in many cases smaller than the absorption (one group "Planck") opacity, which is dominated by spectral peaks near photon energy close to 2.8k B T [44-46, 65, 66].In those realistic cases, an equivalent transport problem cannot be defined.Finally, the boundary condition for the transport problem is defined naturally by an incident radiation field for incoming directions µ > 0, which is given by a black body radiation bath: where the time dependent bath temperature is given by solution of the Diffusion problem via equation 37, which is obtained from the Marshak (Milne) boundary condition, as detailed in Sec.III C. We note that the Marshak boundary condition (equation 34) is obtained by an angular integration of the exact boundary condition given by equation 58, and by assuming a first order spherical harmonics expansion of the angular flux at the boundary, which is equivalent to the diffusion approximation (see also equations 36-38 in Ref. [62] and Sec.II in Ref. [47]).We note that since the diffusion limit is reached in optically thick problems, when T ≫1, it is expected that in this case transport simulations will agree with diffusion simulations (and the self-similar solutions).Independently, we expect the radiation and matter to be out of equilibrium for absorption thin problems where T a ≲ 1.We conclude that it is possible to construct optically thick problems with dominant scattering (k a ≪ k t ), which are absorption thin, for which we expect transport results to agree with diffusion, while the radiation and matter are significantly out of equilibrium.This reasoning is used in the construction of the test cases below.

B. Test cases
We define six benchmarks based on the self-similar solutions and specify in detail the setups for gray diffusion and transport computer simulations.We have performed gray diffusion simulations as well as deterministic discrete-ordinates (S N ) and stochastic implicit Monte-Carlo (IMC) [55,67] transport simulations.The diffusion simulations shown in this subsection were all performed without the application of flux limiters (see also subsection IV C).The S N simulations were performed using the numerical method detailed in [68], while the IMC simulations were performed using the novel numerical scheme which was recently developed by Steinberg and Heizler in Refs.[57,58,69].
The results are compared in figures 4-11, where the temperature profiles are plotted at the final time as well as at the times when the heat front reaches 20% and 60% of the final front position.
The typical scales are as follows: temperatures are in keV, time in nanoseconds and distance in centimeters.The material energy density is given by: The heat front position (equation 26) reads:  x and dimensionless coupling constant (equation 16) is: where T 0 is measured in units of keV/ns 1 α ′ , k 0 in keV α /cm and k ′ 0 in keV α ′ /cm.All tests are run until the final time t end = 1ns, and will have a surface radiation temperature of 1keV at the final time, that is, we take T 0 = 1keV/ns 1/α ′ so that the surface radiation temperature is: 1. TEST 1 In this case we take a heat capacity ratio of ϵ = 0.2 and opacities with α = α ′ = 3 and k 0 = k ′ 0 = 0.1keV 3 /cm, so that the total and absorption opacities are: which models absorption without scattering (k s (T ) = 0 for transport simulations).The surface radiation temperature is: For these parameters equation 61 gives the dimensionless coupling parameter A = 2.9979.The solution of the nonlinear equation 32 gives the boundary temperatures ratio g 1/4 0 = 0.8190643.Since α = α ′ , this case is a non-LTE Henyey wave which propagates linearly in time and has an analytical solution as described in Sec.III D. From equation 41 we find the front coordinate ξ 0 = 0.4747709 which via equation 60 gives the heat front position: From equation 44 we find the bath constant B = 1.029124, so that the bath temperature (equation 43) is: which is used in transport simulations via the incoming bath radiation flux (equation 58) or in diffusion simulations via the Marshak boundary condition (equation 34).We note that diffusion simulations can be run equivalently using the surface temperature boundary condition (equation 6).A comparison of the surface and bath temperatures as a function of time are shown in Figure 5.
The most obvious difference between the two temperatures is that the nearly 20% difference between the surface temperature and the bath temperature in this case.The radiation and matter temperature profiles are given analytically by using equations 39-40 and 21-22: Since α = α ′ , the total and absorption optical depths (see equations [48][49][50][51] are time independent and given by T = T a ≈ 0.864.
In Figure 4 we can see that the transport results disagree with the analytic and numerical diffusion solutions in the temperature profiles, but coincide with the analytic model as to the location of the wavefront and the boundary temperatures in the figure.This overall disagreement between transport and diffusion results is not surprising, since this problem is not optically thick.However, the agreement in the front position occurs since this test is opaque enough such that diffusion theory obeys the freestreaming limit, as discussed below in section IV C.

TEST 2
This case is constructed to give the optically thick limit of Test 1, by significantly increasing the total opacity while keeping the absorption opacity constant.As a result, we expect that the transport calculations will agree with diffusion and therefore with the analytic solutions.A physical mechanism for such increase is the introduction of photon scattering, which alters the total opacity only.We take the same parameters ϵ = 0.2, α = α ′ = 3, k ′ 0 = 0.1 keV 3 /cm as in Test 1, with a total opacity which is increased by a factor of 10 3 so that k 0 = 100 keV 3 /cm.Hence, the total and absorption opacities are: For transport simulations, the scattering opacity (equation 57) is now nonzero and given by the following temperature power law: The surface radiation temperature is the same as in Test 1: Since k ′ 0 is the same as in Test 1, we have also the same dimensionless coupling parameter A. Therefore, since all dimensionless parameters (α, α ′ , ϵ and A) are the same as in Test 1, this test defines the same dimensionless problem as Test 1 with the same g 1/4 0 , ξ 0 and Henyey self similar profiles: keV ( 70) However, the total optical depth, which increases with the total opacity coefficient as k 1/2 0 is now larger by a factor of √ 10 3 ≈ 31.6,so that T ≈ 27.3, while the absorption optical depth is decreased by the same factor so that T a ≈ 0.0273.Therefore, this case indeed defines an opaque heat wave so that exact transport results should coincide with the diffusion approximation.On the other hand, as in Test 1, the wave has a small absorption optical depth, which results in the same significant deviation from equilibrium as given by g 1/4 0 .The heat front position which also depends k −1/2 0 (see equation 60), runs ≈ 31.6 times slower: Similarly, equation 44 gives a smaller bath constant by the same factor, B = 0.032544, so that, as discussed in Sec.III C, the bath temperature is much closer to the surface temperature: In Figure 6 we plot the analytic solution compared with several numerical calculations.We observe that, modulo the stochastic noise in the IMC calculations, all of the calculations agree within the figure scale as expected.

TEST 3
This case is constructed to give the LTE limit of Test 2, by significantly increasing the absorption opacity, while 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 x  keeping the total opacity constant.It is defined with the same parameters ϵ = 0.2, α = α ′ = 3, and with the absorption opacity increased by a factor of 10 3 so that which as in Test 1, models pure absorption without scattering (k s (T ) = 0 for transport simulations).The surface radiation temperature is the same as in Tests 1 and 2: According to equation 61, the dimensionless coupling parameter should increase by a factor of 103 , so that A = 2997.9and ϵA ≈ 600 which is much larger then unity.Therefore, as discussed in Sec.III A, the LTE limit should be reached.The solution of equation 32 gives g 1/4 0 = 0.999446, so that the matter and radiation temperatures are almost equal.From equation 41 we find the front coordinate ξ 0 = 0.471448 which gives the heat front position: From equation 44 we find the bath constant B = 0.05954, so that the bath temperature is: The radiation and matter temperature profiles are again given analytically by: Finally, the total and absorption optical depths are T = T a ≈ 14.93.Hence, the heat wave is optically and absorption thick, which means that the diffusion approximation holds and that the problem is close to LTE.We note that this optical depth is lower than Test 2 by a factor of 0.819 3 ≈ 0.55, since the matter temperature is now almost equal to the radiation temperature, and it was increased by a factor of 1/0.819 relative to Test 2. As in Test 2, all of the numerical solutions agree with the analytic diffusion solution, as shown in Figure 7. Upon zooming into the detail of the solutions at a small value of x (see Figure 8), we observe that there is a small difference between the radiation and material temperatures predicted by the analytic solution.This difference is captured in the diffusion numerical solution as well.

TEST 4
We define another non-LTE Henyey wave similar to Test 2, that is, optically thick and absorption thin, but with a weaker opacity temperature dependence.The resulting solution profile will approach zero more gradually at the wave front than in previous cases.We take ϵ = 0.25 and smaller opacity powers α = α ′ = 1.k 0 = 10 keV 1.5 /cm, k ′ 0 = 0.1 keV 1.5 /cm, so that the total and absorption opacities are, respectively: For transport simulations, the scattering opacity (equation 57) is: The surface radiation temperature is: As k ′ 0 is the same as in tests 1 and 2, we have the same dimensionless coupling parameter A = 2.9979.The resulting value of the boundary temperatures ratio is g 1/4 0 = 0.7431154 and the front coordinate is ξ 0 = 0.877244, so that the heat front position reads: The resulting bath constant is B = 0.129865, so that the bath temperature is: The radiation and matter temperature profiles are given analytically by: The problem is indeed optically thick with a total optical depth of T ≈ 13.69, while it is absorption thin with an absorption optical depth small by a factor of k 0 /k ′ 0 = 100 (see equation 52) so that T a ≈ 0.1369.
From Figure 9 we observe the more gradual approach to zero of the solution near the wavefront.We also observe that in this problem there is increased noise in the IMC solution.This is likely due to the relatively small amount of scattering and absorption/emission behind the wavefront in this problem.

TEST 5
This case defines an optically thick and absorption thin wave (as in Test 2 and 4), but with α > α ′ .This means that in contrast to all previous cases, the heat front will not propagate at a constant speed and as it is not a Henyey solution, the self-similar profiles must be calculated numerically.We take ϵ = 1, a total opacity with α = 3.5, k 0 = 10 keV 3 /cm:  I.The temperature similarity profiles f 1/4 and g 1/4 resulting from the numerical solutions of the ODE system 23-24, as a function of ξ/ξ0, for Test 5 (ξ0 = 0.5006965) and Test 6 (ξ0 = 0.7252338).It is evident that the temperatures ratio g 1/4 (ξ) /f 1/4 (ξ) depends on ξ and is monotonically decreasing/increasing for Test 5 and Test 6, respectively.and an absorption opacity with α ′ = 1.75 and k ′ 0 = 0.1keV 3 /cm: For transport simulations, the scattering opacity (equation 57) is now a difference between two power laws: The surface radiation temperature is: The scattering opacity 84 is positive for all temperatures in the problem for short enough times, and specifically for t ≤ t end = 1 ns (since k s (T s (t end )) > 0).Since k ′ 0 is the same as in tests 1, 2 and 4, we have again A = 2.9979.The solution of the nonlinear equation 32 gives the exact boundary temperatures ratio g 1/4 0 = 0.886692.As mentioned above, since α ̸ = α ′ , the similarity profiles must be solved numerically by integrating the ODE system 23.The resulting temperature similarity profiles are given in table I.It is evident that the tabulated numerical solution agrees with the exact value of g (0) /f (0) and the exact relation 33 for g ′ (0) /f ′ (0).The front coordinate is found to be ξ 0 = 0.5006965 so that the non-linear heat front position is given by: which, since α > α ′ , accelerates in time.The numerical solution has f ′ (0) = −2.217102,and from equation 38 we find the bath constant B = 0.0970622 ns −1/2 , so that the bath temperature, which in this case is not a simple temporal power law, is given by: t ns keV.
(87) Using the self-similar solution 21-22, the radiation and material temperature profiles are given by: T r (x, t) = t ns 4 7 f 1/4 (ξ 0 x/x F (t)) keV, (88) It is evident from the tabulated solution that the ratio T (x, t) /T r (x, t) is a decreasing function of x/x F (t).This is in contrast to the previous Henyey solutions (Tests 1-4) for which this ratio is constant along the heat wave.Finally, we note that since α ̸ = α ′ , the total and absorption optical depths lower bounds (see equations 48-51) are time dependent.Their values at the end of the simulation are T (t end ) ≈ 7.63 and T a (t end ) ≈ 0.062, so that the wave is indeed optically thick and absorption thin.In figure 10 the self-similar solution is compared to numerical simulations.

TEST 6
This case defines a similar wave as in Test 5, but with α < α ′ .We take ϵ = 4, a total opacity with α = 3 and k 0 = 100keV  For transport simulations, the scattering opacity is: The surface radiation temperature is: keV.
Since α < α ′ , the scattering opacity is always negative at low enough temperatures.The coefficients k 0 , k ′ 0 were chosen such that the scattering opacity in equation 90 will become negative only at temperatures lower than 1eV.Therefore, transport simulations are initialized with a "cold" material temperature of 1eV, which is much lower than the final radiation surface temperature of 1keV,and therefore, could only cause a negligible difference between those simulations and the analytic solution.Using equation 61, we find A = 0.0029979 for this case and the solution of equation 32 gives the boundary temperatures ratio g 1/4 0 = 0.6163.As in Test 5, the similarity temperature profiles are solved numerically and given in table I, and the front coordinate is found to be ξ 0 = 0.7252338 so that the heat front position is given by: Since α < α ′ , the heat front decelerates in in time.
The numerical solution has f ′ (0) = −2.038806,and from equation 38 we find the bath constant B = 0.010065 ns 1/5 , and the bath temperature is given by: T bath (t) = 1 + 0.010065 t ns keV.
(93) The radiation and material temperature profiles are given by: The total and absorption optical depths at the end of the simulation are T (t end ) ≈ 97.9 and T a (t end ) ≈ 0.00026, so that the wave is highly optically thick and highly absorption thin, which is not surprising since in this test the scattering opacity is very large and the absorption opacity is very small.In Figure 11 we observe an agreement between the selfsimilar solution and all of the simulations, with the slight discrepancy between the S N results and the other methods at the wavefront.The S N results have a wave speed that is slightly slower than the results from the other methods.Evidence from numerical experiments indicates that this discrepancy can be mitigated by decreasing the time step size used in the calculation.Nevertheless, there is a limitation in this problem because there is numerical minimum for the initial temperature imposed to keep the opacities positive, as discussed above.

C. Flux limited diffusion
In this section we discuss how a common numerical treatment of radiation diffusion, known as flux-limited diffusion [62,70], will behave on the heat waves presented above.We begin with the radiation flux as given by Fick's law: This flux does not necessarily obey the free-streaming (causality) limit that says there cannot be a bigger flux of energy than the total amount of radiation energy present times the speed of light: Flux limiters seek to rein in the flux specified by Fick's law.One can quantify the amount of limiting required via a flux limiter parameter given by:   If R ≤ 3, we have |F | ≤ cE and no limiting is required.In order to obey the free-streaming limit, the flux limiter function λ (R) defines a flux limited diffusion flux F FL in terms of a corrected diffusion coefficient D: Using our analytical solution (equations 17,36) one finds an exact explicit expression: This expression is proportional to the inverse of the total optical depth (see equation 48): which is to be expected, since for optically thick waves a flux limiter is not needed and transport results agree with diffusion results.It is interesting to note that in the case α = α ′ , for which we know the similarity profiles analytically (see section III D), the flux limiter parameter is time and space independent: In Figure 12 we compare the results of gray diffusion simulations using several flux limiters with the analytic diffusion solution and transport simulations, for Test 1.The limiters compared are the Larsen [71] fluxlimiter and its variants, the Minerbo flux limiter [72], the Levermore-Pomraning-Zimmerman limiter [62,73], and the Kershaw limiter [72].The Larsen limiter writes the flux-limited diffusion coefficient as: which is equivalent to the flux limiter function: The value of n is a user-defined parameter and informs how the flux-limiter transitions away from standard diffusion.The limiter with n = 1 is called the "sum" limiter, and the limit as n → ∞ is known as the "max" limiter, which strictly enforces the free-streaming limit.We note that all limiters mentioned above, except the max limiter, limit the flux even for R < 3, which is not required by the free-streaming limit, but lead to better results in certain scenarios.One benefit of Test 1 is that it demonstrates how different flux limiters can give heat fronts that do not agree with either the diffusion or transport solutions.In Figure 12 we see that the Larsen with n = 4 and the max limiters agree with the transport results for the heat front position.None of the other flux-limited solutions captures the behavior of the transport solutions, and result in retarded fronts, due to the nonphysical limiting for R < 3. Interestingly, the Larsen limiter with n = 2 has been used as a default setting comme il faut [74] given that this value preserves the asymptotic diffusion limit as k t → ∞ [71] and is smoother than the max limiter.Because this problem has a large, but finite value for k t , this limiter still has an injurious effect on the heat front.This result is consistent with Ref. [75] in which radiative Marshak waves experiments were analyzed via classical diffusion, flux-limited diffusion and transport simulations.
The reason for these discrepancies can be seen by looking at how the limiters behave as a function of R. In Figure 13 we compare the strength of the limiting, λ(R), with the actual values of R for the Tests 1-4.We observe that, except for the max limiter, all of the limiters affect the solution in the test problems, with the sum limiter having the strongest effect and the Larsen with n = 4 limiter having the least, an order that is consistent with the results in Figure 12.We see that for Tests 2-4, which are optically thick, the effect of all flux limiters other than the sum limiter is negligible, as expected.Moreover, the fact that Test 1 has R < 3, means that no flux limiting is needed in the diffusion simulation, which results in the agreement of the heat front position with the transport results, shown in figure 4.However, since Test 1 defines an optically thin wave, none of the diffusion simulations (with or without flux limiters) agree with the shape of the transport Marshak wave profile.We also note that small deviations from the surface boundary temperature appear near the origin in figure 12 for the flux-limited diffusion simulations.These also result from flux-limiting, and the use of the Marshak boundary condition (equation 34).These deviations do not appear, by construction, in flux-limited diffusion simulations which are performed with a prescribed surface temperature boundary condition (equation 6), instead of the Marshak boundary condition.

V. CONCLUSION
In this work we have developed analytic and semianalytic self-similar solutions to a nonlinear nonequilibrium supersonic Marshak wave problem in the diffusion limit of radiative transfer.The solutions exist under the assumption of a material model with power law temperature dependent opacities and a material energy density which is proportional to the radiation energy density, as well as a temporal power-law surface radiation temperature drive.The solutions are a generalization of the widely used Pomraning-Su-Olson [39,40] non-equilibrium linear Marshak wave problem to the nonlinear regime.
The solutions are analyzed in detail, including a study of the LTE limit and the non-LTE optically thick and thin limits.By inspecting the solution near the origin, it is shown that the ratio between the radiation and material temperatures and derivatives near the origin, can be calculated by the root of a simple nonlinear equation.Moreover, it is shown that for the special case for which the absorption and total opacities have the same temper-ature exponents, the similarity profiles have a simple exact analytic solution, which is essentially a generalization of the well known Henyey LTE Marshak wave [33,47], to the non-LTE regime.
We constructed a set of six non-equilibrium Marshak wave benchmarks for supersonic non-equilibrium radiative heat transfer.These benchmarks were compared in detail to implicit Monte-Carlo and discrete-ordinate radiation transport simulations as well flux limited gray diffusion simulations.The first benchmark which is not optically thick, resulted, as expected, in a good quantitative agreement with the diffusion simulation and only a qualitative agreement with transport simulations.All other benchmarks were defined to be optically thick, and resulted in a very good agreement with transport simulations as well.All benchmarks except Test 3 were defined to be absorption thin, resulting in a substantial state of non-equilibrium, with a large difference between the material and radiation temperatures.This demonstrates the usefulness of the solution developed in this work as a nontrivial but straightforward to implement code verification test problem for non-equilibrium radiation diffusion and transport codes.

TFigure 3 .
Figure 3.The matter-radiation temperature ratio at the origin, as a function of the absorption opacity temperature power α ′ and ϵA.

Figure 4 .
Figure 4. Radiation and material temperature profiles for Test 1. Results are shown at times t = 0.2, 0.6 and 1ns, as obtained from a gray diffusion simulation and from Implicit-Monte-Carlo (IMC) and discrete ordinates (SN ) transport simulations, and are compared to the analytic solution of the gray diffusion equation (given in equations 65,67-68).

Figure 5 .
Figure 5.A comparison between the surface (red line) and bath (dashed line) driving temperatures for Test 1 (given in equation 66).

Figure 8 .
Figure 8.A close view near origin of the temperature profiles (from diffusion simulations and analytic solution), for Test 3 (see Fig.7).

Figure 9 .
Figure 9. Radiation and material temperature profiles for Test 4. Results are shown at times t = 0.2, 0.6 and 1ns, as obtained from a gray diffusion simulation and from Implicit-Monte-Carlo (IMC) and discrete ordinates (SN ) transport simulations, and are compared to the analytic solution of the gray diffusion equation (given in equations 80, 82-83).

Figure 10 .
Figure 10.Radiation and material temperature profiles for Test 5. Results are shown at times t = 0.341995, 0.711379 and 1ns, as obtained from a gray diffusion simulation and from Implicit-Monte-Carlo (IMC) and discrete ordinates (SN ) transport simulations, and are compared to the semi-analytic solution of the gray diffusion equation (given in equations 86, 88-89 and tableI).

Figure 12 .
Figure 12.A comparison of the radiation temperature profiles at the final time of Test 1, between flux limited Diffusion simulations with various flux limiters and the SN and IMC simulations (upper figure).The middle figure is a close up view of the heat front.The bottom figure is a close up view near the origin (excluding the noisy IMC result).

Figure 13 .
Figure 13.Various flux limiter functions λ (R).The vertical dashed black line represents the free-streaming limit R = 3.The blue vertical lines represent the (space and time independent) values of R for Tests 1-4.
Figure 11.Radiation and material temperature profiles for Test 6. Results are shown at times t = 0.133748, 0.528067 and 1ns, as obtained from a gray diffusion simulation and from Implicit-Monte-Carlo (IMC) and discrete ordinates (SN ) transport simulations, and are compared to the semi-analytic solution of the gray diffusion equation (given in equations 92, 94-95 and tableI).