HIGH-TEMPERATURE PROCESSING OF SOLIDS THROUGH SOLAR NEBULAR BOW SHOCKS: 3D RADIATION HYDRODYNAMICS SIMULATIONS WITH PARTICLES

In this section, we show the results of two bow shocks, each of which was run with a fixed adiabatic index γ = 1.46 and a wind V_w = 9 km s⁻¹. One simulation explores the shock in two dimensions (N_x, N_y = 300, 200 cells), while the other explores the fully 3D case (N_x, N_y, N_z = 300, 200, 200 cells). A Cartesian grid is used in both simulations, as has been done in some previous 2D studies. The 2D case effectively simulates a bow shock around a cylinder, while the 3D case simulates a bow shock around a sphere. The results are shown in Figure 1.

Zoom In Zoom Out Reset image size

The overall morphologies of the shocks are similar. Each has a strong bow shock, with very high temperatures directly in front of the planetoid. There is also a strong tail shock with a very high temperature (low density) wake. Nonetheless, there are two fundamental differences that will have an impact on characterizing solid processing. (1) The standoff distance in the 2D shock (∼0.5R_e) is much larger than in the 3D shock (∼0.2R_e). Simulations in two dimensions overestimate the retention of solids, as the available stopping distance is much larger than in 3D simulations. Standoff distances in simulations using a proper equation of state will be discussed in the next section. (2) The opening angle of the bow shock is larger in two dimensions than in three dimensions. This difference will overestimate the available cross section for strong shocks and high temperatures. Thus, the amount of solids that can be processed by a bow shock will be overestimated. Moreover, the 2D shock morphology gives a false impression that significant processing can occur at impact parameters larger than R_e, which is not seen in 3D simulations.

Both of these effects seen in 2D shocks (larger standoff distance and larger opening angle) will overestimate the efficiency of solid alteration through bow shocks. Unless significant care is taken in accounting for the 3D geometry in a 2D code, a 2D bow shock simulation will give unreliable results. For the remaining calculations in this manuscript, we focus on fully 3D simulations.

The 3D simulation in the previous section used a fixed adiabatic index for the equation of state. While γ = 1.46 takes into account the rotational states of molecular hydrogen, it does not include the vibrational states or dissociation. These internal processes are a tremendous thermal buffer. Here, we explore the morphology of three adiabatic bow shocks using a full equation of state, with V_w = 7, 8, and 9 km s⁻¹. Solids are tracked through these simulations and thermal histories of solids that pass through the shock are recorded, assuming perfect thermal coupling between the gas and the solids.

Figure 2 shows the structure of a 7 km s⁻¹ bow shock. Density, temperature, and pressure are all shown as a slice through the shock and center of the planetoid. Particles that are within the cells of the given slice are shown in blue. For one of the temperature plots (Figure 3), trajectories for 20 particles that have impact parameters between 0 and R_e are shown as blue curves. As in the simulations with a fixed adiabatic index, strong bow and tail shocks form, as well as a high temperature, low density wake. Peak gas temperatures are just over 2000 K, and peak pressures are ∼500 μBar. The standoff distance is about 0.17R_e, enabling most solids to escape accretion onto the body. Surviving solids, i.e., solids that are not swept up by the planetoid, are diverted from the wake, creating a dust-free region. We will return to implications of this dust-free region in the Section 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We deconstruct the shock into five regions: the pre-shock, the bow shock, the post-bow shock, the tail shock, and the post-tail shock. The dust-free wake could be considered a sixth region, but we exclude it for the following discussion as solids do not enter it. Particles and gas enter the simulation at the positive x boundary. Besides a modest acceleration toward the planetoid due to the planetoid's gravity, the pre-shock region is relatively undisturbed. This acceleration will be much more important for planetary-mass objects. The temperature and pressure of the gas jumps suddenly at the bow shock. Near the centroid of the bow shock, the temperature and pressure profiles continue to rise as the solids approach the planetoid's surface. The curves that end abruptly result from solids that strike the planetoid.

After a few minutes in the post-shock region, the gas temperature and pressure drop abruptly (see Figure 3). This drop is solely due to adiabatic expansion, and is not captured in 1D models. The cooling rate after the shock is about 1000 to 5000 K hr⁻¹, depending on the impact parameter of the solid and the time after the shock. These high cooling rates (a few 1000 K hr⁻¹) last for about 10 to 20 minutes. Some curves show oscillations, which are due to the discretization of the planetoid onto the grid, causing a series of minor shocks. These oscillations, which are artifacts, are most severe for particles that barely miss the surface of the planetoid, but the general trend is indicative of the real behavior.

After the particles have passed through the bow shock and have cooled in the post-bow shock region, they encounter the tail shock, which is produced as gas is diverted around the planetoid toward the wake. For this particular case, some particles drop below 1400 K, roughly the lower limit of the crystallization temperature, for a minute or two before they are heated by the tail shock back to or above 1500 K. The cooling rate in the post-tail shock region then becomes significantly lower (∼100 K hr⁻¹) than in the post-bow shock. While the initial cooling is very fast, the post-tail shock region can potentially set very protracted periods of cooling for the melt, in the adiabatic limit.

