RADIATION TRANSPORT FOR EXPLOSIVE OUTFLOWS: A MULTIGROUP HYBRID MONTE CARLO METHOD

Applying the IMC discretization to Equations (7) and (9) and expressing the re-balanced equations in differential form (Densmore et al. 2012) gives

$$\begin{eqnarray} C_{v}\frac{DT}{Dt} &=& f_{n}\left(\int _{4\pi }\int _{0}^{\infty } \sigma _{0,\nu _{0},a,n} I_{0,\nu _{0}}d\nu _{0}d\Omega _{0}-\sigma _{P,n}acT_{n}^{4}\right)\nonumber \\ &&-g_{0,s}^{(0)} \end{eqnarray} \tag{ 10 }$$

and

$$\begin{eqnarray} &&(1+\hat{\Omega }_{0}\cdot \boldsymbol {U}/c)\frac{1}{c}\frac{DI_{0,\nu _{0}}}{Dt}+ \hat{\Omega }_{0}\cdot \nabla I_{0,\nu _{0}}-\frac{\nu _{0}}{c}\hat{\Omega }_{0} \cdot \nabla \boldsymbol {U}\cdot \nabla _{\nu _{0}\hat{\Omega }_{0}}I_{0,\nu _{0}}\nonumber \\ &&\quad\quad + \frac{3}{c}\hat{\Omega }_{0}\cdot \nabla \boldsymbol {U}\cdot \hat{\Omega }_{0} I_{0,\nu _{0}}+(\sigma _{0,\nu _{0},s,n}+\sigma _{0,\nu _{0},a,n})I_{0,\nu _{0}}\nonumber \\ &&\quad= \frac{f_{n}}{4\pi }\sigma _{0,\nu _{0},a,n}b_{0,\nu _{0},n}acT_{n}^{4}+ \frac{b_{0,\nu _{0},n}\sigma _{0,\nu _{0},a,n}}{4\pi \sigma _{P,n}}(1-f_{n})\nonumber\\ &&\quad\quad\times\int _{4\pi } \int _{0}^{\infty }\sigma _{0,\nu _{0}^{\prime },a,n}I_{0,\nu _{0}^{\prime }}d\nu _{0}^{\prime }d\Omega _{0}^{\prime }+ \int _{4\pi }\int _{0}^{\infty }\frac{\nu _{0}}{\nu _{0}^{\prime }}\sigma _{0,s,n} \nonumber \\ &&\quad\quad\times(\boldsymbol {r},\nu _{0}^{\prime }\rightarrow \nu _{0}, \hat{\Omega }_{0}^{\prime }\cdot \hat{\Omega }_{0})I_{0,\nu _{0}^{\prime }}d\nu _{0}^{\prime }d\Omega _{0}^{\prime }, \end{eqnarray} \tag{ 11 }$$

where the integer subscript n denotes quantities evaluated at the beginning of a time step. The value $g_{0,s}^{(0)}$ has been lumped entirely into the material equation since physical scattering generally admits direct treatment in MC transport. P, R, and g subscripts indicate Planck, Rosseland, or grouped quantities, respectively. The opacities are evaluated in the comoving frame. The value $b_{0,\nu _{0}}$ is the frequency-normalized Planck function in the comoving frame, and the Fleck factor f_n (Fleck & Cummings 1971) is

$$\begin{equation} f_{n}=\frac{1}{1+\alpha \beta _{n}\sigma _{P,n}c\Delta t_{n}}. \end{equation} \tag{ 12 }$$

The value α ∈ [0, 1] is a time centering control parameter (often set to 1), $\beta _{n}=4aT_{n}^{3}/C_{v,n}$ and Δt_n is the physical time step size for time step n. The second term on the right-hand side of Equation (11) is the source due to effective scattering (Fleck & Cummings 1971; Densmore et al. 2012). The differential effective scattering opacity, $(1-f_{n})b_{0,\nu _{0},n}\sigma _{0,\nu _{0},a,n} \sigma _{0,\nu _{0}^{\prime },a,n}\nu _{0}^{\prime }/4\pi \sigma _{P,n}\nu _{0}$, is separable in ν₀ and $\nu _{0}^{\prime }$; so the new frequency of a photon undergoing effective scattering is probabilistically independent of the old frequency (Mihalas & Mihalas 1984, p. 327).

Equation (11) could in principle be replaced with the fully relativistic comoving transport equation described by Mihalas & Mihalas (1984, p. 434). Here, we consider only first-order relativistic effects but note that higher order effects could be incorporated with relative ease in MC codes (Abdikamalov et al. 2012). The typical outflow speed of SNe Ia is approximately U_max = 30, 000 km s⁻¹, so (U_max/c)² ≈ 0.01. The lab frame Eulerian coordinate is related to the velocity through

$$\begin{equation} \boldsymbol {r}=\boldsymbol {U}(t_{n}+t_{\min }), \end{equation} \tag{ 13 }$$

where $t_{n}\, {=}\, \sum _{n^{\prime }=1}^{n-1}\Delta t_{n}$ and t_min is the starting time of the expansion.

The right-hand side of Equation (11) implies source MC particles may be generated, scattered, or absorbed in the comoving frame in the same manner as the static material IMC method (Abdikamalov et al. 2012; Lucy 2005). To seed interaction locations between events in a particle's history, the particle's properties may be mapped to the lab frame and streamed by Equation (1). Neglecting all terms of order U²/c² or higher, the frequency and direction transformations are (Castor 2004, p. 103)

$$\begin{equation} \nu =\nu _{0}\left(1+\frac{\hat{\Omega }_{0}\cdot \boldsymbol {U}}{c}\right) \end{equation} \tag{ 14 }$$

and

$$\begin{equation} \hat{\Omega }=\frac{\hat{\Omega }_{0}+\boldsymbol {U}/c}{1+\hat{\Omega }_{0}\cdot \boldsymbol {U}/c} \end{equation} \tag{ 15 }$$

respectively. The cumulative effect of these transformations in each particle history accounts for Doppler shifting and aberration in the MC process (Lucy 2005). According to (Pomraning 1973, p. 156) and (Castor 2004, p. 104), the opacity in the lab frame is

$$\begin{equation} \sigma _{{\rm lab}} = \frac{\nu _{0}}{\nu }\sigma _{{\rm cmf}}, \end{equation} \tag{ 16 }$$

where "cmf" stands for comoving frame. A computational convenience may be taken by virtue of Equation (13). Following Kasen, we track both IMC and DDMC particles in velocity space. Thus, the grid itself is logically unchanging and velocity acts as a Lagrangian coordinate. IMC particles move over "velocity distances" denoted with a lower case u. Corresponding to physical distance in standard IMC (Fleck & Cummings 1971), there is a velocity to a cell boundary u_b, a velocity to collision u_col, and a velocity distance to census at the end of a time step u_cen. Through Equation (13), the velocity to a boundary can be calculated with the same formula as the physical distance and is consequently dependent on grid geometry. For the u_col and u_cen, the formulae are

$$\begin{equation} u_{{\rm col}}=\frac{-\ln (\xi)}{(t_{n}+t_{\min })(1-\hat{\Omega }_{p} \cdot {{\boldsymbol U}_{p}}/c)(\sigma _{s,n}+\sigma _{a,n})} \end{equation} \tag{ 17 }$$

and

$$\begin{equation} u_{{\rm cen}}=c\left(\frac{t_{n}+\Delta t_{n}-t_{p}}{t_{n}+t_{\min }}\right), \end{equation} \tag{ 18 }$$

where t_p, U_p, and $\hat{\Omega }_{p}$ are a particle's time, "velocity position," and direction in the lab frame, respectively. The value ξ ∈ (0, 1] is a uniformly sampled random number. The minimum u = min (u_col, u_cen, u_b) indicates which event occurs at each iteration of a particle's history.

Effective absorption can be alternatively calculated with implicit capture (Lewis & Miller 1993, p. 332). If implicit capture is used to reduce the variance of an MC particle tally, the collision velocity only includes scattering opacities (Fleck & Cummings 1971),

$$\begin{equation} u_{{\rm col}}=\frac{-\ln (\xi)}{(t_{n}+t_{\min })(1-\hat{\Omega }_{p} \cdot {{\boldsymbol U}_{p}}/c)(\sigma _{s,n}+(1-f_{n})\sigma _{a,n})}. \end{equation} \tag{ 19 }$$

The energy of a particle in the lab frame is reduced by $E_{p}\rightarrow E_{p}e^{-f_{n}\frac{\nu _{0,p}}{\nu _{p}} \sigma _{a,n}u(t_{n}+t_{\min })}$ where ν_p and ν_{0, p} are the particle's lab and comoving frame frequency, respectively.

For a one-dimensional (1D), spherically symmetric shell geometry,

$$\begin{equation} u_{b}=\left\lbrace \begin{array}{@{}l}|\left(U_{j-1/2}^{2}-\left(1-\mu _{p}^{2}\right)U_{p}^{2}\right)^{1/2}+\mu _{p}U_{p}| \\\ms \;\;\;\;{\rm if }\mu _{p}<-\sqrt{1-(U_{j-1/2}/U_{p})^{2}}\\\ms \left(U_{j+1/2}^{2}-\left(1-\mu _{p}^{2}\right)U_{p}^{2}\right)^{1/2}-\mu _{p}U_{p}\\\ms \;\;\;\;{\rm otherwise} \end{array},\right. \end{equation} \tag{ 20 }$$

where μ_p is the projection of $\hat{\Omega }_{p}$ along the radial coordinate. The index j ∈ {1...J} denotes a velocity zone. For a geometry in which the inner most cell has U_1/2 = 0, the second case in Equation (20) must be applied for the innermost cell, j = 1, and μ_p ∈ [ − 1, 1]. IMC particle position in velocity space must be updated to have its physical position unchanged. The censused position is simply maintained with

$$\begin{equation} \boldsymbol {U}_{p,{\rm new cen}}(t_{n}+\Delta t_{n}+t_{\min }) = \boldsymbol {U}_{p,{\rm old\; cen}} (t_{n}+t_{\min })=\boldsymbol {r}_{p,{\rm cen}}, \end{equation} \tag{ 21 }$$

where r_p is the implied IMC particle position at the end of a time step. Equation (21) can either be implemented before or after the routine that advances the particle set. We find that introducing a time centering parameter, α₂, to split the velocity position shift before and after transport is generally preferable. At particle advance, the algorithm is

• 1.  
  $\boldsymbol {U}_{p,{\rm * cen}}(t_{n}+\alpha _{2}\Delta t_{n}+t_{\min }) = \boldsymbol {U}_{p,{\rm old\; cen}}(t_{n}+t_{\min })$.
• 2.  
  IMC: $\boldsymbol {U}_{p,{\rm * cen}} \rightarrow \boldsymbol {U}_{p,*}$.
• 3.  
  $\boldsymbol {U}_{p,{\rm new\; cen}}(t_{n}+\Delta t_{n}+t_{\min }) = \boldsymbol {U}_{p,*}(t_{n}+\alpha _{2}\Delta t_{n}+t_{\min })$.

If we were transporting massive particles, then the above steps would ensure that the particle does not move in space if it has zero velocity with respect to the lab frame.

In many transport problems, the exact dependence of opacity on frequency is not necessarily known. If the fully continuous opacities are known, they may not be practical to implement in analytic form. One might then a priori posit a group structure for the opacities that are relevant to the transport problem and formulate Equation (11) in terms of such a group structure. We may then set the group Rosseland and Planck opacities, σ_{P, g, n} = σ_{R, g, n} = σ_{a, g, n}. But if each velocity cell has its own frequency grouping and number of groups, we must constrain the resulting equation to one particular fluid cell. For a cell j, group index g, and a total number of (transport) groups G_j,

$$\begin{eqnarray} &&(1+\hat{\Omega }_{0}\cdot \boldsymbol {U}/c)\frac{1}{c}\frac{DI_{0,g}}{Dt}+\hat{\Omega }_{0} \cdot \nabla I_{0,g}+\frac{1}{c}\hat{\Omega }_{0}\cdot \nabla \boldsymbol {U}\cdot \hat{\Omega }_{0} \nonumber \\\ms &&\quad\quad\times(I_{0,g}-\nu _{g-1/2}I_{0,\nu _{g-1/2}}+\nu _{g+1/2}I_{0,\nu _{g+1/2}})\nonumber \\\ms &&\quad\quad-\frac{1}{c}\hat{\Omega }_{0}\cdot \nabla \boldsymbol {U}\cdot (\mathbf {I}-\hat{\Omega }_{0} \hat{\Omega }_{0})\cdot \nabla _{\hat{\Omega }_{0}}I_{0,g}\nonumber \\\ms &&\quad\quad+\frac{3}{c}\hat{\Omega }_{0} \cdot \nabla \boldsymbol {U}\cdot \hat{\Omega }_{0}I_{0,g}+(\sigma _{s,g,n}+\sigma _{a,g,n})I_{0,g}\nonumber\\\ms &&\quad=\frac{f_{n}}{4\pi }\gamma _{g,n} \sigma _{P,n}acT_{n}^{4}+\frac{1}{4\pi }\gamma _{g,n}(1-f_{n})\nonumber \\\ms &&\quad\quad\times\sum _{g^{\prime }=1}^{G_{j}}\int _{4\pi } \sigma _{a,g^{\prime },n}I_{0,g^{\prime }}d\Omega _{0}^{\prime }\nonumber \\\ms &&\quad\quad+\sum _{g^{\prime }=1}^{G_{j}}\int _{4\pi }\sigma _{s,n}(g^{\prime }\rightarrow g,\hat{\Omega }_{0}^{\prime } \cdot \hat{\Omega }_{0})I_{0,g^{\prime }}d\Omega _{0}^{\prime }, \end{eqnarray} \tag{ 22 }$$

where g ∈ {1...G_j}, U ∈ [U_{j − 1/2}, U_{j + 1/2}] and $\gamma _{g,n}= \int _{\nu _{g+1/2}}^{\nu _{g-1/2}}\sigma _{a,n}b_{0,\nu _{0},n}d\nu _{0}/\sigma _{P,n}$. Additionally, $I_{0,g}=\int _{\nu _{g+1/2}}^{\nu _{g-1/2}}I_{0,\nu _{0}}d\nu _{0}$ and a higher group index corresponds to a lower frequency or equivalently a higher wavelength. I is the identity matrix. We have made sure to apply the group grid in the comoving frame. Equation (22) makes use of Castor's division of the Doppler and aberration corrections from the photon momentum gradient. The third term on the left-hand side of Equation (22) indicates that Doppler shifting will cause particles to leak between groups. Moreover, the leakage process will depend on the anisotropy of the radiation field as well as the velocity gradient. We track frequency or wavelength continuously. Hence when an IMC particle crosses a cell boundary, its group in the subsequent cell environment can be determined without ambiguity. If groups in adjacent cells do not align, group intersections determine transitions in phase space for IMC and DDMC particles (Densmore et al. 2012). In-cell group transitions from streaming in velocity may be modeled directly by computing a velocity distance to Doppler shift, u_Dop, and taking u = min (u_col, u_cen, u_b, u_Dop). The form of the velocity to Doppler shift in a spherically symmetric outflow is $u_{{\rm Dop}}=|c(1-\nu _{g+1/2}/\nu _{p})-\boldsymbol {U}_{p}\cdot \hat{\Omega }_{p}|$ where ν_p is particle lab frame frequency and use has been made of Equations (14) and (15).

Integrating Equation (22) over comoving solid angle, using the simplifications described in Castor (2004, pp. 102–111), and assuming comoving isotropic opacities,

$$\begin{eqnarray} &&\frac{1}{c}\frac{\partial \phi _{0,g}}{\partial t}+\frac{1}{c} \nabla \cdot (\boldsymbol {U}\phi _{0,g})+\nabla \cdot \boldsymbol {F}_{0,g}\nonumber \\ &&\quad\quad+(\mathbf {P}_{0,g}-\nu _{g-1/2}\mathbf {P}_{0,\nu _{g-1/2}}+\nu _{g+1/2} \mathbf {P}_{0,\nu _{g+1/2}}):\nabla \boldsymbol {U}\nonumber \\ &&\quad\quad+\left[\sigma _{s,g,n}+f_{n}\sigma _{a,g,n}+(1-\gamma _{g,n})(1-f_{n})\sigma _{a,g,n}\right]\phi _{0,g}\nonumber \\ &&\quad=f_{n}\gamma _{g,n}\sigma _{P,n}acT_{n}^{4}+\gamma _{g,n}(1-f_{n})\sum _{g^{\prime }\not=g}^{G_{j}} \sigma _{a,g^{\prime },n}\phi _{0,g^{\prime }}\nonumber \\ &&\quad\quad+\sum _{g^{\prime }=1}^{G_{j}}\sigma _{s,n}(g^{\prime }\rightarrow g)\phi _{0,g^{\prime }}, \end{eqnarray} \tag{ 23 }$$

where

$$\begin{equation*} \phi _{0,g}=\int _{4\pi }I_{0,g}d\Omega _{0},\nonumber \end{equation*}$$

$$\begin{equation*} \boldsymbol {F}_{0,g}=\int _{4\pi }\hat{\Omega }_{0}I_{0,g}d\Omega _{0},\nonumber \end{equation*}$$

$$\begin{equation*} \mathbf {P}_{0,g}=\frac{1}{c} \int _{4\pi }\hat{\Omega }_{0}\hat{\Omega }_{0}I_{0,g}d\Omega _{0},\nonumber \end{equation*}$$

$$\begin{equation*} \mathbf {P}_{0,\nu _{g\pm 1/2}}=\frac{1}{c} \int _{4\pi }\hat{\Omega }_{0}\hat{\Omega }_{0}I_{0,\nu _{g\pm 1/2}}d\Omega _{0},\nonumber \end{equation*}$$

and the expression W: V is the trace of the matrix product WV^T, i.e., $\mathbf {W}:\mathbf {V}=\sum _{k}\sum _{k^{\prime }}W_{k,k^{\prime }}V_{k,k^{\prime }}$ (Castor 2004, p. 111). If the physical scattering is elastic, then the last term on the right-hand side of Equation (23) cancels with σ_{s, g, n}ϕ_{0, g}.

Equation (23) is a starting point for describing multigroup DDMC (Abdikamalov et al. 2012). The frequency groups for DDMC do not necessarily have to be the same as those for IMC. If group lumping is employed over regimes of energy and velocity cells that are amenable to a diffusion approximation, then the DDMC group structure will be different (less resolved) than the IMC group structure. The transport and diffusion groups however must complement each other over the prescribed (user defined) energy grid. Following Abdikamalov, we operator split Equation (23) into a transport component, a Doppler component, and an advection-expansion component (Abdikamalov et al. 2012):

$$\begin{eqnarray} &&\frac{1}{c}\left(\frac{\partial \phi _{0,g}}{\partial t}\right)_{{\rm Transport}}+\nabla \cdot \boldsymbol {F}_{0,g} \nonumber \\ &&\quad\quad+\left[\sigma _{s,g,n}+f_{n}\sigma _{a,g,n}+(1-\gamma _{g,n})(1-f_{n})\sigma _{a,g,n}\right] \phi _{0,g} \nonumber \\ &&\quad=f_{n}\gamma _{g,n}\sigma _{P,n}acT_{n}^{4}+\gamma _{g,n}(1-f_{n}) \sum _{g^{\prime }\not=g}^{G_{j}}\sigma _{a,g^{\prime },n}\phi _{0,g^{\prime }}\nonumber \\ &&\quad\quad+\sum _{g^{\prime }=1}^{G_{j}}\sigma _{s,n}(g^{\prime }\rightarrow g)\phi _{0,g^{\prime }}, \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} &&\frac{1}{c}\left(\frac{\partial \phi _{0,g}}{\partial t}\right)_{{\rm Doppler}}+\mathbf {P}_{0,g}:\nabla \boldsymbol {U}\nonumber \\ &&\quad=(\nu _{g-1/2}\mathbf {P}_{0,\nu _{g-1/2}}-\nu _{g+1/2}\mathbf {P}_{0,\nu _{g+1/2}}):\nabla \boldsymbol {U}, \end{eqnarray} \tag{ 25 }$$

and

$$\begin{equation} \left(\frac{\partial \phi _{0,g}}{\partial t}\right)_{{\rm Adv/Exp}}+\nabla \cdot (\boldsymbol {U}\phi _{0,g}) = 0. \end{equation} \tag{ 26 }$$

Equations (23)–(26) were obtained by taking the zeroth moment of Equation (22) in solid angle. Multiplying Equation (22) by $\hat{\Omega }_{0}$, integrating over solid angle, and using Buchler's analysis to drop insignificant terms (Buchler 1983), the operator split first moment equations are (Castor 2004; Abdikamalov et al. 2012)

$$\begin{equation} \frac{1}{c}\left(\frac{\partial \boldsymbol {F}_{0,g}}{\partial t}\right)_{{\rm Transport}} +c\nabla \cdot \mathbf {P}_{0,g}=-(\sigma _{s,g,n}+\sigma _{a,g,n})\boldsymbol {F}_{0,g} \end{equation} \tag{ 27 }$$

and

$$\begin{equation} \left(\frac{\partial \boldsymbol {F}_{0,g}}{\partial t}\right)_{{\rm Adv/Exp}}+\nabla \cdot (\boldsymbol {U}\boldsymbol {F}_{0,g})=0, \end{equation} \tag{ 28 }$$

where now terms with the Fleck factor are no longer present due to isotropy. Equations (24) and (27) are merely the P₁ equation set over the Eulerian coordinate corresponding to multigroup IMC transport. If the intensity is only linearly anisotropic and the flux varies slowly with respect to photon mean free time, then Equation (27) reduces to Fick's Law,

$$\begin{equation} \boldsymbol {F}_{0,g}=\frac{-1}{3(\sigma _{s,g,n}+\sigma _{a,g,n})}\nabla \phi _{0,g}. \end{equation} \tag{ 29 }$$

Additionally, Equation (25) becomes

$$\begin{eqnarray} &&\frac{\partial \phi _{0,g}}{\partial t}+\frac{\nabla \cdot \boldsymbol {U}}{3} \phi _{0,g}=\frac{\nabla \cdot \boldsymbol {U}}{3}(\nu _{g-1/2}\phi _{0,\nu _{g-1/2}}- \nu _{g+1/2}\phi _{0,\nu _{g+1/2}}).\nonumber\\ \end{eqnarray} \tag{ 30 }$$

The relevant equations for the DDMC method are Equations (24), (26), (28), (29), and (30). From the MC perspective, Equations (26) and (28) are both solved by advecting the particles along with the fluid. But we are transporting particles over the space of velocities in an outflow. So to solve Equation (26) or (28), we merely leave the DDMC particles in the cells where they are censused. Since IMC does not call for explicit advection of MC particles, explicit advection is DDMC specific.

Equation (30) can be solved in several ways. One might approximate $\phi _{0,\nu _{g-1/2}}\approx \phi _{0,g-1}/\Delta \nu _{g-1}$ if ∇ · U is positive and the local radiation field is cumulatively redshifting, or $\phi _{0,\nu _{g-1/2}}\approx \phi _{0,g}/\Delta \nu _{g}$ if ∇ · U is negative and the local radiation field is blueshifting (Δν_g = ν_{g − 1/2} − ν_{g + 1/2}). Incorporating the approximation for a redshifting field into Equation (30), the resulting equation for the lowest group index is a homogeneous ordinary differential equation (ODE). The solution to g = 1 ODE can be used to calculate the heterogeneity from between-group Doppler shifting in the remaining group equations. Alternatively, one might apply Abdikamalov's approach of tracking particle frequency continuously for both IMC and DDMC (Abdikamalov et al. 2012). The homogeneous solution to Equation (30) can be used to shift the energy weight and wavelengths of each MC particle. Between-group shifting is then obtained by regrouping particle histories based on the new value of the particle frequency. Regrouping particles after shifting frequency accounts for the right-hand side of Equation (30).

To obtain the DDMC method, Equation (29) must be substituted into Equation (24) and the result must be spatially discretized (Densmore et al. 2008). For gray diffusion, a DDMC cell may be adjacent to an IMC cell, the domain boundary, or another DDMC cell (Densmore et al. 2007). For multigroup diffusion, a cell might use DDMC in one group and IMC in another. Using notation similar to Densmore et al. (2012), a DDMC equation for a cell j in a fully optically thick domain away from any domain boundaries is

$$\begin{eqnarray} &&\frac{1}{c}\frac{\partial \phi _{0,j,g}}{\partial t}+\left [\sum _{j^{\prime }} \sigma _{j\rightarrow j^{\prime },g}+\sigma _{s,j,g,n}+f_{j,n}\sigma _{a,j,g,n}+(1-\gamma _{j,g,n})\right. \nonumber \\ &&\left. \quad\quad\times(1-f_{j,n})\sigma _{a,j,g,n}\right ]\phi _{0,j,g}= f_{j,n}\gamma _{j,g,n}\sigma _{P,j,n}acT_{j,n}^{4}\nonumber \\ &&\quad\quad+\frac{1}{V_{j}}\sum _{j^{\prime }}V_{j^{\prime }}\sigma _{j^{\prime }\rightarrow j,g}\phi _{0,j^{\prime },g}+ \gamma _{j,g,n}(1-f_{j,n})\nonumber \\ &&\quad\quad\times\sum _{g^{\prime }\not=g}^{G_{j}}\sigma _{a,j,g^{\prime },n}\phi _{0,j,g^{\prime }}+\sum _{g^{\prime }=1}^{G_{j}}\sigma _{s,j,n}(g^{\prime }\rightarrow g)\phi _{0,j,g^{\prime }}, \end{eqnarray} \tag{ 31 }$$

where j' is the index of a cell that shares a face with j, $\sigma _{j\rightarrow j^{\prime },g}$ are determined from the finite volume discretization of the divergence of the flux along with Fick's law (so they are grid geometry dependent), and V_j is the volume of cell j. In Equation (31) and in what follows, we have dropped the "Transport" subscript from ∂ϕ_{0, j, g}/∂t for simplicity. Having applied Abdikamalov's operator split, the usual DDMC interpretation may be given to Equation (31). Specifically, the "leakage opacities" $\sigma _{j\rightarrow j^{\prime }}$ determine how likely a DDMC particle will leak from cell j to cell j', f_{j, n}σ_{a, j, g, n} is the effective absorption opacity, and (1 − γ_{j, g, n})(1 − f_{j, n})σ_{a, j, g, n} determines how likely a DDMC particle will scatter out of its current group (Densmore et al. 2012). The source terms on the right-hand side of Equation (31) are, respectively, effective thermal emission, particles leaking into j from adjacent cells, in-scattering from groups of cell j into g, and physical scattering.

A DDMC particle's event may be sampled from a histogram of the opacities multiplying ϕ_{0, j, g} in the second term of Equation (31) (Densmore et al. 2007). Since the diffusion equation is kept continuous in time, each DDMC particle is thought to stream in time (Densmore et al. 2007). If scattering is elastic, the time to a next event is (Densmore et al. 2007, 2012)

$$\begin{equation} \delta t_{p} = -\frac{1}{c}\frac{\ln (\xi)}{\sigma _{{\rm DDMC total}}}, \end{equation} \tag{ 32 }$$

where ξ ∈ (0, 1] is again a uniformly sampled random variable and $\sigma _{{\rm DDMC \; total}}=\sum _{j^{\prime }}\sigma _{j\rightarrow j^{\prime },g}+f_{j,n}\sigma _{a,j,g,n} +(1-\gamma _{j,g,n}) (1-f_{j,n})\sigma _{a,j,g,n}$.