The pressure profile demonstrates that very high pressures are attained at the bow shock (∼500 μBar), but drop quickly to tens of μBar in about 20 minutes. For comparison, the pressure in the pre-shock region is about 10 μBar. While the pressure does increase in the tail shock, it never approaches the pressures seen at the bow shock.

The 8 and 9 km s⁻¹ shocks show very similar morphologies to the 7 km s⁻¹ shock. Figures 4 and 5 show the results for the 8 km s⁻¹ shock, and the same is done for the 9 km s⁻¹ shock in Figures 6 and 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In each case, we can break the morphology into the same five regions (pre-shock, bow shock, post-bow shock, tail shock, and post-tail shock). The post-bow shock region again shows very prodigious cooling rates due only to adiabatic expansion. At the higher bow shock speeds of 8 and 9 km s⁻¹, some particles are prevented from dropping below 1400 K, and the tail shocks are stronger than in the 7 km s⁻¹ case. Again, the post-tail shock region has very protracted cooling.

These simulations reveal the basic morphology of bow shocks and the heating and cooling that solids will encounter. Shock heating and adiabatic cooling make these curves non-trivial, with a striking divergence from truly 1D shocks, such as those that might occur due to, e.g., large-scale spiral arms. We now build on the above complexity by introducing gray radiative transfer.

Having characterized the heating and cooling curves in the 3D adiabatic limit, we now include radiative transfer using the scheme described in Appendices A.4 and A.6. We assume that any radiation that strikes the surface of the planetoid is reradiated in the direction it came from. Line cooling is not included at this time. The only solids that contribute to the opacity are assumed to be chondrule precursors (see Section 2 and Appendix A.6).

We present three realizations of bow shocks with radiative transfer. Each realization is for an 8 km s⁻¹ shock, but with chondrule precursor mass opacities of κ_c = 1, 10, and 100 cm² g⁻¹ of solid, labeled Rad 0.1, Rad 1, and Rad 10, respectively. The naming convention reflects the reduction/enhancement of the opacity relative to the "standard" opacity 10 cm² g⁻¹ of solid, which is based on our assumed precursor concentration. Figures 8 and 9 show the results for Rad 1, where the plots are similar to Figures 2 and 3, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The standoff distances for Rad 1 are between 0.1 and 0.13 R_e compared with 0.17 R_e for the 8 km s⁻¹ adiabatic case. The highest temperatures are reached at the bow shock, and unlike the adiabatic cases, the temperatures drop immediately, even when approaching the planetoid. Figure 8 shows that the opening angle is smaller in Rad 1 than in the adiabatic case due to the radiative case's lower temperatures in all regions behind the bow shock. A tail shock still forms, but it is much weaker than in the adiabatic simulations. The only region behind the bow shock that remains at very high temperatures is the dust-free wake, which cannot radiate due to the lack of dust.

Figure 9 shows the trajectories for 20 solids, as done for Figure 3. The curves show that all solids experience a radiative precursor, but this precursor does not increase temperatures to more than 500 K. The radiative bow shock produces a striking spike in temperature, with a sudden rise and drop in the pre- and post-bow shock regions, respectively. The pressure also shows rapid evolution, but does not drop as quickly as the temperature due to the increase in density in the post-bow shock region. While the tail shock produces some mild heating, the post-tail shock region is not substantially hotter than the pre-shock region, having temperatures between 400 and 600 K.

At the shock front, the cooling rates approach 10,000 K hr⁻¹ and remain above ∼4000 K hr⁻¹ through the crystallization temperatures. Radiative transfer is extremely efficient at removing energy from the bow shock that would otherwise keep the chondrules at about 1400 K for protracted periods.

Next, we examine the results from Rad 10 in Figures 10 and 11. For these simulations, we assumed that the opacity is 100 cm² g⁻¹ of solid, 10 times larger than what we expect for our assumed chondrule concentration. There are two main changes compared with Rad 1. First, the radiative precursor is much more significant. This can be seen in both the temperature plots in Figure 10 and the temperature profiles for 20 solids in Figure 11. However, this radiative precursor also redistributes a significant amount of energy from the shock itself, causing the peak temperatures of the bow shock to reach a value between 1500 and 1600 K, significantly lower than seen in any other run. The temperature profiles show a radiative precursor that reaches ∼1000 K for some solids just before the shock. The pressure spike is very similar to the Rad 1 case, except there is a slight increase in pressure before the bow shock due to the radiative precursor. The cooling rates are smaller for Rad 10 than Rad 1, but are still between 2000 and 4000 K hr⁻¹ immediately behind the bow shock. However, the cooling rates are from a peak temperature below or barely in the crystallization temperature range, even for an 8 km s⁻¹ shock. The post-tail shock region is also cooler than in the Rad 1 case, as this region is able to radiate more efficiently.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While one might expect a higher opacity to prevent fast cooling, this is only true if the optical depths are already very large. The results here suggest that the mass opacity will need to be unrealistically high to reach this limit of long radiative cooling timescales. The higher optical depth in the Rad 10 simulation also promotes an efficient propagation of the bow shock energy into the pre-shock region, limiting the peak temperatures of the shock. However, just as very high optical depths can suppress cooling due to long radiation propagation timescales, very low optical depths can also suppress cooling due to a lack of surface area. The low opacity limit may be more plausible for Solar Nebula bow shocks, and we consider it next.