Having discretized the diffusion equation in space and frequency, we must make note of some important ambiguities that affect the implementation of this hybrid diffusion transport method. The first ambiguity is DDMC particle position which we touched upon in the introduction. An additional ambiguity implied by Equation (31) is DDMC particle frequency or wavelength. Specifically, the scattering opacity only accounts for particles leaving their current group. Densmore's form of multifrequency IMC–DDMC employed a gray DDMC method for particles transporting below a threshold frequency (Densmore et al. 2012). Applying DDMC over arbitrary group distributions is then a generalization to opacities that may be strongly non-monotonic in frequency. But this slight modification to the method does not change the notion that DDMC is essentially gray during in-group propagation. So DDMC particles may have continuous frequency as a property as long as it is re-sampled within the current group before being explicitly used in the transport-diffusion algorithm (Densmore et al. 2012). To not re-sample DDMC frequencies at every reference in the code is to neglect the possibility that multiple scattering processes occurred within the group.

If a particular DDMC cell-group, (j, g), is adjacent to an IMC cell (j', g') where the groups of g and g' are contiguous, then one may use either a Marshak boundary condition or the Habetler & Matkowsky (1975) asymptotic diffusion limit boundary condition. In either case, the boundary condition for the split transport operator may be expressed as (Densmore et al. 2008)

$$\begin{eqnarray} &&\phi _{0,g}(\boldsymbol {r}_{b},t)+\left(\frac{\lambda }{\sigma _{a,g,n}+\sigma _{s,g,n}}\right)\boldsymbol {n}\cdot \nabla \phi _{0,g}(\boldsymbol {r}_{b},t)\nonumber \\ &&\quad=2\int _{\hat{\Omega }_{0}\cdot \boldsymbol {n}<0}\int _{\nu _{g+1/2}}^{\nu _{g-1/2}} W(|\hat{\Omega }_{0}\cdot \boldsymbol {n}|)I_{0,\nu _{0}} (\boldsymbol {r}_{b},\hat{\Omega }_{0},t)d\nu _{0}d\Omega _{0}\quad\quad \end{eqnarray} \tag{ 33 }$$

in the comoving frame over Eulerian coordinates, where λ ≈ 0.7104, n is the unit outward normal vector of a cell surface, and r_b is a point on the cell surface (Densmore et al. 2008). The function W(μ) = 2μ for an isotropic Marshak boundary condition or W(μ) ≈ μ + 3μ²/2 for an approximation to Habetler's asymptotic result (Densmore et al. 2007; Habetler & Matkowsky 1975). We have tacitly assumed that diffusive scattering between groups during a boundary crossing does not occur. The near-equilibrium condition at the cell boundaries is reasonable (Habetler & Matkowsky 1975) when the mean free time is small compared to time step size. A small relative mean free time can be enforced by the same heuristic used to determine if a cell is DDMC compatible.

Numerical experiments in Section 3.3 indicate that Equation (33) might be insufficient for some mean free path thresholds dictating whether a cell-group is in IMC or DDMC. Specifically, at least in spherical geometry we have found that applying DDMC in only cell-groups with large (≳ 10) numbers of radial mean free paths in a problem with a monotonic density gradient and a strong outflow (∼10⁹ cm s⁻¹) may cause an artificial depression in the hybrid radiation energy density profiles where the code is applying Equation (33). It has been noted by Densmore et al. (2008) that Equation (33) does not incorporate effects of curvature. Malvagi & Pomraning (1991) derive an asymptotic diffusion limit boundary condition that expands on the work of Habetler & Matkowsky (1975) by incorporating spatial curvature, spatial variation, and opacity variation at the boundary. For the set of Heaviside source outflow problems along with the range of cell-group coupling heuristics we test, the effect of spatial curvature at IMC–DDMC boundaries is found to be negligible. Nevertheless, the work of Malvagi & Pomraning (1991) and Case (1960) may be extended to incorporate fluid effects. Instead of an additional curvature term in Equation (33), we apply a boundary layer asymptotic analysis to the O(U/c) comoving transport equation with frequency independence to obtain an expansion term that may then be given a MC interpretation. The assumption that scattering does not occur simultaneously with DDMC–IMC boundary interactions then allows for trivial extension to multigroup.

Neglecting O(U²/c²) terms and assuming frequency independent opacity, integrating the comoving transport equation gives

$$\begin{eqnarray} &&\frac{1}{c}\frac{\partial I_{0}}{\partial t}+ \hat{\Omega }_{0}\cdot \nabla I_{0}+\sigma _{t,0}I_{0}+ \frac{\boldsymbol {U}}{c}\cdot \nabla I_{0}-\frac{1}{c}\hat{\Omega }_{0} \cdot \nabla \boldsymbol {U}\nonumber \\ &&\quad\quad\cdot (\mathbf {I}-\hat{\Omega }_{0} \hat{\Omega }_{0})\cdot \nabla _{\hat{\Omega }_{0}}I_{0}+ \frac{4}{c}\hat{\Omega }_{0}\cdot \nabla \boldsymbol {U}\cdot \hat{\Omega }_{0}I_{0} = j_{0},\quad\quad \end{eqnarray} \tag{ 34 }$$

where $I_{0}=\int _{0}^{\infty }I_{0,\nu _{0}}d\nu _{0}$, j₀ is the total frequency integrated source due to scattering events and emission and σ_{t, 0} = σ_{s, 0} + σ_{a, 0} is a total opacity. Supposing there exists some spatial surface denoted by b, we now make use of the homologous outflow Equation (13) to obtain

$$\begin{eqnarray} &&\frac{1}{c}\frac{\partial I_{0}}{\partial t}+ \mu \frac{\partial I_{0}}{\partial z}+\mathcal {L}I_{0}+ (\hat{\Omega }_{0}\cdot \nabla)_{\perp }I_{0}\nonumber \\ &&\quad\quad+\sigma _{t,0}I_{0}+\frac{\boldsymbol {r}}{ct_{f}} \cdot \nabla I_{0}+\frac{4}{ct_{f}}I_{0} = j_{0}, \end{eqnarray} \tag{ 35 }$$

where μ is the projection of the comoving angle onto an axis z aligned orthogonal to a plane tangent to surface b, the linear operator $\mathcal {L}$ accounts for the change in the directional derivative, $\hat{\Omega }_{0}\cdot \nabla$, due to non-trivial coordinate curvature, $(\hat{\Omega }_{0}\cdot \nabla)_{\perp }$ is the projection of the directional derivative orthogonal to z (Malvagi & Pomraning 1991), and the "fluid time" t_f = t_n + t_min upon implementation. Now we take the usual step of postulating a parameter ε ≪ 1 such that (Habetler & Matkowsky 1975; Malvagi & Pomraning 1991)

$$\begin{eqnarray} &c\rightarrow c/\varepsilon, \end{eqnarray} \tag{ 36a }$$

$$\begin{eqnarray} &\sigma _{t,0}\rightarrow \sigma _{t,0}/\varepsilon, \end{eqnarray} \tag{ 36b }$$

$$\begin{eqnarray} &\sigma _{a,0}\rightarrow \varepsilon \sigma _{a,0}, \end{eqnarray} \tag{ 36c }$$

$$\begin{eqnarray} &t_{f}\rightarrow t_{f}/\varepsilon, \end{eqnarray} \tag{ 36d }$$

$$\begin{eqnarray} &\boldsymbol {r}_{b}\rightarrow \boldsymbol {r}_{b}/\varepsilon, \end{eqnarray} \tag{ 36e }$$

$$\begin{eqnarray} &s = z/\varepsilon, \end{eqnarray} \tag{ 36f }$$

$$\begin{eqnarray} &\boldsymbol {U}(\boldsymbol {r},t_{f})\approx \boldsymbol {U}(\boldsymbol {r}_{b},t_{f}), \end{eqnarray} \tag{ 36g }$$

where r_b is a location on surface b and the rescaling in Equation (36) permits treating the parameters as O(1). We have scaled the surface coordinate and characteristic fluid time to be large in Equations (36d) and (36e). If $\hat{r}_{b}$ and $\hat{z}$ are unit vectors of the surface coordinate and z, respectively, then incorporating Equation (36) into Equation (35) gives

$$\begin{eqnarray} &&\frac{\varepsilon ^{2}}{c}\frac{\partial I_{0}}{\partial t}+ \left[\mu +\varepsilon (\hat{r}_{b}\cdot \hat{z})\frac{r_{b}}{ct_{f}}\right] \frac{\partial I_{0}}{\partial s}+\varepsilon \left[\mathcal {L}+ \varepsilon \frac{r_{b}}{ct_{f}}\mathcal {L}_{\hat{r}}\right]I_{0}\nonumber \\ &&\quad\quad+\varepsilon \left[(\hat{\Omega }_{0}\cdot \nabla)_{\perp }+ \varepsilon \frac{r_{b}}{ct_{f}}(\hat{r}_{b}\cdot \nabla)_{\perp }\right]I_{0}\nonumber\\ &&\quad\quad+ \sigma _{t,0}I_{0}+\frac{4\varepsilon ^{2}}{ct_{f}}I_{0} = \varepsilon j_{0} \end{eqnarray} \tag{ 37 }$$

and

$$\begin{equation} j_{0}=\left(\frac{\sigma _{t,0}}{\varepsilon }+\varepsilon \sigma _{a,0}\right) \frac{1}{4\pi }\int _{4\pi }I_{0}d\Omega _{0}+\varepsilon q, \end{equation} \tag{ 38 }$$

where q is external or thermal sources, $r_{b}=|\boldsymbol {r}_{b}|$, and the opacities have been assumed isotropic. We have defined $\hat{r}_{b}\cdot \nabla =(\hat{r}_{b}\cdot \hat{z})\partial /\partial z +(\hat{r}_{b}\cdot \nabla)_{\perp }+\mathcal {L}_{\hat{r}}$ as a means of tracking the change with respect to the ballistic fluid trajectories through curved coordinates in a fashion analogous to the streaming term for photon trajectories. At least for spherical symmetry, Equation (36e) implies $\mathcal {L}\rightarrow \varepsilon \mathcal {L}$ in agreement with intuition. Incorporating Equation (38) into Equation (37) and grouping O(ε²) terms on the right-hand side,

$$\begin{eqnarray} &&\left[\mu +\varepsilon (\hat{r}_{b}\cdot \hat{z})\frac{r_{b}}{ct_{f}}\right] \frac{\partial I_{0}}{\partial s}+\varepsilon \mathcal {L}I_{0} +\varepsilon (\hat{\Omega }_{0}\cdot \nabla)_{\perp }I_{0}\nonumber \\ &&\quad\quad+ \sigma _{t,0}I_{0}=\frac{\sigma _{t,0}}{4\pi }\int _{4\pi }I_{0}d\Omega _{0}^{\prime } +O(\varepsilon ^{2}). \end{eqnarray} \tag{ 39 }$$

It is interesting to note that Equation (36) prescribes scalings that make changes in I₀ from curvature and surface variation along the ballistic fluid parcel trajectories at surface b an O(ε²) effect. This is useful, as we only consider O(ε) effects in the construction of the boundary condition. We now neglect surface variations and curvature in the following analysis as these conditions have been thoroughly examined by Malvagi & Pomraning (1991). Following prior authors, we also separate I₀ into a boundary layer solution, I_{0, b}, and an interior solution, I_{0, i} such that I₀ = I_{0, b} + I_{0, i} and lim_{s → ∞}I_{0, b} = 0 (Habetler & Matkowsky 1975; Malvagi & Pomraning 1991). For I_{0, b}, Equation (39) then reduces further to

$$\begin{eqnarray} &&\left[\mu +\varepsilon (\hat{r}_{b}\cdot \hat{z})\frac{r_{b}}{ct_{f}}\right] \frac{\partial I_{0,b}}{\partial s}+\sigma _{t,0}I_{0,b}\nonumber \\ &&\quad\quad=\frac{\sigma _{t,0}}{4\pi }\int _{4\pi }I_{0,b}d\Omega _{0}^{\prime } +O(\varepsilon ^{2}). \end{eqnarray} \tag{ 40 }$$

The boundary and interior solutions may be expanded in the small parameter ε as $I_{0,(b,i)}=\sum _{m=0}^{\infty }I_{0,(b,i)}^{(m)} \varepsilon ^{m}$. Incorporating the ε-expansion into Equation (40), balancing ε⁰ and ε¹ coefficients, and integrating over the azimuthal angle about z, the O(1) and O(ε) equations are

$$\begin{equation} \mu \frac{\partial \tilde{I}_{0,b}^{(0)}}{\partial s}+\sigma _{t,0} \left(\tilde{I}_{0,b}^{(0)}-\frac{1}{2}\int _{-1}^{1}\tilde{I}_{0,b}^{(0)} (\mu ^{\prime })d\mu ^{\prime }\right)=0 \end{equation} \tag{ 41 }$$

$$\begin{eqnarray} &&\mu \frac{\partial \tilde{I}_{0,b}^{(1)}}{\partial s}+\sigma _{t,0} \left(\tilde{I}_{0,b}^{(1)}-\frac{1}{2}\int _{-1}^{1}\tilde{I}_{0,b}^{(1)} (\mu ^{\prime })d\mu ^{\prime }\right)\nonumber \\ &&\quad\quad=-(\hat{r}_{b}\cdot \hat{z}) \frac{r_{b}}{ct_{f}}\frac{\partial \tilde{I}_{0,b}^{(0)}}{\partial s}, \end{eqnarray} \tag{ 42 }$$

where $\tilde{I}_{0,(b,i)}^{(m)}=\int _{0}^{2\pi }I_{0,(b,i)}^{(m)}d\omega$ and ω is the azimuthal angle. Instead of the curvature and spatial variation terms, the heterogeneity of the O() equation is an expansion term. Supposing the boundary intensity at s = 0 is known, matching the asymptotic orders gives (Malvagi & Pomraning 1991)

$$\begin{equation} \tilde{I}_{0,b}^{(0)}(s=0,\mu,t)=F(\boldsymbol {r}_{b},\mu,t)-\frac{1}{2} \phi _{0,i}^{(0)}(\boldsymbol {r}_{b},t) \end{equation} \tag{ 43 }$$

and