For the Rad 0.1 case, we reduce our standard opacity by a factor of ten, resulting in κ_c = 1 cm² g⁻¹ of solid. While it makes the gas more optically thin than in the Rad 1 simulation, it also prevents the gas from radiating efficiently, as the standard opacity is already in the low optical depth regime (see Section 4). Figure 12 shows the shock morphology, and it is closer to the adiabatic behavior than either the Rad 1 or Rad 10 simulation. The opening angle and standoff distances are similar to those for the 8 km s⁻¹ adiabatic case. The tail shock once again becomes very apparent. Nevertheless, Figure 13 shows that radiative cooling causes a consequential deviation from the adiabatic results. Cooling through the crystallization temperature range is about 6000 K hr⁻¹. The tail shock does cause strong heating, but the temperatures do not increase above 1500 K. Although lower opacity simulations are not presented here, it is straightforward to infer that a lower opacity will recover the tail shock heating and post-tail shock cooling seen in the adiabatic limit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the next section, we discuss some of the implications of these results.

Boxzy Hydro is a Cartesian, second-order accurate hydrodynamics code that has the following capabilities.

• 1.  
  Hydrodynamics fluxing using either linear or angular momentum.
• 2.  
  A detailed equation of state that includes the rotational and vibrational states of H₂ as well as H₂ dissociation.
• 3.  
  Self-gravity.
• 4.  
  A radiation algorithm that is valid without directional bias and couples the thick and thin regimes.
• 5.  
  Particle integration using particle-in-cell and particle–particle integration. Drag terms with feedback on the gas are included.
• 6.  
  Inclusion of arbitrary structures in the grid.

Each of these items will be discussed in turn. First, the hydrodynamics equations are solved in the following form

$$\begin{eqnarray} \frac{\partial \rho }{\partial t} + \nabla \cdot \left({\bf v}\rho \right) &=& 0 \quad \mbox{(mass continuity)} \nonumber \\ \frac{\partial {\bf v}\rho }{\partial t} + \nabla \cdot \left({\bf v} {\bf v} \rho \right) &=& {-}\nabla P - \rho \nabla \Phi\quad \mbox{(momentum equation)}\nonumber \\ \frac{\partial {E}}{\partial t} + \nabla \cdot \left({\bf v} \left(E+P\right) \right) &=& {-}\rho {\bf v}\cdot \nabla \Phi ,\quad \mbox{(energy equation)},\nonumber \\ \end{eqnarray} \tag{ A1 }$$

where E = (1/2)ρ|v|² + for gas density ρ, velocity v, and internal energy density . In the absence of a gravitational potential field Φ, the energy equation is written in conservative form, allowing total energy conservation to machine precision. However, this degree of precision is rarely reached in astrophysical studies due to gravitational fields (see, however, recent developments by Jiang et al. 2013). The pressure P is described in more detail below. For the moment, focus on the momentum equation. Note that while the equation is conservative for the linear momentum, it makes no such guarantees for the angular momentum. This can have dire consequences for simulations of rotating objects on Cartesian grids. For example, rotating polytropes in a code fluxing linear momentum can either gain or lose significant angular momentum. Any user can verify this behavior in their code of choice by running an n = 1 polytrope with modest initial rotation. In the Boxzy Hydro, angular momentum tends to decay if the equations are solved as given above. To address this, the following two methods are available for momentum fluxing. (1) The linear momentum, as shown above, can be solved in each direction, as is standard in Cartesian grid codes. (2) As an alternative, the z direction can be solved using linear momentum, but the x and y momenta can be transformed into radial and angular momenta and fluxed accordingly. For this second option, although the grid is Cartesian, angular momentum in the z direction can be conserved in the absence of potential field errors. As a result of the transformation for fluxing, the force calculation (sourcing) requires additional terms such that the right-hand side of the momentum equation becomes

$$\begin{eqnarray} \rho v_x & += & (-\partial _x P -\rho \partial _x \Phi + \rho (x v_y - y v_x)^2 x/r^4) \Delta t, \end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray} \rho v_y & += & ({-}\partial _y P -\rho \partial _y \Phi + \rho (x v_y - y v_x)^2 y/r^4) \Delta t,\ {\rm and}\nonumber \\ \end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray} \rho v_z & += & \left(-\partial _z P -\rho \partial _z \Phi \right)\Delta t. \end{eqnarray} \tag{ A4 }$$

The major trade off is that a singularity is now introduced in the grid, so the choice of fluxing is problem-dependent and must be used with caution. At this time, the viscous stress tensor is not included, but is in development.

Fluxing is performed using an approximate Riemann solver (HLLE) with a generalized MINMOD flux limiter (e.g., Toro 2009). In practice, the calculation of the pressure force is incorporated into the fluxing routine. Consider flux variables $\bf {F}$ and conserved variables $\bf {U}$ in, for example, the x direction, where