$$\begin{equation} \tilde{I}_{0,b}^{(1)}(s=0,\mu,t)=-\frac{1}{2} \left(\phi _{0,i}^{(1)}(\boldsymbol {r}_{b},t)-\frac{\mu }{\sigma _{0,t}} \left.\frac{\partial \phi _{0,i}^{(0)}}{\partial z}\right|_{\boldsymbol {r}_{b}} \right) \end{equation} \tag{ 44 }$$

as boundary conditions, where μ > 0 is the magnitude of the angular projection into the diffusive domain along axis z, $\phi _{0,(b,i)}^{(m)}=\int _{4\pi }I_{0,(b,i)}^{(m)}d\Omega _{0}$ and $F(\boldsymbol {r}_{b},\mu,t)=\int _{0}^{2\pi }I_{0}(\boldsymbol {r}_{b}, \hat{\Omega }_{0},t)d\omega$. Equation (41) along with Equation (43) is a form of the standard half-space albedo problem examined by Case (1960) and Larsen & Habetler (1973). The solution to Equation (41) is (Malvagi & Pomraning 1991)

$$\begin{equation} \tilde{I}_{0,b}^{(0)}=k_{+}^{(0)}+\int _{0}^{1}k^{(0)}(\varpi) \varphi _{\varpi }(\mu)e^{-\sigma _{t,0}s/\varpi }d\varpi, \end{equation} \tag{ 45 }$$

where $k_{+}^{(0)}$ is a constant, ϖ is an eigenvalue of the singular eigenfunction φ_ϖ(μ) (the eigenfunction is formally a distribution), and k⁽⁰⁾(ϖ) is a function determined by the orthogonalities, Equations (48) and (49) below (Habetler & Matkowsky 1975). The constant, $k_{+}^{(0)}$, constitutes the "discrete" component of the solution (Habetler & Matkowsky 1975) and is the limiting behavior of a more general discrete eigenvalue solution having taken σ_{0, a} ≪ σ_{0, s} with ε. The distribution, φ_ϖ(μ), is given by

$$\begin{equation} \varphi _{\varpi }(\mu)=\frac{\varpi }{2}P:\frac{1}{\varpi -\mu } +\lambda (\varpi)\delta (\varpi -\mu), \end{equation} \tag{ 46 }$$

where δ(ϖ − μ) is the Dirac delta function and λ(ϖ) = 1 − ϖtanh ⁻¹(ϖ) ensures a normalization of $\int _{-1}^{1}\varphi _{\varpi }(\mu)d\mu =1$ for all ϖ ∈ [0, 1]. The P: is merely a notational device to indicate the principal value is taken upon integration (Habetler & Matkowsky 1975; Case 1960). Case (1960) rigorously proves that a function, H(μ), may be found such that

$$\begin{equation} \int _{0}^{1}\mu H(\mu)\varphi _{\varpi }(\mu)d\mu = 0 \end{equation} \tag{ 47 }$$

when $\tilde{I}_{0,b}^{(0)}$ satisfies a Hölder condition. It turns out H(μ) is Chandrasekhar's H-function (Malvagi & Pomraning 1991). Making use of the Poincaré–Bertrand formula (Hang & Jiang 2009) and Equation (47), the orthogonalities (Malvagi & Pomraning 1991)

$$\begin{equation} \int _{-1}^{1}\mu \varphi _{\varpi ^{\prime }}(\mu)\varphi _{\varpi }(\mu)d\mu = N(\varpi)\delta (\varpi -\varpi ^{\prime }) \end{equation} \tag{ 48 }$$

and

$$\begin{equation} \int _{0}^{1}\mu H(\mu)\varphi _{\varpi ^{\prime }}(\mu)\varphi _{\varpi }(\mu) d\mu = N(\varpi)H(\varpi)\delta (\varpi -\varpi ^{\prime }), \end{equation} \tag{ 49 }$$

where N(ϖ) = ϖ(λ(ϖ)² + (πϖ/2)²), must hold. Incorporating Equation (46) into Equation (45) and Equation (45) into the boundary condition Equation (43), Equations (47) and (49) indicate that multiplying the result by μH(μ) or $\mu H(\mu)\varphi _{\varpi ^{\prime }}(\mu)$ and integrating over μ ∈ [0, 1] yields

$$\begin{equation} k_{+}^{(0)}=\frac{\sqrt{3}}{2}\int _{0}^{1}\mu H(\mu)F(\boldsymbol {r}_{b}, \mu)d\mu - \frac{1}{2}\phi _{0,i}(\boldsymbol {r}_{b}) \end{equation} \tag{ 50 }$$

or

$$\begin{equation} k^{(0)}(\varpi)=\frac{1}{N(\varpi)H(\varpi)}\int _{0}^{1}\mu H(\mu) \varphi _{\varpi }(\mu)F(\boldsymbol {r}_{b},\mu)d\mu, \end{equation} \tag{ 51 }$$

respectively (Malvagi & Pomraning 1991; $\int _{0}^{1}\mu H(\mu)=2/\sqrt{3}$ and we have dropped t as an argument). As s → ∞, Equation (42) becomes homogeneous; then as s → ∞, the O(ε) solution $\tilde{I}_{0,b}^{(1)}$ must tend to some constant, $k_{+}^{(1)}$, as well. Now using the boundary condition (45), $k_{+}^{(1)}$ may be found in a similar manner to $k_{+}^{(0)}$ as

$$\begin{eqnarray} k_{+}^{(1)} &=&\frac{\sqrt{3}}{2}\int _{0}^{1}\mu H(\mu) \left[-\frac{1}{2}\left(\phi _{0,i}^{(1)}(\boldsymbol {r}_{b})- \frac{\mu }{\sigma _{0,t}}\left.\frac{\partial \phi _{0,i}^{(0)}}{\partial z} \right|_{\boldsymbol {r}_{b}}\right)\right]d\mu \nonumber \\ &&-\int _{0}^{\infty }\int _{-1}^{1}\beta (s,\mu)(\hat{r}_{b}\cdot \hat{z}) \frac{r_{b}}{ct_{f}}\frac{\partial \tilde{I}_{0,b}^{(0)}}{\partial s}d\mu \,ds, \end{eqnarray} \tag{ 52 }$$

where a considerable variational analysis by Malvagi & Pomraning (1991) gives

$$\begin{equation} \beta (s,\mu)=1+\frac{3}{2}(\mu +\sigma _{0,t}s)-\left(1+\frac{3}{2}\mu \right)e^{-\sigma _{0,t}s/\mu }\Theta (-\mu), \end{equation} \tag{ 53 }$$

and Θ is the Heaviside function. Roughly speaking, to obtain Equation (53), Malvagi & Pomraning (1991) construct a linear functional for $k_{+}^{(1)}$ with Lagrange multipliers as an estimate to $k_{+}^{(1)}$, find an adjoint transport solution (Lewis & Miller 1993, p. 47) and μ times the adjoint solution as the appropriate multipliers, and incorporate constants as trial functions for $\tilde{I}_{b}^{(1)}$ and the adjoint solution.

Now we may use the boundary layer constraint, $\lim _{s\rightarrow \infty }\tilde{I}_{0,b}=0$; $k_{+}^{0}$ and $k_{+}^{(1)}$ must vanish. Equations (52) becomes (Malvagi & Pomraning 1991)

$$\begin{eqnarray} &&\phi _{0,i}^{(1)}(\boldsymbol {r}_{b})+\frac{\lambda }{\sigma _{t}}\hat{z}\cdot \nabla \phi _{0,i}^{(0)}=2(\hat{r}_{b}\cdot \hat{z})\frac{r_{b}}{ct_{f}} \int _{0}^{\infty }\ldots\left(\int _{-1}^{1}\beta (\eta,\mu)\right. \nonumber \\ &&\left. \quad\times\int _{0}^{1} \frac{1}{\varpi }k(\varpi)\varphi _{\varpi }(\mu)e^{-\eta /\varpi } d\varpi \,d\mu \right)\,d\eta,\quad \end{eqnarray} \tag{ 54 }$$

where $\lambda =\sqrt{3}\int _{0}^{1}\mu ^{2}H(\mu)d\mu /2$ and σ_{t, 0}s = η. We have incorporated the eigenvalue form of $\tilde{I}_{0,b}^{(0)}$ into Equation (54). Multiplying Equation (54) by ε, adding the result to $\phi _{0,i}^{(0)}=\sqrt{3}\int _{0}^{1}\mu H(\mu) F(\boldsymbol {r}_{b},\mu)d\mu$, reintroducing the interior solution $\phi _{0,i}=\phi _{0,i}^{(0)}+\varepsilon \phi _{0,i}^{(1)}+ O(\varepsilon ^{2})$, and reverting the ε-scalings from Equation (36) gives

$$\begin{eqnarray} &&\phi _{0,i}(\boldsymbol {r}_{b},t)+\frac{\lambda }{\sigma _{t,0}}\hat{z}\cdot \nabla \phi _{0,i}=2(\hat{r}_{b}\cdot \hat{z})\frac{r_{b}}{ct_{f}} \int _{0}^{\infty }\ldots \nonumber \\ &&\quad\times\left(\int _{-1}^{1}\beta (\eta,\mu)\int _{0}^{1} \frac{1}{\varpi }k(\varpi)\varphi _{\varpi }(\mu)e^{-\eta /\varpi } d\varpi \,d\mu \right)\,d\eta \nonumber \\ &&\quad+\sqrt{3}\int _{0}^{1}\mu H(\mu)F(\boldsymbol {r}_{b},\mu,t)d\mu, \end{eqnarray} \tag{ 55 }$$

which should be correct to O(ε²). We find

$$\begin{eqnarray} &&\int _{0}^{\infty }\int _{-1}^{1}\beta (\eta,\mu)\varphi _{\varpi }(\mu) e^{-\eta /\varpi }d\mu d\eta \nonumber \\ &&\quad=\varpi \left(1+\frac{3}{2}\varpi \right) -\frac{\varpi ^{2}}{2}(1+3\varpi)\ln \left(\frac{1+\varpi }{\varpi }\right) \nonumber \\ &&\quad\quad+\frac{\varpi ^{2}}{2(1+\varpi)}\left(\frac{5}{2}+3\varpi \right)\nonumber \\ &&\quad\equiv \varpi h(\varpi), \end{eqnarray} \tag{ 56 }$$

so

$$\begin{eqnarray} &&\int _{0}^{1}\frac{k^{(0)}(\varpi)}{\varpi } \int _{0}^{\infty }\int _{-1}^{1}\beta (\eta,\mu)\varphi _{\varpi }(\mu) e^{-\eta /\varpi }\,d\mu \,d\eta \,d\varpi \nonumber \\\ms\ms &&\quad\quad=\int _{0}^{1}h(\varpi)k^{(0)}(\varpi)d\varpi. \end{eqnarray} \tag{ 57 }$$

Incorporating Equation (51), using $\sqrt{3}\mu H(\mu)\approx 2W(\mu)$ which is a corollary of the variational derivation of β(s, μ) by Malvagi & Pomraning (1991), and defining

$$\begin{equation} G_{U}(\mu)=1+\frac{2}{c}\hat{z}\cdot \boldsymbol {U}(\boldsymbol {r}_{b},t_{f}) \int _{0}^{1}\frac{h(\varpi)\varphi _{\varpi }(\mu)}{(2+3\varpi)N(\varpi)} d\varpi, \end{equation} \tag{ 58 }$$

Equation (55) becomes

$$\begin{equation} \phi _{0,i}(\boldsymbol {r}_{b},t)+\frac{\lambda }{\sigma _{t,0}}\hat{z} \nabla \phi _{0,i}= 2\int _{0}^{1}W(\mu)G_{U}(\mu)F(\boldsymbol {r}_{b},\mu,t)d\mu. \end{equation} \tag{ 59 }$$

Notwithstanding the factor G_U, the form of Equation (59) is fortunately similar to Equation (33). To evaluate the integral in Equation (59), we incorporate Equation (46) for φ_ϖ(μ) to obtain

$$\begin{eqnarray} &&\int _{0}^{1}\frac{h(\varpi)\varphi _{\varpi }(\mu)}{(2+3\varpi)N(\varpi)}d\varpi = \lambda (\mu)\frac{h(\mu)}{(2+3\mu)N(\mu)}\nonumber \\\ms\ms &&\quad\quad+\frac{1}{2}\int _{0}^{1}\left(\varpi \frac{h(\varpi)}{(2+3\varpi)N(\varpi)}- \mu \frac{h(\mu)}{(2+3\mu)N(\mu)}\right)\nonumber\\\ms\ms &&\quad\quad\times\frac{d\varpi }{(\varpi -\mu)}+ \frac{\mu }{2}\frac{h(\mu)}{(2+3\mu)N(\mu)}\ln \left(\frac{1-\mu }{\mu }\right), \end{eqnarray} \tag{ 60 }$$

where, following Malvagi & Pomraning (1991), we have converted the principal value integration to a nonsingular form amenable to quadrature.

Equation (60) along with the form of the functions λ(μ), N(μ), and h(μ) reveal the leading order behavior of the angular dependence in G_U. For the case of μ → 0: λ(μ) → 1, N(μ) → μ and h(μ) → 1. Hence the first term on the right-hand side of Equation (60) tends to 0.5/μ as μ → 0. From the last term on the right-hand side of Equation (60), the next order of divergence as μ → 0 is logarithmic. The remaining terms tend to a bounded integral over ϖ as μ → 0. We find that the behavior in μ of Equation (60) is well approximated by C₁/μ − C₂μ where C₁ and C₂ are positive constants. In Section 3.3, we test $0.5c(G_{U}-1)/(\hat{z}\cdot \boldsymbol {U})=0.55/\mu -1.25\mu$ for a three mean free path threshold between IMC and DDMC and $0.5c(G_{U}-1)/(\hat{z}\cdot \boldsymbol {U})=0.6/\mu -1.25\mu$ for a ten mean free path threshold between IMC and DDMC. These choices are seen to be good approximations to the quadrature expressed in Equation (60).

Finally, in Equation (59) we replace $\hat{z}$ with $\boldsymbol {n}$, drop the subscript i, replace μ with $|\hat{\Omega }_{0}\cdot \boldsymbol {n}|$, and extend to multigroup to obtain