$$\begin{equation} U_x =\left(\begin{array}{@{}c@{}}\rho \\ \rho v_x\\ \rho v_y\\ \rho v_z\\ E\\ \end{array}\right), \end{equation} \tag{ A5 }$$

and

$$\begin{equation} F_x =\left(\begin{array}{@{}c@{}}\rho v_x\\ \rho v_x v_x +P \\ \rho v_y v_x\\ \rho v_z v_x\\ (E+P) v_x\\ \end{array}\right). \end{equation} \tag{ A6 }$$

All variables are cell-centered on the Cartesian grid. The hydrodynamics is solved for each direction (x, y, and z) independently, but the conserved variables are only updated after the contribution from each direction has been calculated. The hydrodynamics time step Δt_CFL is limited according to the Courant–Friedrichs–Lewy condition Δt_CFL[α_x + α_y + α_z] < 1, where α_x = (c_s + |v_x|)/Δx for adiabatic sound speed c_s. We typically set the time step to be Δt_CFL = 1/(2[α_x + α_y + α_z]).

The pure hydrodynamics equations (no gravity or radiation) are evolved by advancing the initial state Uⁿ over a full time step, which gives a predicted state Uⁿ*. The solution is then corrected using a second flux calculation based on the predicted state. The scheme is second order in time and space. The following example shows how this method is applied in the 1D limit:

$$\begin{eqnarray} {U}^{n*}&=&{ U}^n+\lambda [F(U^{n})_{+1/2}-F(U^{n})_{-1/2}] \end{eqnarray} \tag{ A7 }$$

$$\begin{eqnarray} {U}^{n+1} &=&(U^{n}+{ U}^{n*} + \lambda [ F(U^{n*})_{+1/2}-F(U^{n*})_{-1/2}])/2, \nonumber \\ \end{eqnarray} \tag{ A8 }$$

where λ = Δt_CFL/Δx for CFL-limited time step Δt_CFL. The quantities F(U)_+1/2 − F(U)_−1/2 represent the flux differences for any given cell based on the fluxes at the + and − cell faces, respectively. These flux values are determined using the approximate HLLE Riemann solver, in combination with a second-order upwind interpolation scheme for the hydrodynamics variables. Slopes for the upwind interpolations are limited (flux-limiter) by the generalized MINMOD algorithm. Unless otherwise stated, the generalized MINMOD limiter is set to the standard MINMOD limiter (diffusion parameter θ = 1). Please see Toro (2009) for additional details. Gravity and radiation transfer are included by sourcing Uⁿ for one-half of a time step before advancing the hydrodynamics (before Equations (A7) and (A8)). After the hydrodynamics step is completed, the gravitational potential is updated using the U^{n + 1} solution for the density. Using the U^{n + 1} state, a final gravity and radiation update over one-half of a time step is performed. This strategy for including gravity and radiation is analogous to a kick-drift-kick method.

The numerical method is tested using Sod (1978) shock tubes. For this test, a perfect gas with equation of state p = (γ − 1) is placed in a control volume. The gas is given left and right states, with one state at a higher pressure than the other. When the system is allowed to evolve, a shock, contact discontinuity, and rarefaction wave develop, which have analytic solutions. We show a typical example in Figure 15, where we use ρ_l = 1 and p_l = 1 for the left states and ρ_r = 0.125 and p_r = 0.1 for the right states, each in code units. The test uses a 1D grid with 200 grid cells. For this test, we explore diffusion parameter values of θ = 1 and 2. While θ = 2 leads to sharper transitions (less diffusive), it also leads to a slight dip in density without a compensating dip in pressure, which is a spike in entropy. For this reason, we typically set θ = 1. No other form of artificial diffusion is used, and no artificial viscosity is employed.

Zoom In Zoom Out Reset image size

The equation of state is an extension of that used in Boley et al. (2007) and Galvagni et al. (2012). A mixture of metals, helium, and hydrogen is used to compute tables for the relations between the gas's internal energy density, adiabatic index, and temperature. With these relations, an ideal gas law can be used that takes into account the partitioning of the gas between translational and internal degrees of freedom, i.e., p = ρR_gT(, ρ)/μ(, ρ) for gas constant R_g and mean molecular weight per proton mass μ. Because ionization is not yet included, only translational degrees of freedom are used for computing helium and metal partial pressures. The contribution from H₂ is much more complex. First, the equation of state is derived directly from the partition function, so the effects of heating and cooling due to partitions between translational and internal molecular degrees of freedom are directly included. The internal degrees of freedom include the rotational and vibrational states of molecular hydrogen for different para- and orthohydrogen ratios (see Boley et al. 2007) and molecular hydrogen dissociation. This behavior of H₂ for a range of thermodynamic conditions is shown in Figures 16 and 17. At very low temperatures, rotational states are inaccessible, and a pure molecule hydrogen gas will have an adiabatic index γ = 5/3. At temperatures T ∼ 300 K, rotation adds two degrees of freedom, and the gas behaves as γ = 7/5. At even higher temperatures, vibrational modes are activated, which also add two degrees of freedom (partitioned between kinetic and potential "spring" energies), giving γ → 9/7. When T ≳ 2000 K, depending on gas density, H₂ begins dissociating. This creates an effective heat sink, relative to a fixed-γ gas, where internal energy goes into breaking apart molecule H₂ instead of providing pressure support for the gas, ultimately driving γ → 1. After dissociation is complete, the gas again behaves like a γ = 5/3 gas. These tables allow accurate modeling of the compressibility of the gas and naturally take into account the effective heating and cooling rates due to internal molecular processes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Self-gravity is included by combining multi-grid techniques with a tree algorithm. First, the mass on the finest grid is averaged down onto subsequent levels. For each node of the tree, the center of mass and total mass of its eight children/leaf cells are included. As in a standard tree algorithm (Barnes & Hut 1986), we use an opening angle to determine whether a node's contribution to a given local potential is sufficient, or whether the cell needs to be opened, i.e., use the finer levels as nodes. For the algorithm here, we define the opening angle to be θ = (dx² + dy² + dz²)^1/2/r, where r is the distance between the point of interest and the center of the parent cell that has sizes dx, dy, and dz. The opening angle θ is adjustable according to the specific problem, but θ = 0.7 is typically used in the simulations presented here. For any given node's contribution to the potential, we use ϕ_node = −M_node/|r − r_com|, where M_node is the node's mass, r is the distance to the point where the potential is needed, and r_com is the distance to the given node's center of mass (COM). Using the COM of the node automatically includes the dipole moment of the cell in the potential calculation. At this time, the quadrupole and higher moments are not included, as the tree is only needed to calculate an approximate solution.