$$\begin{eqnarray} &&\phi _{0,g}(\boldsymbol {r}_{b},t)+\left(\frac{\lambda }{\sigma _{a,g,n}+\sigma _{s,g,n}}\right)\boldsymbol {n}\cdot \nabla \phi _{0,g}(\boldsymbol {r}_{b},t)\nonumber \\\ms\ms &&\quad=2\int _{\hat{\Omega }_{0}\cdot \boldsymbol {n}<0}\int _{\nu _{g+1/2}}^{\nu _{g-1/2}} W(|\hat{\Omega }_{0}\cdot \boldsymbol {n}|)G_{U}(|\hat{\Omega }_{0}\cdot \boldsymbol {n}|) \ldots \nonumber \\\ms\ms &&\quad\quad\times I_{0,\nu _{0}}(\boldsymbol {r}_{b},\hat{\Omega }_{0},t)d\nu _{0}d\Omega _{0}. \end{eqnarray} \tag{ 61 }$$

when $\boldsymbol {U}=0$, the standard diffusion limit boundary condition Equation (33) is obtained. As is seen in Section 3.3, some form of the new factor, G_U, must be implemented at IMC–DDMC boundaries to furnish satisfactory agreement between the hybrid method and pure IMC over a grid moving at relativistic speed.

We now turn to the discussion of hybridizing IMC and DDMC in velocity space (cells) and frequency space (groups). To incorporate Equation (61) into the DDMC routine, the second term on the left-hand side can be finite differenced and incorporated into the discretized form of Equation (29). The resulting expression for $\boldsymbol {F}_{0,g}$ may then be incorporated into the discrete form of Equation (24) to yield

$$\begin{eqnarray} &&\frac{1}{c}\frac{\partial \phi _{0,j,g}}{\partial t}+\left [\sum _{j^{\prime \prime }\not=j^{\prime }} \sigma _{j\rightarrow j,^{\prime \prime }g}+\sigma _{b(j,j^{\prime }),g}+\sigma _{s,j,g,n}+f_{j,n}\sigma _{a,j,g,n}\right. \nonumber &&\left. \quad\quad+(1-\gamma _{j,g,n})(1-f_{j,n})\sigma _{a,j,g,n}\right ]\phi _{0,j,g}\nonumber\\\ms\ms &&\quad= f_{j,n}\gamma _{j,g,n}\sigma _{P,j,n}acT_{j,n}^{4}+\frac{1}{V_{j}}\sum _{j^{\prime \prime }\not=j}V_{j^{\prime \prime }}\sigma _{j^{\prime \prime }\rightarrow j,g}\phi _{0,j,^{\prime \prime }g}\nonumber \\\ms\ms &&\quad\quad+ \gamma _{j,g,n}(1-f_{j,n})\sum _{g^{\prime }\not=g}^{G_{j}}\sigma _{a,j,g^{\prime },n}\phi _{0,j,g^{\prime }}\nonumber \\\ms\ms &&\quad\quad+\frac{1}{V_{j}}\int _{A_{b(j,j^{\prime })}}\int _{\hat{\Omega }_{0}\cdot \boldsymbol {n}<0} \int _{\nu _{g+1/2}}^{\nu _{g-1/2}}\ldots \nonumber \\\ms\ms &&\quad\quad\times P_{b(j,j^{\prime })}(|\hat{\Omega }_{0}\cdot \boldsymbol {n}|) |\hat{\Omega }_{0}\cdot \boldsymbol {n}|I_{0,g}d\nu _{0}d\Omega _{0}d^{2}\boldsymbol {r}\nonumber \\\ms\ms &&\quad\quad+\sum _{g^{\prime }=1}^{G_{j}}\sigma _{s,j,n}(g^{\prime }\rightarrow g)\phi _{0,j,g^{\prime }}, \end{eqnarray} \tag{ 62 }$$

where the b(j, j') subscript indicates the boundary between DDMC cell j and IMC cell j', $\sigma _{b(j,j^{\prime }),g}$ is the modified leakage opacity resulting from Equation (61), j'' denotes cells adjacent to j with a diffusive group g, $A_{b(j,j^{\prime })}$ is the area of the boundary between cells j and j', and $P_{b(j,j^{\prime })}$ may be interpreted as a transmission probability for an IMC particle in cell j' incident on cell j (Densmore et al. 2007). We have only converted one cell-group adjacent to (j, g) to IMC in Equation (62) but in principle any of the leakage opacities and sources could be replaced by $\sigma _{b(j,j^{\prime \prime })}$ and IMC boundary transmission sources. In spherically symmetric geometry over a velocity grid, the forms of $\sigma _{j\rightarrow j,^{\prime \prime }g}$, $\sigma _{b(j,j^{\prime }),g}$, and $P_{b(j,j^{\prime })}$ are (Abdikamalov et al. 2012)

$$\begin{eqnarray} &&\sigma _{j\rightarrow j,^{\prime \prime }g}=\nonumber \\ &&\left\lbrace \begin{array}{@{}l}\displaystyle \frac{2U_{j-1/2}^{2}}{(t_{n}+t_{\min })^{2}\Delta (U^{3})_{j}} \left(\frac{1}{\sigma _{j-1/2,g,n}^{-}\Delta U_{j-1}+ \sigma _{j-1/2,g,n}^{+}\Delta U_{j}}\right) \\\ms\ms \,,\;\;j^{\prime \prime }=j-1\\\ms\ms \displaystyle \frac{2U_{j+1/2}^{2}}{(t_{n}+t_{\min })^{2}\Delta (U^{3})_{j}} \left(\frac{1}{\sigma _{j+1/2,g,n}^{-}\Delta U_{j}+ \sigma _{j+1/2,g,n}^{+}\Delta U_{j+1}}\right) \\\ms\ms \,,\;\;j^{\prime \prime }=j+1 \end{array},\right. \nonumber \\ \end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray} &&\sigma _{b(j,j^{\prime }),g}=\nonumber \\ &&\left\lbrace \begin{array}{@{}l}\displaystyle \frac{2U_{j-1/2}^{2}}{(t_{n}+t_{\min })\Delta (U^{3})_{j}} \left(\frac{1}{\sigma _{j-1/2,g,n}^{+}(t_{n}+t_{\min })\Delta U_{j}+ 2\lambda }\right),\nonumber \\\ms\ms \,,\;\;j^{\prime }=j-1\\\ms\ms \displaystyle \frac{2U_{j+1/2}^{2}}{(t_{n}+t_{\min })\Delta (U^{3})_{j}} \left(\frac{1}{\sigma _{j+1/2,g,n}^{-}(t_{n}+t_{\min })\Delta U_{j}+ 2\lambda }\right), \\\ms\ms \,,\;\;j^{\prime }=j+1 \end{array},\right. \nonumber \\ \end{eqnarray} \tag{ 64 }$$

and

$$\begin{equation} P_{b(j,j^{\prime }),g}(\mu)= \left\lbrace \begin{array}{@{}l}\displaystyle \frac{4(1+3\mu /2)G_{U_{j-1/2}}(\mu)}{3\sigma _{j-1/2,g,n}^{+}(t_{n}+t_{\min })\Delta U_{j}+6\lambda }\\\ms \,,\;\;j^{\prime }=j-1\\\ms \displaystyle \frac{4(1+3\mu /2)G_{U_{j+1/2}}(\mu)}{3\sigma _{j+1/2,g,n}^{-}(t_{n}+t_{\min })\Delta U_{j}+6\lambda }\\\ms \,,\;\;j^{\prime }=j+1 \end{array},\right. \end{equation} \tag{ 65 }$$

respectively, where ΔU_j = U_{j + 1/2} − U_{j − 1/2} is the radial velocity width, $\Delta (U^{3})_{j}=U_{j+1/2}^{3}-U_{j-1/2}^{3}$. Following Densmore et al. (2007), the σ_{j ± 1/2, g, n} = σ_{j ± 1/2, a, g, n} + σ_{j ± 1/2, s, g, n} are evaluated with cell-edge temperature while the + (−) superscript indicates the remaining material properties are evaluated on the outer (inner) side of the cell edge with respect to index j. The asymptotic diffusion-limit W(μ) has been used to determine $P_{b(j,j^{\prime }),g}(\mu)$. It can be shown that the expressions in Equations (63)–(65) tend to the appropriate planar geometry forms when the inner radius of cell j is large with respect to the width of cell j.

Since DDMC tracks comoving radiation energy, the interface conditions are in the IMC–DDMC comoving frame as well. Spatial hybridization of IMC and DDMC in the comoving frame is evident from Equation (62). With a hybridized IMC–DDMC method, the DDMC particles are advected along with the material while the IMC particles are not. Since an IMC particle may be moved over the velocity grid past cell bounds, a DDMC zone may eventually advect into it. In our outflow tests, we find that the best treatment of this occurrence is to stop the IMC particle at the DDMC cell boundary and allow the diffuse albedo condition described by Equation (65) to determine the particle's admission into the DDMC region. So the split velocity position shift algorithm for IMC particle p is now:

• 1.  
  Find current fluid cell j from $\boldsymbol {U}_{p,{\rm old\; cen}}$ and expected fluid cell j' from $\boldsymbol {U}_{p,{\rm * cen}} (t_{n}+\alpha _{2}\Delta t_{n}+t_{\min }) = \boldsymbol {U}_{p,{\rm old\; cen}}(t_{n}+t_{\min })$.
• 2.  
  IMC–DDMC: $\boldsymbol {U}_{p,{\rm * cen}}\rightarrow \boldsymbol {U}_{p,*}$
• 3.  
  If p has become a DDMC particle, the position update is finished.
• 4.  
  Otherwise, use $\boldsymbol {U}_{p,{\rm new\; cen}}(t_{n}+\Delta t_{n}+t_{\min }) = \boldsymbol {U}_{p,*}(t_{n}+\alpha _{2}\Delta t_{n}+t_{\min })$ to find fluid cell j'.

Having discussed hybridization of IMC with DDMC at velocity grid boundaries, it remains to discuss the treatment of hybrid scattering and the effect non-uniform group structuring at grid boundaries. The hybrid scattering follows simply from both Abdikamalov and Densmore's multifrequency IMC–DDMC methods. We may replace one of the DDMC scattering group source terms in Equation (31) or Equation (62) with an IMC scattering source:

$$\begin{eqnarray} &&\gamma _{j,g,n}(1-f_{j,n})\sigma _{a,j,g^{\prime },n}\phi _{0,j,g^{\prime }}\rightarrow\frac{\gamma _{j,g,n}(1-f_{j,n})}{V_{j}} \nonumber \\ &&\quad\quad\times\int _{V_{j}}\int _{4\pi } \int _{\nu _{g^{\prime }+1/2}}^{\nu _{g^{\prime }-1/2}}\sigma _{a,j,n}I_{0,\nu _{0}}d\nu _{0}d\Omega _{0}d^{3}\boldsymbol {r} \end{eqnarray} \tag{ 66 }$$

for effective scattering and

$$\begin{eqnarray} &&\sigma _{s,j,n}(g^{\prime }\rightarrow g)\phi _{0,j,g^{\prime }} \rightarrow\frac{1}{V_{j}}\int _{V_{j}}\int _{4\pi }\int _{4\pi }\int _{\nu _{g+1/2}}^{\nu _{g-1/2}} \int _{\nu _{g^{\prime }+1/2}}^{\nu _{g^{\prime }-1/2}}\frac{\nu _{0}}{\nu _{0}^{\prime }}\ldots \nonumber \\ &&\quad\times\sigma _{s,j,n}(\nu _{0}^{\prime }\rightarrow \nu _{0},\hat{\Omega }_{0}^{\prime }\cdot \hat{\Omega }_{0}) I_{0,\nu _{0}^{\prime }}d\nu _{0}^{\prime }d\nu _{0}d\Omega _{0}^{\prime }d\Omega _{0}d^{3}\boldsymbol {r} \end{eqnarray} \tag{ 67 }$$

for physical scattering. For the non-uniform group structuring at grid boundaries, we extend the Planck averaging formula presented by Densmore et al. (2012) to

$$\begin{eqnarray} \tilde{\sigma }_{j\rightarrow j^{\prime },g} &=& \left(\frac{\sum _{g_{D}^{\prime }}b_{j,g\leftrightarrow g_{D}^{\prime },n}}{b_{j,g,n}}\right) \sigma _{j\rightarrow j^{\prime },g}\nonumber \\ &&+\left(\frac{\sum _{g_{T}^{\prime }}b_{j,g\leftrightarrow g_{T}^{\prime },n}}{b_{j,g,n}}\right) \sigma _{b(j,j^{\prime }),g} \end{eqnarray} \tag{ 68 }$$

where $g_{D}^{\prime }$ is a DDMC group index in cell j', $g_{T}^{\prime }$ is an IMC group index in cell j', $b_{j,g,n}=\int _{\nu _{g+1/2}}^{\nu _{g-1/2}}b_{0,\nu _{0}}(T_{j,n})d\nu _{0}$, and

$$\begin{eqnarray} &&b_{j,g\leftrightarrow g^{\prime },n}=\nonumber \\ &&\left\lbrace \begin{array}{@{}l}\displaystyle \int _{\min ([\nu _{g+1/2},\nu _{g-1/2}]\cap [\nu _{g^{\prime }+1/2},\nu _{g^{\prime }-1/2}])} ^{\max ([\nu _{g+1/2},\nu _{g-1/2}]\cap [\nu _{g^{\prime }+1/2},\nu _{g^{\prime }-1/2}])}b_{0,\nu _{0}}(T_{j,n})d\nu _{0},\\\ms \;\;[\nu _{g+1/2},\nu _{g-1/2}]\cap [\nu _{g^{\prime }+1/2},\nu _{g^{\prime }-1/2}]\not=\emptyset \\\ms 0\,,\;\;[\nu _{g+1/2},\nu _{g-1/2}]\cap [\nu _{g^{\prime }+1/2},\nu _{g^{\prime }-1/2}]=\emptyset. \end{array}\right. \end{eqnarray} \tag{ 69 }$$

Equation (68) presumes the radiation field near the cell boundary is well approximated by a Planck function at the diffusive cell temperature. With $\sigma _{b(j,j^{\prime }),g}$ on the right-hand side of Equation (68), it is clear that changing $\sigma _{j\rightarrow j^{\prime },g}\rightarrow \tilde{\sigma }_{j\rightarrow j^{\prime },g}$ requires a complimentary modification to the right-hand side of our DDMC equation that balances the field at boundaries (Densmore et al. 2012).

In Figure 1, group bounds are given by horizontal lines at different values along the vertical frequency axis and cells are enumerated along the horizontal axis. The circles are possible (j, g) locations of IMC–DDMC particles and the dashed lines represent transition possibilities for particles at those phase space locations. If all group bounds are aligned, the method reduces in complexity to the standard multigroup approach and the set of (j, g)-transition graphs become completely reducible. If a leakage is sampled, the probability of a DDMC to DDMC leakage transition is $(b_{j,g\leftrightarrow g_{D}^{\prime },n}\sigma _{j\rightarrow j^{\prime },g}) /(b_{j,g,n}\tilde{\sigma }_{j\rightarrow j^{\prime },g})$ and the probability of an DDMC to IMC leakage transition is $(b_{j,g\leftrightarrow g_{T}^{\prime },n}\sigma _{b(j,j^{\prime }),g}) /(b_{j,g,n}\tilde{\sigma }_{j\rightarrow j^{\prime },g})$. So the transition values for DDMC are these probabilities that sum to 1 for each (j, g) DDMC region. If (j, g) is treated with IMC, then the group in the subsequent cell is determined by particle frequency in the comoving frame of the boundary.

Zoom In Zoom Out Reset image size

The fully hybridized comoving DDMC equation with elastic physical scattering is

$$\begin{eqnarray} &&\frac{1}{c}\frac{\partial \phi _{0,j,g}}{\partial t}+\left (\sum _{j^{\prime }}\tilde{\sigma }_{j\rightarrow j^{\prime },g}+(1-\gamma _{j,g,n})(1-f_{j,n})\sigma _{a,j,g,n}\nonumber \right. \\ &&\left. \quad\quad+f_{j,n}\sigma _{a,j,g,n}\right)\phi _{0,j,g}=f_{j,n}\gamma _{j,g,n}\sigma _{P,j,n}acT_{j,n}^{4}\nonumber \\ &&\quad\quad+\frac{1}{V_{j}}\sum _{j^{\prime }}V_{j^{\prime }}\sum _{g_{D}^{\prime }} \frac{b_{j,g\leftrightarrow g_{D}^{\prime },n}}{b_{j,g,n}}\sigma _{j^{\prime }\rightarrow j,g_{D}^{\prime }}\phi _{0,j,g_{D}^{\prime }}\nonumber \\ &&\quad\quad+\frac{1}{V_{j}}\sum _{j^{\prime }}\sum _{g_{T}^{\prime }}\int _{A_{b(j,j^{\prime })}}\int _{\hat{\Omega }_{0}\cdot \boldsymbol {n}<0} \int _{g\leftrightarrow g_{T}^{\prime }}\ldots \nonumber \\ &&\quad\quad\times P_{b(j,j^{\prime })}(|\hat{\Omega }_{0}\cdot \boldsymbol {n}|)|\hat{\Omega }_{0}\cdot \boldsymbol {n}|I_{0,\nu _{0}}d\nu _{0}d\Omega _{0}d^{2}\boldsymbol {r}\nonumber \\ &&\quad\quad+\frac{\gamma _{j,g,n}(1-f_{j,n})}{V_{j}}\ldots \nonumber \\ &&\quad\quad\times\sum _{g_{T}}\int _{V_{j}}\int _{4\pi }\int _{\nu _{g_{T}+1/2}}^{\nu _{g_{T}-1/2}} \sigma _{a,j,n}I_{0,\nu _{0}}d\nu _{0}d\Omega _{0}d^{3}\boldsymbol {r}\nonumber \\ &&\quad\quad+\gamma _{j,g,n}(1-f_{j,n})\sum _{g_{D}}\sigma _{a,j,g_{D},n}\phi _{0,j,g_{D}}, \end{eqnarray} \tag{ 70 }$$

where

$$\begin{eqnarray} &&\int _{g\leftrightarrow g^{\prime }}f(\nu)d\nu =\nonumber \\ &&\left\lbrace \begin{array}{@{}l}\displaystyle \int _{\min ([\nu _{g+1/2},\nu _{g-1/2}]\cap [\nu _{g^{\prime }+1/2},\nu _{g^{\prime }-1/2}])} ^{\max ([\nu _{g+1/2},\nu _{g-1/2}]\cap [\nu _{g^{\prime }+1/2},\nu _{g^{\prime }-1/2}])}f(\nu)d\nu,\\\ms \;\;[\nu _{g+1/2},\nu _{g-1/2}]\cap [\nu _{g^{\prime }+1/2},\nu _{g^{\prime }-1/2}]\not=\emptyset \\\ms 0,\;\;[\nu _{g+1/2},\nu _{g-1/2}]\cap [\nu _{g^{\prime }+1/2},\nu _{g^{\prime }-1/2}]=\emptyset. \end{array}\right. \end{eqnarray} \tag{ 71 }$$

Groups represented by sub-indices g_D or g_T can be thought of as the complimentary sets over frequency that are used during the advance of particles. In other words, an optimization over frequency may be applied such that the comoving group structure in a time step is distinct from the user-prescribed group structure. This regrouping must not degrade the accuracy of the observables of interest (for instance, a spectral tally of particles leaving the domain). Densmore's approach to multifrequency DDMC lumps adjacent groups into a gray DDMC group (Densmore et al. 2012). Indices (j, g) and (j, g + 1) are combined if their opacities multiplied by a characteristic cell length are greater than a threshold number of mean free paths per cell (Densmore et al. 2012; Abdikamalov et al. 2012). Additionally, this combined group has distinct Planck and Rosseland opacities,

$$\begin{equation} \sigma _{P,j,g\cup g+1,n}=\frac{b_{j,g,n}\sigma _{a,j,g,n}+b_{j,g+1,n}\sigma _{a,j,g+1,n}}{b_{j,g,n}+b_{j,g+1,n}}, \end{equation} \tag{ 72 }$$

and

$$\begin{equation} \frac{1}{\sigma _{R,j,g\cup g+1,n}}=\frac{b_{j,g,n}/\sigma _{a,j,g,n}+b_{j,g+1,n}/\sigma _{a,j,g+1,n}}{b_{j,g,n}+b_{j,g+1,n}}, \end{equation} \tag{ 73 }$$

respectively, where g∪g + 1 denotes a DDMC group union. Equations (72) and (73) can be used in place of σ_{a, j, g, n} in the formulae for leakage opacity and in absorption sampling. The increase in efficiency is obtained from

$$\begin{equation} \gamma _{j,g\cup g+1,n}=\gamma _{j,g,n}+\gamma _{j,g+1,n}. \end{equation} \tag{ 74 }$$

For either g or g', γ_{j, g∪g + 1, n} used in the DDMC effective out-scattering coefficient, (1 − f_{j, n})(1 − γ_{j, g∪g + 1, n}), in place of the values on the left-hand side of Equation (74). For Densmore's simulations (Densmore et al. 2012), only opacities that are monotonic in frequency are considered, so the group lumping method is conflated with a threshold frequency. Nothing in Densmore's asymptotic analysis suggests that group lumping cannot be employed in frequency regions away from ν = 0. Densmore indeed notes the possibility of this extension (Densmore et al. 2012).

Our first verification is in a static material with a multifrequency structure that is amenable to analytic solution. Specifically, we use the Su & Olson (1999) picket-fence opacity structure. We apply this picket-fence opacity distribution to the static P₁ equations with thermal coupling. McClarren et al. (2008a) have solved the thermally coupled P₁ equations in gray materials for several 1D geometries. The P₁ method is higher order than diffusion, but still approximate. This verification is therefore performed in an optically thick material where DDMC, IMC, and P₁ should show quantitative agreement for a range of spatial grids.

The picket-fence opacity dependence on frequency can be constructed as the limit of a discrete multigroup distribution. Partitioning the frequency grid into regular Δν intervals, a portion p₁ of Δν is attributed an opacity σ₁ while the remainder p₂ = 1 − p₁ is attributed an opacity σ₂. The limit as Δν → 0 of this alternating grouping is the picket-fence opacity σ(ν). Each picket at some ν is dense over the real number line of ν; meaning both picket values are in any nonempty, open interval over the real number line for ν.

The picket-fence distribution has the nice property that integrals of I_ν and B_ν over the dense groupings simplify when these integrands are assumed to be smooth (Su & Olson 1999). Su and Olson solve the transport equations with the picket-fence opacity to obtain a semi-analytic result. For specific values of σ₁, p₁, σ₂, and p₂, we may develop a simple generalization of McClarren's P₁ solution that includes a rudimentary test of multifrequency for SuperNu.

Neglecting scattering, taking the zeroth and first angular moments of Equation (1), and integrating the result over the set of frequencies that only yield a contribution from picket g ∈ {1, 2}, the thermal picket-fence P₁ equations in planar 1D geometry are

$$\begin{equation} \frac{1}{c}\left(\frac{\partial E_{g}}{\partial t}+\frac{\partial F_{g}}{\partial z}\right) =\sigma _{g}(p_{g}aT^{4}-E_{g})+p_{g}S, \end{equation} \tag{ 75 }$$

$$\begin{equation} \frac{1}{c}\frac{\partial F_{g}}{\partial t}+\frac{c}{3}\frac{\partial E_{g}}{\partial z} =-\sigma _{g}F_{g}, \end{equation} \tag{ 76 }$$

$$\begin{equation} \frac{C_{v}(T)}{c}\frac{\partial T}{\partial t}=\sum _{g^{\prime }=1}^{2}\sigma _{g}E_{g}-\bar{\sigma } aT^{4}, \end{equation} \tag{ 77 }$$

where z is the spatial coordinate, E_g = ϕ_g/c is the radiation energy density, $\bar{\sigma }=p_{1}\sigma _{1}+p_{2}\sigma _{2}$, and S is an external source. The system of equations is linearized with C_v(T) = a'T³ (McClarren et al. 2008a; Su & Olson 1999). Furthermore, we apply the usual non-dimensionalizations for convenience (McClarren et al. 2008a; Su & Olson 1999):

$$\begin{equation} x=\bar{\sigma }z,\;\;\epsilon =\frac{4a}{a^{\prime }},\;\;\tau =\epsilon c\bar{\sigma }t,\;\; w_{g}=\frac{\sigma _{g}}{\bar{\sigma }} \end{equation} \tag{ 78 }$$

and

$$\begin{equation} \mathcal {E}_{g}=\frac{E_{g}}{aT_{r}^{4}},\;\;\mathcal {F}_{g}=\frac{F_{g}}{aT_{r}^{4}}, \;\;\mathcal {M}=\frac{T^{4}}{T_{r}^{4}},\;\;Q=\frac{S}{\bar{\sigma }aT_{r}^{4}}, \end{equation} \tag{ 79 }$$

where T_r is a reference temperature. Incorporating Equations (78) and (79), Equations (75)–(77) become

$$\begin{equation} \epsilon \frac{\partial \mathcal {E}_{g}}{\partial \tau }+\frac{1}{c} \frac{\partial \mathcal {F}_{g}}{\partial x}=w_{g}(p_{g}\mathcal {M}-\mathcal {E}_{g})+p_{g}Q, \end{equation} \tag{ 80 }$$

$$\begin{equation} \epsilon \frac{\partial \mathcal {F}_{g}}{\partial \tau }+\frac{c}{3}\frac{\partial \mathcal {E}_{g}}{\partial x}=-w_{g}\mathcal {F}_{g}, \end{equation} \tag{ 81 }$$

$$\begin{equation} \frac{\partial \mathcal {M}}{\partial \tau }=\sum _{g^{\prime }=1}^{2}w_{g}\mathcal {E}_{g}-\mathcal {M}. \end{equation} \tag{ 82 }$$

If Equations (80)–(82) are Laplace transformed over time and reduced to equations for only energy density, what results are two fourth-order linear differential equations in space with coefficients that are algebraically irrational functions of the Laplace variable, s. Denoting transformed quantities with a tilde, the equations can be solved to obtain $\tilde{\mathcal {E}}_{g}(x,s)$ and subsequently inverse Laplace-transformed to $\mathcal {E}_{g}$. The inverse Laplace transform of $\tilde{\mathcal {E}}_{g}$ is not known analytically, but it may be performed numerically. The material temperature is then proportional to a temporal convolution of the opacity weighted sum of the radiation energy densities, $\mathcal {E}_{g}$, and e^−τ (McClarren et al. 2008a).

To leverage the work done by McClarren, we constrain the g = 1 picket with w₁/w₂ ≪ 1 or w₁ ≈ 0. The g = 2 picket opacity is constrained to w₂ = . The relation between the non-dimensional optically thick picket and the non-dimensional heat capacity does not have physical justification but is merely a device to make certain expressions Laplace invertible. With p₂w₂ ≈ 1, p₂ ≈ 1/ and p₁ ≈ 1 − 1/, the value of p₁ can be comparable to p₂ despite the disparity in opacity strength. To our knowledge, the approach of relating w₂ with is novel.

With the above constraints and Q = δ(x)δ(τ), the Fourier-Laplace transformed system of equations is

$$\begin{equation} (k^{2}+3\epsilon ^{2}s^{2})\skew3\tilde{\tilde{\mathcal {E}}}_{1}=3\epsilon s \left(1-\frac{1}{\epsilon }\right), \end{equation} \tag{ 83 }$$

$$\begin{equation} (k^{2}+3\epsilon ^{2}s(s+2))\skew3\tilde{\tilde{\mathcal {E}}}_{2}=3(s+1), \end{equation} \tag{ 84 }$$

$$\begin{equation} (s+1)\skew3\tilde{\tilde{\mathcal {M}}}\approx w_{2}\skew3\tilde{\tilde{\mathcal {E}}}_{2}, \end{equation} \tag{ 85 }$$