The tree potential calculation is used to find the grid boundary values, including the value at the finest level, as well as the potential solution everywhere for the penultimate grid level (or coarser if desired). The approximate solution on the penultimate finest level is then successively relaxed until a minimum residual for the grid, set by the user, is obtained. Here, the residual $\mathcal {R}=(\nabla \phi \cdot {\bf dS}-4\pi \rho dx dy dz)/L$, where L = 2(dydz/dx + dxdz/dy + dxdy/dz). For each iteration i in the relaxation, $\phi _{i+1}=\phi _i+\alpha \mathcal {R}$, where we found α < 1 provides stability. The relaxed solution on the penultimate finest level is then prolongated to the finest level, except the boundary cells. Relaxation is then performed on the finest grid as done on the previous level. For the simulations shown here, the maximum $\mathcal {R}$ on the grid is typically restricted to be less than 10⁻⁴. This method is reasonably fast and accurate.

To demonstrate the accuracy of the potential solution, we present two tests. First we place an n = 1 gas polytrope on the grid, which has the analytic density solution ρ = ρ₀sin (πx)/(πx) for x = R/R_surface. The polytrope's radius is set to 1 AU, the maximum density to 2 × 10⁻³, and the grid resolution Δx = Δy = Δz = 0.025 AU. We set the opening angle to θ = 0.7 and set the maximum residual $\mathcal {R}_{\rm max}<10^{-4}$. The potential is first calculated using the tree+relaxation (TR) scheme described above and is then compared to the solution based on direct summation, where each cell is treated as a point mass. Figure 18 shows that the result has a maximum deviation of 0.1%. The actual solution will be slightly different from both the TR and direct summation solutions due to discretization effects. We did compare the TR solution against an analytic solution for a constant density sphere, and for the same parameters, the solution was better than 0.1% everywhere.

Zoom In Zoom Out Reset image size

The second test uses a hydrodynamics calculation with self-gravity. An analytic polytrope with n = 1 is placed on the grid and then allowed to evolve, with the equation of state set to P = Kρ² for polytropic constant K. For this test, we set R = 1 AU and the central density ρ_c = 1.19 × 10⁻⁹ g cm⁻³ (2 × 10⁻³ in code units), which gives the clump a mass of about 2.5 M_J. The grid is set up to have 96³ cells, with grid spacing Δx = Δy = Δz = 0.025 AU. Figure 19 shows the results after about nine dynamical times, where we take the free-fall time ($t_{\rm ff}=(3\pi /[32 G \bar{\rho }])^{1/2}$) to be representative of the dynamical time.

Zoom In Zoom Out Reset image size

Radiative transfer plays a fundamental role in heating and cooling astrophysical objects, which in turn has a substantial impact on an object's evolution. Unfortunately, integrating radiation transport into simulations can be a technical challenge. Diffusion schemes are accurate in regions of high optical depth, but ultimately solve the wrong equation near a system's photosphere, where energy is actually lost. Ray tracing schemes do well for the optically thin limit and will capture energy loss near the photosphere, but can become problematic at high optical depths where temperature gradients, which are used in the diffusion limit, are more informative of the flow of radiation. In addition, and regardless of the scheme, there is significant potential to have wildly different radiation and hydrodynamics timescales. In this case, an explicit scheme can fail severely if Courant time steps are taken. One could choose instead to use an implicit scheme, which will be numerically stable but will not guarantee accuracy.

Here, we seek to describe a new radiation transport algorithm that combines elements of diffusion, ray tracing, and Monte Carlo methods. Our goal is to develop a radiation scheme that works well in the optically thick and thin limits, as well as in the transition regime between these limits. We also want to attempt to resolve the radiation hydrodynamics timescale. With these considerations in mind, we will first discuss the diffusion limit, followed by the optically thin limit, and then describe how we have combined these cases to form what we call the Monte Carlo-Radiative-Transfer-Flux-Limited-Diffusion (MCRTFLD) method.

Whenever the optical depth through a region of gas is very large (Δτ ≳ 10), radiative transfer is well approximated by diffusion. Consider the radiative flux in the direction x. In the diffusion limit,

$$\begin{equation} F_x = \frac{16 \sigma \beta T^3}{\rho \kappa }\frac{\Delta T}{\Delta x}, \end{equation} \tag{ A9 }$$

where

$$\begin{equation} \beta = \frac{2+y}{6+3y+y^2} \end{equation} \tag{ A10 }$$

for

$$\begin{equation} y=\frac{4}{\rho \kappa T} \left \vert \frac{\Delta T}{\Delta x}\right \vert. \end{equation} \tag{ A11 }$$

The β term is the radiative flux limiter, which permits extension of the radiative diffusion treatment into the optically thin regime by ensuring that the divergence of the flux never exceeds the value corresponding to free-streaming. In this way, the flux limiter is an interpolation between the optically thick and thin regimes, where the above β is from Bodenheimer et al. (1990). However, this reliance on an interpolation scheme is a potential problem as well. Radiative diffusion by itself does not remove energy from the system; it just diffuses it. Radiative losses, which take place near the photosphere, are determined by the flux limiter, so caution needs to be taken when applying a flux-limited diffusion scheme. Moreover, because radiation is allowed to diffuse, some temperature structures can get washed away. The severity of this effect is dependent on the desired accuracy and the problem at hand.

Now consider the optically thin regime (Δτ ≪ 1). In this case, a gas can emit as much radiation over 4π as its emissivity allows. For a source function S = σT⁴/π, the luminosity per volume of gas is ℓ = 4ρκσT⁴. If the contributions of the incident irradiation, j, on this same volume are taken into account, then ∇ · F = j − ℓ. For this relation, we have already folded in the solid angle components, so j and l are flux-like quantities summed over the area of the volume of gas.

Many astrophysical objects span both the optically thin and thick regimes. Thus, a method for coupling the two limits is highly desired. To reiterate, most of the cooling will occur in the transition region between these limits. We build this coupling by readdressing radiative transport in the thin regime. If we have a large volume of gas that represents, for example, a cell in a simulation, then assuming constant density ρ, temperature T, and opacity κ across a cell width Δx, the intensity at the cell face due to the gas in the cell alone is

$$\begin{equation} I=\frac{\sigma T^4}{\pi } \left(1-\exp {\left[-\Delta \tau \right]}\right), \end{equation} \tag{ A12 }$$

where we have assumed black body emission and Δτ = ρκΔx. The amount of energy carried away from the cell by emission can be calculated by discretizing the cell into six directions, one for each cell face. Each cell face then corresponds to a portion of the total emission of the cell over the entire sky, i.e., 2π/3. The luminosity per cell volume V that is carried away in any single direction is then δℓ = (2πIA/3V), where A is the area of a given cell face. In the limit that Δτ ≪ 1, then the luminosity per volume ℓ = 6δℓ → 4ρκσT⁴, which is the free-streaming limit.

So far, we have not discussed how the energy is transferred from one cell to other cells using the above assumptions, nor have we discussed how j is ultimately determined. Both can be calculated using particles to represent energy packets. At the face of each cell, we launch particles in a random direction, constrained to be within δθ = 60° of the normal to the cell face. Allowing random angles ensures that directional bias is removed, although it does permit the occurrence of noise. Each particle is then allowed to propagate a distance of $\Delta s=0.5 \sqrt{\Delta x^2+\Delta y^2 +\Delta z^2}$ along the launch angle. During that propagation, the energy is absorbed by δℓexp (− Δτ), where the optical depth is based on the current cell conditions such that Δτ = ρκΔs. The local value of j is likewise determined by summing the energy that is absorbed by a local cell for all particles that enter the cell. For simplicity, only the nearest node is used for determining the energy absorbed during any one single propagation.

A possible disadvantage to this method is that it can easily become time consuming. However, this is a tunable feature of the algorithm. If a poorer treatment of radiative transfer can be tolerated, then combinations of maximum propagation limits, emissivity thresholds, and computational volume limits can be used to increase performance. This portion of the algorithm also tends to be stable for large time steps, as energy is deposited over many cells unless an energy packet propagates into a very high optical depth region. This takes us back to the FLD approximation.

When the optical depth across a cell becomes very high, the variation of the source function along a cell becomes an important component for understanding energy transfer, and the above assumption of a constant temperature for a cell leads to inaccuracies in computing the cell-to-cell coupling. In the limit Δτ ≫ 1, radiative diffusion is an accurate description of energy transfer in the gas. In this spirit, we combine the FLD and the energy packet propagation algorithms by setting

$$\begin{eqnarray} &&\nabla \cdot {\bf F}=j-(\ell \exp [-3 (\Delta \tau)^2] + \ell _{\rm FLD}(1-\exp [-3 (\Delta \tau)^2])). \nonumber \\ \end{eqnarray} \tag{ A13 }$$