where double tilde indicates a Fourier-Laplace transformed quantity. The remainder of the derivation follows closely from McClarren et al. (2008a) and is not repeated here. Solving the above equations yields the following kernel planar solutions for the constrained picket-fence distribution:

$$\begin{equation} \mathcal {E}_{1}=\left(1-\frac{1}{\epsilon }\right)\frac{\sqrt{3}}{2}\delta (\tau -\epsilon \sqrt{3}|x|), \end{equation} \tag{ 86 }$$

$$\begin{eqnarray} \mathcal {E}_{2}&=&\frac{1}{\epsilon }\frac{\sqrt{3}}{2}e^{-\tau } \left [\frac{\tau I_{1}(\sqrt{\tau ^{2}-3\epsilon ^{2}x^{2}})}{\sqrt{\tau ^{2}-3\epsilon ^{2}x^{2}}}\Theta (\tau -\epsilon \sqrt{3}|x|)\right. \nonumber \\ &&\left. +I_{0}(\sqrt{\tau ^{2}-3\epsilon ^{2}x^{2}})\delta (\tau -\epsilon \sqrt{3}|x|)\right ], \end{eqnarray} \tag{ 87 }$$

and

$$\begin{equation} \mathcal {M}= \frac{\sqrt{3}}{2}e^{-\tau }I_{0}(\sqrt{\tau ^{2}-3\epsilon ^{2}x^{2}})\Theta (\tau -\epsilon \sqrt{3}|x|), \end{equation} \tag{ 88 }$$

where Θ is the Heaviside function, and I₀ (I₁) is the zero-order (first-order) modified Bessel function. Setting = 1 yields McClarren's result (McClarren et al. 2008a). The properties of the P₁ equations allow for the application of a simple planar to spherical Green's function mapping (McClarren et al. 2008a),

$$\begin{equation} G_{{\rm point}}(r,\tau)=-\frac{1}{2\pi r}\left.\frac{\partial G_{{\rm plane}}}{\partial x}\right|_{x=r},\;\;r>0, \end{equation} \tag{ 89 }$$

where G_point and G_plane are the Green functions for 1D spherically symmetric and 1D planar geometry, respectively. Since Equations (86)–(88) were derived with a kernel source, $\mathcal {E}_{g}$ or $\mathcal {M}$ may be substituted into the right-hand side of Equation (89). The unitless spherically symmetric material temperature kernel is

$$\begin{eqnarray} \mathcal {M}(r,\tau)&=&\frac{3\sqrt{3}}{4\pi }e^{-\tau }\epsilon ^{2} \frac{I_{1}(\sqrt{\tau ^{2}-3\epsilon ^{2}r^{2}})}{\sqrt{\tau ^{2}-3\epsilon ^{2}r^{2}}} \Theta (\tau -\epsilon \sqrt{3}r) \nonumber \\ &&+\frac{3}{4\pi }e^{-\tau }I_{0}(\sqrt{\tau ^{2}-3\epsilon ^{2}r^{2}})\epsilon \frac{\delta (\tau -\epsilon \sqrt{3}r)}{r}.\quad\quad \end{eqnarray} \tag{ 90 }$$

The picket-fence opacity is simple to implement in IMC–DDMC. Giving the nature of our picket-fence constraint, we do not need group lumping or non-uniform particle transition probabilities. In this case, $\sigma _{P,j,n}=\bar{\sigma }$, $\gamma _{j,g,n}=p_{g}\sigma _{g}/\bar{\sigma }=p_{g}w_{g}$, and σ_{a, j, g, n} = σ_{P, j, g, n} = σ_{R, j, g, n}. We experiment with = 2, 15, or 50 cells over a spherical domain of 100 mean free paths, 400,000 source particles per time step over 100 time steps from τ = 300 to τ = 600. An MC particle propagates with DDMC when there are greater than three mean free paths per the particle's (j, g) coordinate. In order to simulate a problem with an instantaneous point source, we initialize the material temperature using the analytic solution and instantiate the initial MC particle field accordingly. Figure 2 contains 15 cell IMC–DDMC material temperature data plotted against Equation (90). The 50 cell case should be more converged at each time relative to the results in Figure 2. Indeed, in Figure 3 we see that IMC is now applied everywhere in the spectrum and the solution is more converged.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We next measure the error of IMC and DDMC relative to the P₁ temperature solution for different numbers of spatial cells. A grid convergence study for a stochastic method must have enough particles per cell to make the statistical error negligible. Another consideration involves a peculiarity of IMC specifically. Much work in IMC pertaining to the effect of the spatial grid on the continuous transport has been to formally characterize a pathology referred to as teleportation error (McKinley et al. 2003; Densmore 2011). Teleportation error occurs in IMC when there are many absorption mean free paths per cell but not many absorption mean free times per time step. Particle energies absorbed at one location in a cell may be reemitted in a subsequent time step many mean free paths away from the absorption locations (McKinley et al. 2003). Densmore (2011) demonstrates formally that a piecewise constant representation of scalar flux along with a time step that resolves a mean free time does not asymptotically converge to a correct discretization of the diffusion equation.

Despite these complications, there does exist literature to indicate that the IMC temperature error scales with Δr where Δr is a typical cell size of the simulation (if not the actual cell length for a 1D simulation). In investigating a method to emit particles at sub-cell deposition locations, Irving et al. (2011) measure total relative error of standard IMC for a 1D Marshak wave problem in planar geometry to a converged solution at several temporal and spatial resolutions. While not explicitly stated, their findings appear to yield an total relative error of about 9/J, particularly between J = 10 and J = 30 cells, for several time step sizes. Cheatham (2010) plots relative errors of IMC with respect to a gray form of the Su & Olson (1999) solutions in two cells of a simulation. This IMC error also appears to roughly scale linearly with cell width for each cell. Having set each simulation to have well over 13,000 particles generated per cell per time step and testing in a regime where IMC teleportation should be minimal, Figure 4 indicates our code indeed achieves an approximately linear scaling in L₂ relative error for both IMC and DDMC as well.

Zoom In Zoom Out Reset image size

Our results indicate that DDMC computes temperature with higher fidelity at low cell resolution for this particular problem. Having obtained good agreement with a static grid analytic solution, we present and test a manufactured solution the next section that allows for outflow and a group structure that includes scattering and absorption.

The MMS (Oberkampf & Roy 2010, p. 219) provides an avenue of code verification that is useful for problems with governing equations that are not amenable to a direct solution with a prescribed source. As the name MMS suggests, one postulates, or manufactures, a solution (McClarren et al. 2008b; Warsa & Densmore 2010). The next step simply involves incorporating this postulated answer into the system of equations to see what additional terms are produced through calculus and algebra (McClarren et al. 2008b). The additional terms can be included in an appropriate code as artificial sources. If the routines in the code function and interface correctly, the numerical experiment should reproduce the manufactured solution. Manufactured solutions have been developed by Warsa & Densmore (2010) for discrete ordinate codes modeling diffusive problems and by McClarren et al. (2008b) for planar, gray radiation-hydrodynamics problems in the optically thick and optically thin limits. We draw from these solutions as well as analytic forms presented by (Mihalas & Mihalas 1984, p. 474), and by Pinto & Eastman (2000) for high-velocity outflow problems.

Our manufactured solution has two groups and the temperature and radiation fields are constant in time and over the space of the expansion. As will be shown, in order to achieve an ostensibly simple solution, the code must have a time dependent source that can be distinct in form for each group. The source must supply energy to counteract the non-trivial adiabatic cooling of the field. The higher energy group has pure elastic scattering and the remaining group has both elastic scattering and absorption. These groups must couple through Doppler shifting in a way that appropriately balances the supply of energy to each group. Moreover, the solutions must exhibit an invariance with respect to several wavelength or frequency grid values.

The constrained material properties are

$$\begin{eqnarray} & \boldsymbol {U}(r,t)=\frac{U_{\max }\boldsymbol {r}}{R(t)}, \end{eqnarray} \tag{ 91a }$$

$$\begin{eqnarray} & \rho (r,t)=\rho (t)=\frac{3M}{4\pi R(t)^{3}}, \end{eqnarray} \tag{ 91b }$$

where $\boldsymbol {U}(r,t)$ is fluid velocity, ρ(r, t) is density, U_max is the maximum outflow speed, R(t) = R(0) + U_maxt is the outer expansion radius, r is the Eulerian position vector, and M is the total outflow mass. The frequency integrated manufactured radiation intensity and temperature are

$$\begin{eqnarray} & I_{0}(r,t) = \frac{\phi _{m}}{4\pi }, \end{eqnarray} \tag{ 92a }$$

$$\begin{eqnarray} & T(r,t) = T_{m}=(\phi _{m}/(ca))^{1/4}, \end{eqnarray} \tag{ 92b }$$

respectively. The subscript m denotes a manufactured value. For a general frequency grid, ν_{G + 1/2} < ν_{G − 1/2} < ⋅⋅⋅ < ν_1/2, and group index, g ∈ {1...G}, we constrain the absorption opacity to be pre-grouped as

$$\begin{equation} \sigma _{0,\nu _{0},a} = \sigma _{g}(t) = \kappa _{g}\rho (t), \end{equation} \tag{ 93 }$$

where κ_g is constant within group g. The differential scattering opacity has the form

$$\begin{equation} \sigma _{0,s}(\nu _{0}^{\prime }\rightarrow \nu _{0},\hat{\Omega }_{0}^{\prime }\cdot \hat{\Omega }_{0}) =\frac{\sigma _{s}(t)}{4\pi }\delta (\nu _{0}^{\prime }-\nu _{0})= \rho (t)\frac{\kappa _{s}}{4\pi }\delta (\nu _{0}^{\prime }-\nu _{0}), \end{equation} \tag{ 94 }$$

where κ_s is a constant. We now may express the manufactured multifrequency solution as an opacity-dependent superposition of the normalized Planck function and a piecewise constant function,

$$\begin{eqnarray} \varphi _{m,\nu _{0}} &=& \left(\frac{\sigma _{s}}{\sigma _{g}+\sigma _{s}}\right) \frac{\varphi _{m,s,g}}{\Delta \nu _{g}}+ \left(\frac{\sigma _{g}}{\sigma _{g}+\sigma _{s}}\right)b_{\nu _{0}},\nonumber \\ &&\times\nu _{0}\in [\nu _{g+1/2},\nu _{g-1/2}], \end{eqnarray} \tag{ 95 }$$

where $I_{0,\nu _{0}}=I_{0}\varphi _{m,\nu _{0}}$, φ_{m, s, g} is the uniform non-thermal contribution to g and $b_{\nu _{0}}$ is the normalized Planck function. With particular choices of opacity, frequency grid, and φ_{m, s, g}, the integral of $\varphi _{\nu _{0}}$ may be constrained to 1.

Incorporating the form of the opacities into the transport and temperature equations, the system to be solved with manufactured sources is

$$\begin{eqnarray} &&\left(1+\hat{\Omega }_{0}\cdot \frac{\boldsymbol {U}}{c}\right) \frac{1}{c}\frac{D I_{0,\nu _{0}}}{Dt}+ \hat{\Omega }_{0}\cdot \nabla I_{0,\nu _{0}}- \frac{1}{c}\hat{\Omega }_{0}\cdot \nabla \boldsymbol {U}\nonumber \\ &&\quad\quad\cdot \hat{\Omega }_{0}\nu _{0} \frac{\partial I_{0,\nu _{0}}}{\partial \nu _{0}}-\frac{1}{c}\hat{\Omega }_{0}\cdot \nabla \boldsymbol {U}\cdot (\mathbf {I}-\hat{\Omega }_{0}\hat{\Omega }_{0})\cdot \nabla _{\hat{\Omega }_{0}}I_{0,\nu _{0}}+\frac{3}{c}\hat{\Omega }_{0}\nonumber \\ &&\quad\quad\cdot \nabla \boldsymbol {U}\cdot \hat{\Omega }_{0}I_{0,\nu _{0}}= \frac{1}{4\pi }\sigma _{g}b_{\nu _{0}}(T)acT^{4}-(\sigma _{g}+\sigma _{s})I_{0,\nu _{0}}\nonumber \\ &&\quad\quad+\frac{\sigma _{s}}{4\pi }\int _{4\pi }I_{0,\nu _{0}}(\boldsymbol {r},\hat{\Omega }_{0}^{\prime },t)d\Omega _{0}^{\prime } + \frac{S_{m,\phi,\nu _{0}}(r,t)}{4\pi }, \end{eqnarray} \tag{ 96 }$$

and

$$\begin{eqnarray} C_{v}\frac{DT}{Dt} &=&\sum _{g=1}^{G}\int _{4\pi }\int _{\nu _{g+1/2}}^{\nu _{g-1/2}}\sigma _{g} (I_{0,\nu _{0}}-ac T^{4}b_{0,\nu _{0}})d\nu _{0}d\Omega _{0}\nonumber \\ &&+S_{m,T}(r,t), \end{eqnarray} \tag{ 97 }$$

where $S_{m,\phi,\nu _{0}}$ and S_{m, T} are the manufactured radiation and material sources to be determined. At implementation, S_{m, T} can be treated with the usual Fleck factor re-balance of material source terms. With the manufactured solutions specified, the source terms are found to be

$$\begin{eqnarray} \frac{S_{m,\phi,\nu _{0}}}{\phi _{m}} &=& \frac{U_{\max }}{cR(t)}\left(3\varphi _{m,\nu _{0}}- \nu _{0}\frac{\partial \varphi _{m,\nu _{0}}}{\partial \nu _{0}}\right)\nonumber\\ &&+\sigma _{g}(\varphi _{m,\nu _{0}}-b_{\nu _{0}}(T_{m})), \end{eqnarray} \tag{ 98 }$$

and

$$\begin{equation} \frac{S_{m,T}}{\phi _{m}}=\sum _{g=1}^{G}\sigma _{g}(b_{g}(T_{m})-\varphi _{m,g}), \end{equation} \tag{ 99 }$$