The FLD term is calculated using Equation (A9) using the following method: For each cell, the standard FLD flux is calculated for each cell face. If the flux is outward from the cell, the value is kept. If it is inward, the value is set to zero. The flux is converted into a luminosity per volume using F_xA/V. The energy is then propagated and absorbed in the same way as that used in the low optical depth limit. This limits exactly to the FLD formalism when the cells are at high optical depth, but allows smooth, directionally debiased propagation of energy in regions where the optical depth becomes low. The optical depth-based interpolation ensures both a smooth and sharp transition between the thick and thin regimes and the corresponding algorithm, where the form of the interpolation was ultimately chosen from trial and error with the tests given below. We finally note that the above description assumes constant volume between cells. For codes where the physical cell size changes, the algorithm needs to be reformulated using luminosity instead of luminosity per volume, which is a straightforward modification.

We build our first test case by considering a spherical distribution of gas with a luminosity point source at r = 0 of L = 10⁻⁴L_☉. The temperature gradient stellar structure equation can be expressed as

$$\begin{equation} \frac{dT}{dr} = -\frac{3 \rho \kappa L}{64 \pi \sigma r^2 T^3}. \end{equation} \tag{ A14 }$$

If ρκ is parameterized, then a solution for the temperature profile can readily be calculated. A function for ρκ that varies rapidly in r will test the algorithm's behavior in both the thick and thin regimes, and will also test the interpolation between the limits. For this reason, the first set of tests uses ρκ = ρ₀κ₀((r/r_c)⁴ + 1)⁻¹, which yields the solution

$$\begin{equation} T(r)^4=T_0^4+\frac{3 \rho \kappa L}{128 \pi \sigma r_c }\left(f\left(r\right)-f\left(r_0\right)\right), \end{equation} \tag{ A15 }$$

where

$$\begin{eqnarray} f\left(r\right)&=&\sqrt{2}\ln \left[ \frac{ r^2+\sqrt{2} r r_c + r_c^2}{ r^2-\sqrt{2} r r_c + r_c^2}\right] + 2\sqrt{2}\nonumber \\ &&\times \left(\arctan \left[1-\sqrt{2}\frac{r}{r_c}\right]-\arctan \left[1+\sqrt{2}\frac{r}{r_c}\right]\right)-8\frac{r}{r_c}. \nonumber \\ \end{eqnarray} \tag{ A16 }$$

T₀ and r₀ correspond to values at the sphere's surface.

In this test, the fluid is not permitted to move on the grid, isolating the radiation transfer algorithm. We use ρ₀κ₀ = 10⁻⁷ and 10⁻⁸ cm⁻¹ for the cored power-law profile. The sphere is truncated at r = 4 AU and r_c = 0.1 AU. The former tests how well the algorithm transitions from optically thick to thin, while the latter creates a mostly optically thick sphere. For the analytic solutions, we assume the radiative zero solution, i.e., we take T = 0 at the surface. The solution should not be expected to extend to regions of very low optical depth, as energy transport becomes very non-local. We should expect the analytic and the simulation solutions to converge rapidly near a radially integrated optical depth τ_r ∼ 1, for integrations going from r = ∞ to 0. The simulations are run until the average net energy input on the grid is zero. The results are shown in Figure 20, which show very good agreement with the analytic solution for the ρ₀κ₀ = 10⁻⁷ cm⁻¹ case. The analytic solution matches the ρ₀κ₀ = 10⁻⁸ cm⁻¹ solution well until the atmosphere drops to very low optical depths. At these depths the gas temperature becomes hotter than expected in the analytic solution because the emissivity of the gas becomes small, preventing the gas from cooling efficiently.

Zoom In Zoom Out Reset image size

Particles are included using a particle-in-cell method, where each particle represents a swarm of smaller particles. We refer to the particles that are directly integrated on the grid (the swarm) as super-particles and any member of that swarm as a sub-particle. All sub-particles have the same mass, density, and size, and are not allowed to evolve in time. The particles are connected to the grid using a cloud-in-cell (CIC) interpolation and prolongation, described as follows. Consider a 2D Cartesian grid as depicted in Figure 21. The intersections of lines represent grid nodes, and the yellow smiley face represents a super-particle. To interpolate values from the nodes a–d to the super-particle, the areas A–D are used for weights, where A is the weight from node a, etc. Prolongation is achieved using the same weighting scheme. In a 3D scheme, the nearest eight nodes are used for interpolation and prolongation in a volume-weighted scheme analogous to the area-weighted scheme shown Figure 21.

Zoom In Zoom Out Reset image size

Each super-particle has a position, velocity, and mass. For self-gravity calculations, the super-particle's contribution to the total grid mass is included by setting the super-particle's mass density to ρ_sp = m/V, where m is the super-particle's mass and V is the volume of a single cell on the grid. This mass density is then prolongated to the nearest eight cells, which are used when calculating the total gravitational potential field on the grid. The CIC algorithm with volume-weighted prolongation ensures that gravitational self-forces are absent. The temperature, density, and pressure of the gaseous fluid are tracked using CIC interpolation. In addition, the gaseous fluid's velocity is interpolated to the particle's position to calculate drag.