where $\varphi _{m,g}=\int _{\nu _{g+1/2}}^{\nu _{g-1/2}}\varphi _{m,\nu _{0}}d\nu _{0}$. Integration of Equation (98) over group g yields

$$\begin{eqnarray} \frac{S_{m,\phi,g}}{\phi _{m}}&=&\frac{U_{\max }}{cR(t)}\left(4\varphi _{m,g}-\nu _{g-1/2}\varphi _{m,\nu _{g-1/2}}+ \nu _{g+1/2}\varphi _{m,\nu _{g+1/2}}\right)\nonumber \\ &&+\sigma _{g}(\varphi _{m,g}-b_{g}(T_{m})). \end{eqnarray} \tag{ 100 }$$

We now construct a two group instantiation of Equations (99) and (100) that is simple to implement but is a good test of energy balance and group coupling in the code. A scattering opacity σ_s = 0.1ρ(t) is applied along with absorption opacities σ₁ = 0 and σ₂ = 0.1ρ(t). The sub-group profile of the radiation energy density in g = 1 is constant. Thus, particle frequency may be sampled uniformly in g = 1. If the group domains are implemented in wavelength, then the MC source particle wavelengths in g = 1 must be calculated as the reciprocal of a uniform sampling between the reciprocals of the group wavelength bounds. Since radiation only redshifts, the upwind group approximation (Mihalas & Mihalas 1984, p. 475) for the group-interface terms in Equation (100) must be a reasonable approach for Equation (95) and the MC process described. Equation (100) consequently may be expressed as

$$\begin{equation} \frac{S_{m,\phi,1}}{\phi _{m}}= \frac{U_{\max }}{cR(t)}\left(4\varphi _{m,1}+\nu _{3/2}\varphi _{m,\nu _{3/2}}^{+}\right), \end{equation} \tag{ 101 }$$

and

$$\begin{eqnarray} \frac{S_{m,\phi,2}}{\phi _{m}} &=& \frac{U_{\max }}{cR(t)}\left(4\varphi _{m,2}-\nu _{3/2}\varphi _{m,\nu _{3/2}}^{+}\right)\nonumber \\ &&+\sigma _{2}(\varphi _{m,2}-b_{2}(T_{m})), \end{eqnarray} \tag{ 102 }$$

where the plus superscript denotes evaluation on the right side of the frequency bound. We exploit our ability to choose a frequency grid that further simplifies the form of the source terms. Using Equation (95), if ν_5/2 = 0, φ_{m, s, 1} = φ_{m, s, 2} = 1/2 and the integral of $b_{\nu _{0}}$ over ν₀ ∈ [0, ν_3/2] is 1/2, then φ₂ = b₂ = 1/2 and the source terms become

$$\begin{equation} \frac{S_{m,\phi,1}}{\phi _{m}}= \frac{U_{\max }}{cR(t)}\left(2+\frac{1}{2}\frac{\nu _{3/2}}{\Delta \nu _{1}}\right), \end{equation} \tag{ 103 }$$

$$\begin{equation} \frac{S_{m,\phi,2}}{\phi _{m}}= \frac{U_{\max }}{cR(t)}\left(2-\frac{1}{2}\frac{\nu _{3/2}}{\Delta \nu _{1}}\right), \end{equation} \tag{ 104 }$$

and

$$\begin{equation} S_{m,T} = 0. \end{equation} \tag{ 105 }$$

Newton iteration yields hν_3/2/kT_m ≈ 3.503 to obtain b₂(T_m) = 1/2. To satisfy Equation (105), the MC process must deposit the correct energy in g = 2 directly from radiation generated in g = 2 and indirectly from redshifting radiation originating in g = 1. The strength of the group coupling is quantified with ν_{3, 2}/Δν₁.

The numerical results for our manufactured solution include strong and weak Doppler coupling for pure IMC, IMC–DDMC, and pure DDMC. For the IMC–DDMC hybrid, IMC is employed in the pure scattering group and DDMC is employed in the thermally coupled group. We note there are several ways to implement the Doppler shifting in the pure DDMC test that are equivalent. For instance, the group bound terms on the right-hand side of Equation (30) can be implemented with the upwind approximation as probabilities of redshift during the DDMC process if the scattering is elastic. For the following results, each DDMC particle has its wavelength or frequency sampled according to the sub-group distribution (which is constant in this case) after transport. The frequency value is redshifted by the same formula that lowers the particle's energy weight. If the new value of frequency is located in a new group, the particle is transferred to that group for the next time step.

For the first test, we set ν_3/2/Δν₁ ≈ 0.036 which of course is small relative to 4. The other domain quantities are set in a manner that makes adiabatic cooling of the trapped radiation field non-negligible: U_max = 10⁹ cm s⁻¹, R(0) = 1.728 × 10¹⁴ cm, t ∈ [172, 800, 181, 440] s (or 2–2.1 days), M = 10³³ g, C_v = 2 × 10⁷ρ, T_m = 1.1602 × 10⁷ K, and a wavelength grid of {λ_1/2, λ_3/2, λ_5/2} = {1.239 × 10⁻⁹, 3.542 × 10⁻⁸, 1.2398 × 10⁻³} cm.

The pertinent computational quantities are: 10 time steps, J = 10 spatial cells, 400,000 source particles generated per time step, and 400,000 initial particles. For the specification provided, we expect the code to produce reasonable agreement to the manufactured profiles. Figures 5(a) and (b) have pure IMC data for radiation energy density and temperature at several times, respectively.

Zoom In Zoom Out Reset image size

Figures 5(c) and (d) have pure DDMC data for radiation energy density and temperature at several times, respectively. The DDMC results indicate the method does not predict as much leakage at the outer cell as IMC through the observed times. DDMC appears to produce less noise near the origin for the problem depicted.

We do not show the IMC–DDMC results for this problem here but note that they too reproduce the manufactured profiles over the computed time scale. With the computation time of pure DDMC scaled to 1, Table 1 has the relative computation times of IMC–DDMC and IMC for the weak coupling manufactured source test.

We now change the coupling term to be comparable to four with a value of ν_3/2/Δν₁ ≈ 3.0. This is a more strenuous test on the code's ability to tally the correct redshift rates because the manufactured source in g = 2 is much reduced. The modified wavelength grid has λ_1/2 = 2.656838745 × 10⁻⁸ cm to implement this coupling strength. Otherwise, all quantities are the same from the weak coupling test. The IMC–DDMC radiation energy density and temperature results for this test are shown in Figures 5(e) and (f). Results for pure IMC and pure DDMC are similar.

The DDMC scaled computation times for the strong redshift coupling test are shown in Table 2.

We ensure the pure scattering group, g = 1, indirectly sustains the temperature profile's steady state by removing the radiation in that group. The profile should steadily drop relative to the solution that includes Doppler shifting. Evidence for this effect can be found in Figures 6(a) and (b).

Zoom In Zoom Out Reset image size

Finally, we construct a multigroup outflow problem that may be used to test both code efficiency and quality of the operator split treatment of grid motion in IMC–DDMC. The problem consists of a Heaviside spherical source out to 0.8U_max of strength 4 × 10²⁴/(t_n + t_min)³ erg cm⁻³ s⁻¹ in a purely absorbing fluid. The opacity is over 10 groups logarithmically spaced in wavelength from 1.238 × 10⁻⁹ cm to 1.238 × 10⁻³ cm. As in the manufactured solution, we alternate small and large opacities across the wavelength grouping where odd groups have the larger opacity. Again, the speed U_max = 10⁹ cm s⁻¹ and the total mass M = 10³³ g. Instead of equal density in each cell, we attribute M/J mass to each cell for J fluid cells. Uniformly partitioning the mass creates a density gradient that couples into the transport process through the macroscopic opacity. All simulations in this section use 50 velocity cells from 0 cm s⁻¹ to U_max, 192 time steps from 2 days to 11 days, and 100,000 source particles per time step. The heat capacity is adopted from Pinto & Eastman (2000); C_v ≈ 2 × 10⁷ρ erg cm⁻³ K⁻¹. For the calculations discussed, we use DDMC for cell j when ΔU_jΔt_nσ_g ⩾ τ_DDMC and IMC otherwise. Following Abdikamalov et al. (2012), we have labeled the threshold mean free path number between IMC and DDMC as τ_DDMC. For the first two tests discussed, depicted in Figures 7(a) and 8, τ_DDMC = 3. For the tests of the new boundary condition, depicted in Figures 7(b) and 9, τ_DDMC = 3 for Figure 7(b) and τ_DDMC = 10 for Figure 9. As a consequence of the problem's structure, the density gradient creates a "method front" for IMC–DDMC where an outer shell of IMC moves inward over the grid. The method front is heterogeneous in group space, meaning the even groups are converted to IMC sooner (or are already IMC) over the specified problem duration.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The discrepancies between the IMC and IMC–DDMC profiles in Figures 7(a) and 8 arise in part from implementing Equation (33) instead of Equation (61). For mean free path thresholds on the order of two to three per cell, we find that these errors are systematic yet generally minimal. The discrepancy can be made much worse by implementing a more conservative mean free path threshold of about 10 for this problem set. Higher mean free path thresholds require IMC particles emitted from a DDMC spatial surface to propagate through an optically thicker sub-cell environment. Of course this discrete interface is not present for pure IMC; in other words there is not a significant source of IMC radiation originating from one surface in pure IMC since the method interface is not present. Complementarily, the DDMC field is not as sourced by IMC radiation for a high mean free path threshold. This is discernible from the form of $P_{b(j,j^{\prime })}(\mu)$. We demonstrate that incorporating an approximation of the new factor, G_U, from Equation (61) indeed appears to mitigate the over-redshift near the method front in Figures 7(b) and 9. However, we recommend a mean free path threshold in the range of two to five mean free paths. Arguments presented by Pinto & Eastman (2000) indicate that diffusion theory remains valid on large outflow time scales if radiation momentum does not greatly affect fluid momentum.

In the first case tested, we use

$$\begin{equation} \sigma _{g}=\left\lbrace \begin{array}{@{}l}0.13\rho,\;\;g=2k-1\\\ms 0.13\times 10^{-4}\rho,\;\;g=2k. \end{array}\right. \end{equation} \tag{ 106 }$$

Results for Equation (106) are plotted in Figure 7(a). We find IMC–DDMC is faster than IMC by a factor of 3.36.

In the second case tested, we use

$$\begin{equation} \sigma _{g}=\left\lbrace \begin{array}{@{}l}0.13\rho,\;\;g=2k-1\\\ms 0.13\times 10^{-7}\rho,\;\;g=2k. \end{array}\right. \end{equation} \tag{ 107 }$$

Results for Equation (107) are plotted in Figure 8. We find IMC–DDMC is faster than IMC by a factor of 4.59. We have improved the performance of IMC–DDMC relative to IMC by increasing the disparity in adjacent group opacities. If instead the g = 2k − 1 opacities are increased, IMC–DDMC provides even further improvement in the diffusive groups of the spectrum. In both simulations, IMC–DDMC transitions to pure IMC over the 9 day period.

We now describe a possible implementation of Equation (61). There are two readily discernible means of ascribing a MC interpretation to the G_U(μ) factor. In one, the probability that an IMC particle transmits into DDMC when incident on a diffusive region may be taken as P(μ) ∼ (1 + 3μ/2)G_U(μ). However, this would make P(μ) > 1 for some μ regardless of cell material properties. To constrain P(μ) ⩽ 1, some range of angular projections, μ ∈ [0, ε], where ε < 1, would have to have G_U modified. As an alternative, the probability of IMC-to-DDMC transmission may be maintained as its original form P(μ) ∼ 1 + 3μ/2 and the particle weight must then be multiplied by G_U(μ). Since G_U(μ) is of O(1/μ), the energy current across the IMC–DDMC interfaces is bounded. Physically, it is supposed that this implies the importance of IMC particles interacting with a surface at a DDMC region to the energy balance in the adjacent cells must still be bounded at small values of μ. In our tests, for an inner DDMC region adjacent to an outer IMC region at a boundary with speed U,

$$\begin{equation} G_{U}(\mu)=1+2\frac{U}{c}\left(\frac{C_{1}}{\mu }-C_{2}\mu \right), \end{equation} \tag{ 108 }$$

where C₁, C₂ > 0 are constants. In passing, we note that a angularly uniform, heuristic G_U may be calibrated from simulation. From phenomenological considerations, it is found that a good form of G_U is

$$\begin{equation} \bar{G}_{U}=1+2\min (0.055\tau _{{\rm DDMC}},1)\frac{U}{c} \end{equation} \tag{ 109 }$$

at least for a range of τ_DDMC ∈ [3, 10]. Applying Equation (108) with C₁ = 0.55 and C₂ = 1.25 for τ_DDMC = ΔU_jΔt_nσ_g = 3 and Equation (106), we see a small improvement in Figure 7(b).

To further demonstrate the potential utility of Equation (108), we apply the G_U factor to a problem where the discrepancy is made very large by setting τ_DDMC = 10. The fitting constants for this test are C₁ = 0.6 and C₂ = 1.25. We observe a significant improvement at all times including when the discrepancy is very large at 1.5 days, as plotted in Figure 9. However, with this improvement comes some additional MC noise due to the weight modifications having a large range of G_U and insufficient sampling for μ → 0. Since the IMC portion of the simulation is in the lab frame, the minimal comoving projection into the DDMC interface cell is U/c. But physically, the comoving projection lower bound is 0. For sampling comoving directions with μ < U/c, the new boundary condition is applied to particles with μ < U/c that are advected onto DDMC surfaces through the velocity position shift algorithm delineated in Section 2.4. We find that these samplings are not important for τ_DDMC = 3 but are important for τ_DDMC = 10 or greater. At larger τ_DDMC, IMC particles that "graze" the DDMC surface are important due to the increased level of angular isotropy near the IMC–DDMC boundary. We note that the proper implementation of the weight modifying factor, G_U, is an open research question. Additionally, our choice of asymptotic scalings and analysis in Section 2.3 is approximate and tailored for homologous outflow.

We caution that what may be thought of as a conservative selection of a mean free path threshold between IMC and DDMC may lead to the types of errors depicted in Figure 9.