We use a drag algorithm that is heavily based on Boley & Durisen (2010, hereafter BD10), which includes the back-reaction of the particle onto the gas in addition to gas drag effects. As done in BD10, the Knudsen number is used to interpolate between the Epstein and Stokes drag regimes, where piecewise, analytic solutions for the velocity differences between the gas and the super-particle in the Epstein and Stokes regimes are calculated based on the Reynolds number. The velocity difference is then transformed to velocities relative to the grid using conservation of momentum. For reference, the Knudsen number Kn = μm_p/(2ρσ_gass), where μ is the mean molecular weight of the gas in proton masses m_p, s is the sub-particle size, ρ is the gas density, and σ_gas is the typical cross section of a gas particle, taken to be 10⁻¹⁶ cm². It represents the ratio of the mean-free-path of a particle embedded in a gas to the particle's size. The Reynolds number ${\rm Re}=3 \sqrt{\pi /8} \mathcal {M}/{\rm Kn}$ for Mach number $\mathcal {M}$, which reflects the ratio of inertial to viscous forces.

While most of the details of the algorithm are given in BD10, there are several important differences. First, BD10 used a nearest node approximation for coupling solid–gas interactions, where we use the CIC algorithm. After the change in gas momentum has been calculated, the total change is prolongated to the nearest eight grid nodes. Second, we include conservation of energy as well as momentum, while BD10 only included conservation of momentum. After the new momenta of the gas and super-particle are calculated, the change in total energy of the gas and super-particle is determined. The difference in total energy is then applied to the gas at the nearest eight cells using prolongation. Third, and perhaps most importantly, we include a method for handling drag whenever the drag equations become stiff relative to the hydrodynamics, i.e., the stopping time of the gas becomes comparable to or shorter than the Courant timescale.

The stopping time of a solid of size s and internal density ρ_s moving through a gas with density ρ and sound speed c_a is $t_s=\sqrt{\pi /8}(\rho _s s)/(\rho c_a)$. The terminal velocity of a solid moving through a gas is gt_s, where g is the local gravitational acceleration. Take the velocity difference between the gas and the super-particle after a full time step using the BD10 scheme as δv'. If the Courant step is greater than ρ_ssexp (− 2)/(ρc_a), then

$$\begin{equation} \delta {\bf v} = \delta {\bf v}^{\prime } w - {\bf g} t_s (1-w), \end{equation} \tag{ A17 }$$

where w = [ρ_ssexp (− 2)/(ρc_a)/t_CFL]², which was found to be an acceptable solution based on numerical experiments.

One potential problem with the above method for handling the stiff regime is that a super-particle will always reach the terminal velocity relative to the gas, even if the gas and the particle are in free fall. We attempt to circumvent this difficulty by using the pressure gradients ∇P in lieu of g. When the gas is in hydrostatic equilibrium, the two uses are equivalent, and using the pressure gradient allows the solids and gas to move together when pressure gradients are negligible.

The time evolution of solids' position and velocity is included by first applying a half kick (changing the velocity of the super-particle) during sourcing. During a half kick, the velocity is first updated without the drag term, and then the velocity is updated again by incorporating drag. After the grid variables are updated from full fluxing, a full drift term is applied (changing the position of the super-particle). A second half kick is then calculated using these updates, including the updated potential. This creates a second-order leap-frog scheme with time-step averaging, i.e., for a given time step, half a kick is applied using the previous potential field and the other half with the most recent potential field.

Tests including settling times in gaseous spheres demonstrate that particles reach terminal velocities as expected. Additional tests show that particles follow the flow of the gas when they are well-coupled and recover the correct stopping times when they are not perfectly coupled.

The radiation algorithm discussed in Appendix A.4 requires some estimate of the local opacity of the gas. The actual values used are independent of the MCRTFLD algorithm itself, and can be altered based on the problem at hand. One possibility is to use a Pollack et al. (1994) opacity table, which assumes some fixed composition of gas and solids, with a set distribution of solids sizes. With such tables, only the gas mass needs to be known, and the opacity can then be set according to the local gas density.

However, integrating super-particles directly into the fluid evolution gives rise to another possibility: using the super-particles themselves as the source of the opacity. As with a lookup table, such a strategy is not applicable to each problem, as many complexities can potentially arise. Nonetheless, there are multiple problems for which such a particle-based opacity method will be advantageous. One such example is the limit where there is no size evolution among sub-particles.

Consider for simplicity a region of some astrophysical object that is populated by dust grains that have the same composition and are of the same size. Further assume that each grain has a size s = 0.03 cm and has an internal density ρ_s = 3 g cm⁻³. A straightforward estimate for the mass opacity of such a population of dust is κ ≈ ξπs²/(4π/3sρ_s) = ξ3/(4sρ_s). For our system of dust grains, κ ≈ 8.3 cm² g⁻¹ of dust. We have assumed that the radiative efficiency of the grain ξ ∼ 1 due to the large size of the grain. A separate array can be used to track the gas density of the dust only, using the same prolongation technique needed for spatially evolving super-particles. With this information, the dust optical depth across a cell of width Δx is then Δτ = ρ_dustκΔx, assuming gray transfer.