The Planetary Accretion Shock. II. Grid of Postshock Entropies and Radiative Shock Efficiencies for Nonequilibrium Radiation Transport

The simulation domain extends from close to the planetary surface to almost out to the "accretion radius" ${R}_{\mathrm{acc}}$ (Bodenheimer et al. 2000; Paper I), near one-third of the Hill radius (${k}_{\mathrm{Lissauer}}=1/3$ Lissauer et al. 2009). A grid with 1000 uniformly spaced cells between the inner edge of the domain, ${r}_{\min }$, and ${r}_{\min }+0.5\,{R}_{{\rm{J}}}$ and 1000 geometrically stretched cells from ${r}_{\min }+0.5\,{R}_{{\rm{J}}}$ to the outer edge of the domain, ${r}_{\max }$, was found to yield good results for several parameter combinations. For the other cases, an increase to 2000+2000 cells sufficed to obtain smooth profiles (e.g., in the local accretion rate $\dot{M}(r)=4\pi {r}^{2}\rho (v)v(r)$ below the shock). We verified that increasing the resolution does not produce different results except for sharpening the Zel'dovich spike (see below) and changing its peak value. Thus there are no qualitative consequences on the pre- or postshock structure.

To avoid numerical issues at early times, we have revised the initial setup. We do not use a hydrostatic atmosphere as was done in Paper I, but rather begin directly with an accretion profile in the density ρ and velocity v. The initial radiation and gas temperatures are set to the nebula temperature ${T}_{\mathrm{neb}}=150$ K (Pollack et al. 1994) at the outer edge ${r}_{\max }$ and increase linearly up to $1.1{T}_{\mathrm{neb}}$ at ${r}_{\min }$. The radiation and gas temperatures are set equal, and the initial pressure is obtained from ρ and ${T}_{\mathrm{gas}}$.

We use a CFL number of 0.8 (we verified that CFL = 0.4 gives identical profiles), and use the tvdlf Riemann solver along with the MinMod slope limiter. This somewhat diffusive combination ensures numerical stability while not significantly influencing the outcome.

The inner wall is reflective and the radiation flux there is set to zero. We let gas fall in freefall from the outer edge under the action of the central potential, enforcing an accretion rate at the outer edge only. The gas piles up at the bottom of the domain, forming a shock that defines the surface of the planet at ${r}_{\mathrm{shock}}={R}_{{\rm{p}}}$ (both variables are used interchangeably throughout, with ${{r}_{\mathrm{shock}}}^{-}$ (${{r}_{\mathrm{shock}}}^{+}$) denoting the downstream (upstream) location). This shock surface often but not always moves outward over time, although at a rate that is slow even compared to the strongly subsonic postshock velocity. The deepest pressure in the atmosphere we simulate is greater than the postshock (ram) pressure by one to several orders of magnitude, depending on the amount of accumulated mass and the local gravitational acceleration $g={{GM}}_{{\rm{p}}}/{r}^{2}$, where r is the radial coordinate in spherical geometry, and G is the gravitational constant.

While we use a fully time-dependent code, our simulations represent steady-state snapshots in the formation-parameter space, as discussed in Paper I. For reference, the profiles for the runs presented here needed on the order of 5 × 10⁶ s to come into equilibrium. This is entirely negligible compared to the timescale on which the protoplanet grows, which is on the order of 10⁴–10⁶ yr. (Even after ${t}_{\mathrm{stop}}=2\times {10}^{7}\,{\rm{s}}$, our usual stopping time, some starting structures were still visible at a position r deep below the shock for which ${t}_{\mathrm{stop}}\lesssim -{\int }_{r}^{{r}_{\mathrm{shock}}}(1/v){dx}$ in at least some simulations; these regions are not of interest here, however.)

A significant improvement since Paper I is the change from equilibrium to nonequilibrium (or two-temperature, 2-T) FLD in the radiation transport routine Makemake (Kuiper et al. 2010, R. Kuiper et al. in preparation) In this approach, the gas and radiation temperatures are not enforced to be equal, and the full energy equation reads (see also Mihalas & Mihalas 1984; Kley 1989; Turner & Stone 2001; R. Kuiper et al. 2010; Commerçon et al. 2011b; Klassen et al. 2014)

$$\begin{eqnarray}&&\displaystyle \frac{\partial \left({E}_{\mathrm{kin}}+{E}_{\mathrm{int}}\right)}{\partial t}+{\rm{\nabla }}\left(\left[{E}_{\mathrm{kin}}+{E}_{\mathrm{int}}+P\right]v\right)=-{\rm{\Lambda }}-\rho {gv},\end{eqnarray} \tag{ 3a }$$

$$\begin{eqnarray}&&\,\displaystyle \frac{\partial {E}_{\mathrm{rad}}}{\partial t}+{\rm{\nabla }}{F}_{\mathrm{rad}}=+{\rm{\Lambda }},\end{eqnarray} \tag{ 3b }$$

omitting the radiation pressure because it is negligible in this problem, and with ${F}_{\mathrm{rad}}$ the radiative flux. In FLD, one writes (Paper I)

$$\begin{eqnarray}&&\,{F}_{\mathrm{rad}}=-{D}_{{\rm{F}}}{\boldsymbol{\nabla }}{E}_{\mathrm{rad}},\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&{D}_{{\rm{F}}}\equiv \displaystyle \frac{\lambda (R)c}{{\kappa }_{{\rm{R}}}\rho },\end{eqnarray} \tag{ 4b }$$

$$\begin{eqnarray}&&\,R\equiv \displaystyle \frac{\parallel {\boldsymbol{\nabla }}\mathrm{ln}{E}_{\mathrm{rad}}\parallel }{{\kappa }_{{\rm{R}}}\,\rho },\end{eqnarray} \tag{ 4c }$$

where ${D}_{{\rm{F}}}$ is the diffusion coefficient, λ the flux limiter (see Section 5.1 below), ${\kappa }_{{\rm{R}}}$ the Rosseland mean, and R the "local radiation quantity" (see Paper I). Thus the radiative flux is assumed to be colinear with and opposite in direction to the gradient of radiation energy density.⁷ In Equation (3), the combined cooling and heating or exchange term is

$$\begin{eqnarray}&&{\rm{\Lambda }}\equiv c{\kappa }_{{\rm{P}}}\rho \displaystyle \frac{4\pi S}{c}-c{\kappa }_{{\rm{E}}}\rho {E}_{\mathrm{rad}}\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&=\,c{\kappa }_{{\rm{P}}}\rho \left(a{{T}_{\mathrm{gas}}}^{4}-{E}_{\mathrm{rad}}\right).\end{eqnarray} \tag{ 5b }$$

We use the second line because we take the source function S to be the Planck function $B=\sigma {{T}_{\mathrm{gas}}}^{4}/\pi$, with $4\sigma ={ac}$ ($\sigma$ being the Stefan–Boltzmann constant and a and c the radiation constant and speed of light). This assumes local thermodynamic equilibrium. Our gray approximation suffices to make ${\kappa }_{{\rm{P}}}$ and ${\kappa }_{{\rm{E}}}$, the Planck (blackbody-weighted) and ${E}_{\mathrm{rad}}$-weighted mean opacities, respectively, be the same. Equations 3(a) and (b) are followed separately to prevent the Λ terms from canceling. This system states that the material (gas or dust, assumed here to have the same temperature) is losing energy at the rate $c{\kappa }_{{\rm{P}}}\rho a{{T}_{\mathrm{gas}}}^{4}$ but absorbing energy (photons) at the rate $c{\kappa }_{{\rm{E}}}\rho {E}_{\mathrm{rad}}$, and conversely for the radiation.

To solve the energy Equation 3(b) in the radiation transport substep, the radiative flux ${F}_{\mathrm{rad}}$ is replaced by its expression from Equation (4), so that the FLD approach in fact does not naturally yield ${F}_{\mathrm{rad}}$ (nor thus the luminosity). While it is possible at a given time to calculate it from the output quantities $\rho (r)$, ${E}_{\mathrm{rad}}(r)$, ${\kappa }_{{\rm{R}}}(r)$ by Equation (4), this is approximate because in our implicit approach the diffusion coefficient ${D}_{{\rm{F}}}$ (Equation 4(b)) is computed from the quantities at the current time while the ${\rm{\nabla }}{E}_{\mathrm{rad}}$ factor is written with ${E}_{\mathrm{rad}}$ at the following timestep. To avoid any numerical noise and inaccuracy, we therefore store ${D}_{{\rm{F}}}$ before the FLD step and combine it with ${E}_{\mathrm{rad}}$ after, thus exactly reconstructing the flux as it is effectively computed (albeit otherwise not explicitly). This was also done in Paper I but was not reported there.

We use the same outer boundary condition for ${E}_{\mathrm{rad}}$ as in Paper I, but verified that changing the outer temperature boundary condition (in particular to a Dirichlet boundary condition) does not affect the temperature or luminosity structure except at the largest radii.

With respect to Paper I, we now always use the dust opacities of Semenov et al. (2003), taking by default their simplest model for the dust grains, the "normal homogeneous spheres" (nrm_h_s). It features more evaporation transitions than in Bell & Lin (1994), which we used originally. We slightly modified the opacity routine⁸ to leave the evaporation of the most refractory component to be handled in the main code and to avoid adding the contribution from the Helling et al. (2000) gas opacities.

Instead, the gas opacity is provided by Malygin et al. (2014), and we take their one-temperature Planck mean, as opposed to the 2-T Planck mean (see their Equation (4)). We found that it is crucial for numerical stability, especially toward higher constant μ and γ, to evaluate the single-temperature ${\kappa }_{{\rm{P}}}$ at the radiation temperature ${T}_{\mathrm{rad}}\equiv {\left({E}_{\mathrm{rad}}/a\right)}^{1/4}$ and not at the gas temperature ${T}_{\mathrm{gas}}$. Fortunately, it is also justified. Indeed, looking at ${\kappa }_{{\rm{P}}}({T}_{\mathrm{rad}},{T}_{\mathrm{gas}})$ for fixed densities, one generally incurs a smaller mistake when using ${\kappa }_{{\rm{P}}}\approx {\kappa }_{{\rm{P}}}({T}_{\mathrm{rad}},{T}_{\mathrm{rad}})$ than with ${\kappa }_{{\rm{P}}}\approx {\kappa }_{{\rm{P}}}({T}_{\mathrm{gas}},{T}_{\mathrm{gas}})$. We also evaluate the dust opacities at ${T}_{\mathrm{rad}}$ but note that this barely makes a difference because we will find that the radiation and gas are always well coupled (${T}_{\mathrm{gas}}={T}_{\mathrm{rad}}$) in the region where the dust opacity matters (${T}_{\mathrm{gas}}\lesssim 1500$ K).

We recall that the Planck opacity is relevant for the coupling between dust/gas and radiation, whereas the Rosseland mean determines whether the radiation can stream freely or has to diffuse. This is discussed in more detail in Section 7.2.

The Planck opacity ${\kappa }_{{\rm{P}}}$ is shown in Figure 2 and is compared to the Rosseland mean ${\kappa }_{{\rm{R}}}$. In the low-temperature, dust-dominated regime ${\kappa }_{{\rm{P}}}$ is similar to ${\kappa }_{{\rm{R}}}$ with ${\kappa }_{{\rm{P}}}\sim {\kappa }_{{\rm{R}}}\,\sim 1$–10 cm² g⁻¹. For comparison, the different models of Semenov et al. (2003) are shown, with the exception of the porous models. We note that the fluffiness or porosity of dust aggregates during their growth remains an open question (see Cuzzi et al. 2014; Kataoka et al. 2014; Kirchschlager et al. 2019; Tazaki et al. 2019, and references in these works). In any case, the model curves of Semenov et al. (2003) for porous composite particles are not markedly different from the ones shown in Figure 2.

Zoom In Zoom Out Reset image size

While the Rosseland mean opacity past dust evaporation (near $T\approx 1500$ K) drops to ${\kappa }_{{\rm{R}}}\approx {10}^{-2}$ cm² g⁻¹ (see Figure 1 of Paper I), the Planck mean does not drop much below ${\kappa }_{{\rm{P}}}\,\approx 1$ cm² g⁻¹ and even increases (between a few thousand to 10⁴ K) to ${\kappa }_{{\rm{P}}}={10}^{2}$–10⁴ cm² g⁻¹, depending on the density.

Note that the Bell & Lin (1994) Rosseland mean opacities are three to six orders of magnitude (!) lower than the Plank average above the dust destruction temperature. This implies that studies using nonequilibrium radiation transport with the Bell & Lin (1994) Rosseland mean as their Plank opacity are effectively assuming much less coupling between opacity carrier (dust or gas) and radiation. As Malygin et al. (2014) found out, a similar word of caution applies to the gas opacities of Helling et al. (2000), which are included in Semenov et al. (2003). Whether using these low opacities would actually lead to a strong disequilibrium between the matter and radiation will also depend on the local density and velocity, however, as briefly discussed in Section 7.2.

For our choice of macrophysical parameters, the shock temperature is always above 1000 K. In Figure 3, simulations with a given $(\dot{M},{M}_{{\rm{p}}})$ but different ${r}_{\min }$ (symbol size) are seen to lead to the same shock temperature at a given radius ${r}_{\mathrm{shock}}\,={R}_{{\rm{p}}}$. This confirms that even though the postshock density structures differ, the shock properties do not depend on our placement of the inner boundary, thus supporting the robustness of our results. We have also verified that varying other numerical parameters such as the resolution or the tolerance in the FLD solver step (R. Kuiper et al. 2019, in preparation) does not modify the simulation outcomes.

An interesting result is that in all cases, the postshock gas and radiation (not shown) are in equilibrium, with ${T}_{\mathrm{gas}}={T}_{\mathrm{rad}}$ to better than 1%, so that one can speak of a single temperature. The same applies to the preshock temperature. (For other choices of γ the equality is less strict but still within a few percent.) In Figure 3(a), the values we obtained are compared to the analytical shock estimate from Equation (27) of Paper I

$$\begin{eqnarray}&&\sigma {T}_{\mathrm{shock}}^{4}=\displaystyle \frac{{\eta }^{\mathrm{kin}}}{4{\rm{\Delta }}{f}_{\mathrm{red}}}\displaystyle \frac{{\rho }_{\mathrm{ff}}{{v}_{\mathrm{shock}}}^{3}}{2}\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&\,\approx \,\displaystyle \frac{{{GM}}_{{\rm{p}}}\dot{M}}{16\pi {{R}_{{\rm{p}}}}^{3}},\end{eqnarray} \tag{ 6b }$$

where ${\eta }^{\mathrm{kin}}\equiv {\rm{\Delta }}{F}_{\mathrm{rad}}/(0.5\rho {v}^{3})\approx {L}_{\mathrm{acc}}/({{GM}}_{{\rm{p}}}\dot{M}/{R}_{{\rm{p}}})$ is the normalized jump in luminosity and ${\rm{\Delta }}{f}_{\mathrm{red}}\equiv {f}_{\mathrm{red}}({{r}_{\mathrm{shock}}}^{+})\,-{f}_{\mathrm{red}}({{r}_{\mathrm{shock}}}^{-})$ the jump in ${f}_{\mathrm{red}}$ at the shock. This equation is revisited in Section 5.3. The freefall density ${\rho }_{\mathrm{ff}}=\dot{M}/(4\pi {r}^{2}{v}_{\mathrm{ff}})$ is evaluated ahead of the shock, with ${v}_{\mathrm{ff}}=\sqrt{2{{GM}}_{{\rm{p}}}/r}$ the approximate freefall velocity Paper I. The reduced flux

$$\begin{eqnarray}&&{f}_{\mathrm{red}}\equiv {F}_{\mathrm{rad}}/{{cE}}_{\mathrm{rad}}\end{eqnarray} \tag{ 7 }$$

is also termed the "streaming factor" (Kley 1989) because it indicates to what extent the radiation is freely streaming (${f}_{\mathrm{red}}\to 1$) or diffusing (${f}_{\mathrm{red}}\to 0$ ). In usual shock terminology, free-streaming regions are called "transmissive" (see Figure 8 of Vaytet et al. 2013b; Drake 2006; we recall in passing that Equation 6(a) is valid for ${L}_{\mathrm{int}}\gtrless {L}_{\mathrm{acc}}$.) The simulations have preshock Mach numbers ${\mathscr{M}}\approx 7$–35 (Figure 4) and therefore ${\eta }^{\mathrm{kin}}=1$ (used in going from Equations 6(a) to (b)) because this is above ${\mathscr{M}}\approx 2.5$ (see the thin line in Figure 4; Paper I). Moreover, we find that ${\rm{\Delta }}{f}_{\mathrm{red}}\approx 1$ for a large part of the parameter space, i.e., the downstream regions are in the diffusion limit with ${f}_{\mathrm{red}}({{r}_{\mathrm{shock}}}^{-})\approx 0.03$–0.05 (equivalently, a flux limiter λ ≈ 1/3), while the preshock region is in the free-streaming regime, with ${f}_{\mathrm{red}}({{r}_{\mathrm{shock}}}^{+})\approx 1$. Thus the shock is a "thick–thin" shock in the classification of Drake (2006).

As a consequence of this, Equation 6(b) holds very accurately for almost all simulations; we are almost always in the limit of Equation (28) of Paper I, discussed as Equation (33) below. There is an exception to this, namely for the lowest mass (1.3 M_J) at a low accretion rate (${10}^{-3}\,{M}_{\oplus }\,{\mathrm{yr}}^{-1}$) and toward larger radii $({R}_{{\rm{p}}}\gtrsim 2.5\,{R}_{{\rm{J}}})$. In this case, the postshock temperature is higher than predicted from Equation 6(b) by ΔT ≈ 50 K. This is because ${f}_{\mathrm{red}}$ is lower ahead of the shock, with, e.g., ${{f}_{\mathrm{red}}}^{+}\,\equiv {f}_{\mathrm{red}}({{r}_{\mathrm{shock}}}^{+})=0.65$ for the ${T}_{\mathrm{shock}}\approx 1100$ K case. Indeed, because ${F}_{\mathrm{rad}}\propto {f}_{\mathrm{red}}{T}_{\mathrm{rad}}^{4}$ and we find ${T}_{\mathrm{gas}}={T}_{\mathrm{rad}}$, a lower ${f}_{\mathrm{red}}$ at a location requires a higher temperature for the same radiation flux (equal to the kinetic energy flux) to flow through that location. As discussed in Paper I, this smaller ${f}_{\mathrm{red}}$ can also be pictured as a slower effective speed of light, so that ${E}_{\mathrm{rad}}$ must increase so as to have the same ${F}_{\mathrm{rad}}={c}_{\mathrm{eff}}{E}_{\mathrm{rad}}$. We return to these points in Section 5.3 and show that another effect is also at play. Note, however, that this is the gas temperature but not the effective temperature, with only the latter setting the spectral shape.

Figure 3(b) shows the luminosity at the outer edge of the grid, which is very nearly equal to the luminosity at the Hill sphere, as we show in Section 6.1.2. The Mach numbers we find here are all ${\mathscr{M}}\gt 2.5$, so that essentially the entire kinetic energy is converted into radiation (see the pale gray curve in Figure 4). Because for the specific choice of $(\mu =2.29,\gamma =1.44)$, the decrease in L from the shock to ${R}_{\mathrm{Hill}}$ is insignificant (see Section 6.1), the entire kinetic energy is transformed into visible radiation, according to the usual expression for the shock luminosity,

$$\begin{eqnarray}&&{L}_{\mathrm{acc}}=\displaystyle \frac{{{GM}}_{{\rm{p}}}\dot{M}}{{R}_{{\rm{p}}}}.\end{eqnarray} \tag{ 8 }$$

This equation is shown as solid lines and seen to match very well (note that the symbols are almost smaller than the line), even though we neglected the finite "accretion radius" ${R}_{\mathrm{acc}}$. This radius is much larger than the planet radius, however, because we are considering the detached phase during giant planet formation (Mordasini et al. 2012b).

We also produced a grid with (μ = 1.1, γ = 1.1) (not shown), which increases the contrast with the present situation. One example is discussed in Section 6.1 below. Over the grid, however, the shock temperatures were the same, as one would expect from Equation (6) because they do not explicitly depend on the EOS. The situation could be different for an ideal but nonperfect EOS, but due to the potential sink of energy (dissociation and ionization). As for the luminosity at the Hill radius, it was different from the $\mu =1.23$ case and was lower by at most ≈10% compared to the immediate shock upstream luminosity. The luminosity profile is discussed in Section 6.1.

Figure 3(c) displays the Rosseland mean optical depth from the shock to the outer edge of the computational domain, ${\rm{\Delta }}{\tau }_{{\rm{R}}}=-{\int }_{{r}_{\max }}^{{r}_{\mathrm{shock}}}\rho {\kappa }_{{\rm{R}}}\,{dr}$. Especially for the lower masses, the contribution from the outer layers is actually significant despite their low density, and ${\rm{\Delta }}{\tau }_{{\rm{R}}}$ does depend on the choice of the domain size (here, ${r}_{\max }=0.7{R}_{\mathrm{acc}}$ as in Paper I). We find that ${\rm{\Delta }}{\tau }_{{\rm{R}}}\approx 0.2$–5, increasing with accretion rate and decreasing with mass. This is qualitatively as expected from Equation (24) of Paper I

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{\tau }_{{\rm{R}}} & \sim & 3\left(\displaystyle \frac{{\kappa }_{{\rm{R}}}}{1\,{{\rm{c}}{\rm{m}}}^{2}\,{{\rm{g}}}^{-1}}\right)\,\left(\displaystyle \frac{\dot{M}}{{10}^{-2}\,{M}_{\oplus }\,{{\rm{y}}{\rm{r}}}^{-1}}\right)\\ & & \times \sqrt{\left(\displaystyle \frac{1\,{M}_{{\rm{J}}}}{{M}_{{\rm{p}}}}\right)\,\left(\displaystyle \frac{2\,{R}_{{\rm{J}}}}{{r}_{{\rm{s}}{\rm{h}}{\rm{o}}{\rm{c}}{\rm{k}}}}\right)},\end{array}\end{eqnarray} \tag{ 9 }$$

shown as dotted lines in Figure 3(c). The agreement with the nominal models is not too rough (to at worse 1 dex) considering that we took a constant ${\kappa }_{{\rm{R}}}=1$ cm² g⁻¹ in Equation (9) in this estimate for all simulations. In the radial profiles, the maximum value of the actual ${\kappa }_{{\rm{R}}}(r)$ is near 7(3) cm² g⁻¹ for $\dot{M}={10}^{-3}$(10⁻²) ${M}_{\oplus }$ yr⁻¹. However, contrary to what was explained in Paper I, what is relevant is the immediate preshock opacity, as we show in Section 5.

The significance of Figure 3(c) can only be appreciated in conjunction with panels (a) and (b). It shows that the total optical depth, at least up to moderate depths of ${\rm{\Delta }}{\tau }_{{\rm{R}}}\sim 10$, does not set the shock temperature or the shock luminosity (which is nearly identical to the Hill-radius luminosity in the present case). In fact, the deviations of ${T}_{\mathrm{shock}}$ from Equation (6) do not occur at the highest optical depths. (What causes these deviations is explained in Section 5.) Furthermore, when, as is often the case, the layers between the ${\rm{\Delta }}{\tau }_{{\rm{R}}}\sim 1$ surface and the shock are not sufficiently diffusive (${f}_{\mathrm{red}}\lesssim 0.1$), it does not hold that ${T}_{\mathrm{shock}}^{4}\approx \tfrac{3}{4}\left({\rm{\Delta }}{\tau }_{{\rm{R}}}(r)+\tfrac{2}{3}\right){T}_{\mathrm{eff}}^{4}$, where ${T}_{\mathrm{eff}}$ is the temperature at the radius where ${\rm{\Delta }}{\tau }_{{\rm{R}}}=2/3$, the photosphere, and where ${\rm{\Delta }}{\tau }_{{\rm{R}}}(r)=-{\int }_{{r}_{\max }}^{r}\rho {\kappa }_{{\rm{R}}}\,{dr}^{\prime}$ is the optical depth measured from the outer radius inward. The temperature rise from the photosphere to the shock is much higher because the optical depth increases more slowly when the radiation is less diffusive, i.e., when $d\tau /{dr}={\kappa }_{{\rm{R}}}\rho$ is small.

Finally, we note that there is nothing special about the ${\rm{\Delta }}{\tau }_{{\rm{R}}}\sim 1$ surface. As Figure 3(c) shows, most simulations at $\dot{M}={10}^{-3}\,{M}_{\oplus }$ yr⁻¹ remain optically thin or barely optically thick above the shock, but their temperature structure (not shown) is qualitatively identical to cases with an optically thick accretion flow. In any case, ${\rm{\Delta }}{\tau }_{{\rm{R}}}$ is not well defined because it depends on the outer integration radius. We return to considerations of shock temperature in Section 5.3.

When calculating a planet structure for given $(\dot{M},{M}_{{\rm{p}}},{R}_{{\rm{p}}})$, only the postshock $({P}_{\mathrm{post}},{T}_{\mathrm{shock}})$ point is used in an approach equivalent to the one used in mesa by Berardo et al. (2017) and Berardo & Cumming (2017). This yields the density ρ; from this and $\dot{M}$, the boundary condition on the velocity v and thus the luminosity profile L(r) in the postshock region follow (Berardo et al. 2017). When only the global energetics are followed (see the review in Section 2.1 of Berardo et al. 2017), we argued in Paper I that one should use ${\eta }^{\mathrm{phys}}$ and not ${\eta }^{\mathrm{kin}}$.

We show in Figure 4 the "global physical efficiency" ${\eta }^{\mathrm{phys}}$ of the shock against the preshock Mach number. The quantity ${\eta }^{\mathrm{phys}}$ is measured as (Paper I, their Equation (18))

$$\begin{eqnarray}&&{\eta }^{\mathrm{phys}}\equiv \displaystyle \frac{\dot{E}({r}_{\max })-\dot{E}({{r}_{\mathrm{shock}}}^{-})}{\dot{E}({r}_{\max })},\end{eqnarray} \tag{ 10 }$$

where ${{r}_{\mathrm{shock}}}^{-}$ is immediately downstream of the shock. The material-energy flow rate is defined as

$$\begin{eqnarray}&&\dot{E}(r)\equiv -| \dot{M}| \left[{e}_{\mathrm{kin}}(r)+h(r)+{\rm{\Delta }}{\rm{\Phi }}(r,{r}_{\mathrm{shock}})\right],\end{eqnarray} \tag{ 11 }$$

where ${e}_{\mathrm{kin}}=\tfrac{1}{2}{v}^{2}$, ${e}_{\mathrm{int}}$, and $h={e}_{\mathrm{int}}+P/\rho$ are the kinetic energy, internal energy density, and the enthalpy per unit mass, respectively, and Φ is the external potential. The ${\rm{\Delta }}{\rm{\Phi }}$ term in Equation (11) accounts for the work done by the potential on the gas down to the shock, with the potential difference from r₀ to r given by

$$\begin{eqnarray}&&{\rm{\Delta }}{\rm{\Phi }}(r,{r}_{0})=-{{GM}}_{{\rm{p}}}\left(\displaystyle \frac{1}{r}-\displaystyle \frac{1}{{r}_{0}}\right).\end{eqnarray} \tag{ 12 }$$

We found that for all simulations here, using the tenth cell below the shock as the postshock location was a robust prescription. The shock itself was identified by the ${dv}/{dr}\lt 0$ and ${\rm{\Delta }}P/\min (P)\gt 5$ criterion in Appendix B of Mignone et al. (2012).

The efficiency ${\eta }^{\mathrm{phys}}$ measured from Equation (10) is compared in Figure 4 as a function of the Mach number to the analytical result in the isothermal limit (Paper I, their Equation (36)),

$$\begin{eqnarray}&&{\eta }_{{\rm{isoth}}}^{\mathrm{phys}}={\eta }_{{\rm{isoth}}}^{\mathrm{kin}}\times {\left(1+\displaystyle \frac{2}{\gamma -1}\displaystyle \frac{1}{{{\mathscr{M}}}^{2}}\right)}^{-1}\end{eqnarray} \tag{ 13a }$$

$$\begin{eqnarray}&&\,=\,\displaystyle \frac{\left({\gamma }^{2}{{\mathscr{M}}}^{4}-1\right)(\gamma -1)}{{\gamma }^{2}{{\mathscr{M}}}^{2}\left[(\gamma -1){{\mathscr{M}}}^{2}+2\right]},\end{eqnarray} \tag{ 13b }$$

with the isothermal "kinetic efficiency" given by (Commerçon et al. 2011a)

$$\begin{eqnarray}&&{\eta }_{{\rm{isoth}}}^{\mathrm{kin}}=1-\displaystyle \frac{1}{{\gamma }^{2}{{\mathscr{M}}}^{4}}.\end{eqnarray} \tag{ 14 }$$

This definition was derived from energy conservation for the 1-T case but it is also meaningful here because shocks are isothermal in the radiation temperature, and we find that the gas and radiation are tightly coupled.

As expected, we find essentially perfect agreement between the measured (Equation (10)) and the theoretical (Equation (13)) efficiencies. This reflects both energy conservation by our radiation-hydrodynamical code and the isothermality—i.e., supercriticality—of the shock.

We also computed the kinetic efficiency,

$$\begin{eqnarray}&&{\eta }^{\mathrm{kin}}\equiv \displaystyle \frac{{\rm{\Delta }}{F}_{\mathrm{rad}}}{\tfrac{1}{2}{\rho }_{+}{{v}_{+}}^{3}},\end{eqnarray} \tag{ 15 }$$

where ${\rm{\Delta }}{F}_{\mathrm{rad}}$ is the jump in radiative flux at the shock and ρ₊ and v₊ are the density and velocity upstream of the shock. This measured ${\eta }^{\mathrm{kin}}$ was found to match Equation (14). As mentioned above, because ${\mathscr{M}}\gtrsim 2.5$ and the shock is isothermal, we have ${\eta }^{\mathrm{kin}}\approx 100$%. Locally at the shock, the whole incoming kinetic energy is thus converted into radiation.

Figure 5 shows the global physical efficiency ${\eta }^{\mathrm{phys}}$ but now as an explicit function of the (possibly observable, macrophysical) formation parameters $(\dot{M},{M}_{{\rm{p}}},{R}_{{\rm{p}}})$. The efficiencies range from 97% at high accretion rate to almost 100% at low $\dot{M}$, increasing both with decreasing radius and increasing mass. The precise range depends on the assumed EOS, but qualitatively our results should be robust. Typically, in core-accretion formation, the radius decreases as the mass grows in the detached phase.⁹ Because the gas accretion onto a planet usually slows down with time, at least in the single-embryo-per-disk simulations of Mordasini et al. (2012b), our simulations clearly suggest that the efficiency ${\eta }^{\mathrm{phys}}$ increases over time. Thus, the accretion of the outer layers is associated with less energy recycling (Paper I) in the accretion flow, and a greater net fraction of the kinetic energy escapes the system.

However, as discussed in Mordasini (2013), there could be a self-amplifying memory effect: an ${\eta }^{\mathrm{phys}}$ that is small early after detachment should lead to the accretion of "hot" (high-entropy) material into the planet. If the planetary accretion timescale is much shorter than its global cooling timescale—i.e., ${t}_{\mathrm{acc}}\equiv {M}_{{\rm{p}}}/\dot{M}\ll {t}_{\mathrm{KH}}\equiv {{GM}}_{{\rm{p}}}^{2}/({R}_{{\rm{p}}}{L}_{{\rm{p}}})$, where ${t}_{\mathrm{acc}}$ and ${t}_{\mathrm{KH}}$ are the accretion and Kelvin–Helmholtz timescales, respectively, with L the luminosity at the radiative–convective boundary (Berardo et al. 2017)—, a large planet radius should ensue. This large radius would in turn keep ${\eta }^{\mathrm{phys}}$ small, so that the planet would always remain in a regime where it is accreting rather high-entropy material. Ultimately, this would lead to a hot (or at least warm) start. This scenario should be tested with dedicated evolutionary simulations coupling the shock efficiency self-consistently with the interior structure. It might also be instructive to compare this with accretion in the context of star formation, where a similar phenomenon is seen (Hosokawa & Omukai 2009; Hosokawa et al. 2013; Kuiper & Yorke 2013). Care should be taken to distinguish in the analysis the local accretion and cooling times of the material immediately below the accretion shock from the global ones ${t}_{\mathrm{acc}}$ and ${t}_{\mathrm{KH}}$ (i.e., for the planet as a whole).

Our detailed calculations of the accretion shock are meant to serve as outer boundary conditions for calculating the structure of accreting planets as in Mordasini et al. (2012b), Berardo et al. (2017), Berardo & Cumming (2017), and Cumming et al. (2018). Because these planet structure calculations always use a full EOS including chemical reactions (dissociation and ionization), we use this to compute the entropy corresponding to the shock temperatures and pressures we obtained above. This is formally not self-consistent given that we assumed a perfect EOS for the shock simulations. However, because ${T}_{\mathrm{shock}}$ should be independent of the EOS, as we argue in Section 6.1, the entropy values are possibly realistic. In any event, they will serve as a comparison point for simulations with a realistic EOS (G.-D. Marleau et al. 2019, in preparation).

To calculate the entropy, we can use the ideal-gas form $P=\rho /(\mu {m}_{{\rm{H}}}){k}_{{\rm{B}}}T$ (but with variable μ) because degeneracy effects start being relevant only at a conservative limit of ρ ∼ 10⁻² g cm⁻³ (see Figure 1 of SCvH), several orders of magnitude above even the highest postshock densities ${\rho }_{\mathrm{post}}\propto {\rho }_{\mathrm{pre}}{{\mathscr{M}}}^{2}$, with ${\rho }_{\mathrm{pre}}\sim {10}^{-13}$–10⁻¹⁰ g cm⁻³ (see Figure 4 of Paper I) and ${\mathscr{M}}\ll 100$, so that ${\rho }_{\mathrm{post}}\ll {10}^{-6}$ g cm⁻³. We use the Saha equation and Sackur–Tetrode formula as implemented¹⁰ in Berardo et al. (2017). The entropy zero-point (see Appendix B of Marleau & Cumming 2014 and footnote 2 of Mordasini et al. 2017) is the same as in the published version of SCvH and mesa (Paxton et al. 2011, 2013, 2015, 2018), and thus the entropy values reported here are higher by $(1-Y)\mathrm{ln}2=0.52\,{k}_{{\rm{B}}}\,{{m}_{{\rm{H}}}}^{-1}$, where Y = 0.25 is the helium mass fraction, than those in, e.g., Mordasini et al. (2017). The zero-point of the entropy is not physically meaningful but does have to be taken into account when comparing entropies from different works.

These effective postshock entropies ${s}_{\mathrm{post}}$ are shown in Figure 5. For the range of shock positions ${r}_{\mathrm{shock}}\approx 1.5$–3 R_J and masses ${M}_{{\rm{p}}}=1.3$–10 M_J shown, the entropies for $\dot{M}\,={10}^{-3}\,{M}_{\oplus }\,{{\rm{yr}}}^{-1}$ are mostly around ${s}_{\mathrm{post}}\approx 12$–14 but go up to ${s}_{\mathrm{post}}\approx 19$ (dropping, also in the following, the usual units of ${k}_{{\rm{B}}}\,{{m}_{{\rm{H}}}}^{-1}$ for clarity). For $\dot{M}={10}^{-2}\,{M}_{\oplus }\,{{\rm{yr}}}^{-1}$, the whole range ${s}_{\mathrm{post}}\approx 13$–20 is covered. These values are high compared to the postformation entropy of planets, which is at most around 10–14 ${k}_{{\rm{B}}}\,{{m}_{{\rm{H}}}}^{-1}$ according to current, though not definitive, predictions (Mordasini 2013; Berardo & Cumming 2017; Berardo et al. 2017; Mordasini et al. 2017). However, we caution and emphasize that this postshock entropy is not the same as the entropy below the postshock settling layer; this latter quantity is most likely the one most relevant in setting the entropy of the planet as it accretes.

Before calculating the (non-self-consistent) postshock entropy analytically, we briefly discuss the ram pressure,

$$\begin{eqnarray}&&{P}_{\mathrm{ram}}=\rho {v}_{\mathrm{ff}}^{2}=\displaystyle \frac{\dot{M}}{4\pi {{R}_{{\rm{p}}}}^{2}}{v}_{\mathrm{ff}}=\displaystyle \frac{\sqrt{2G}}{4\pi }\displaystyle \frac{{M}_{{\rm{p}}}^{1/2}\dot{M}}{{{R}_{{\rm{p}}}}^{5/2}}.\end{eqnarray} \tag{ 16 }$$

We inserted the expressions for freefall from infinity and find that Equation (16) holds very well for all simulations, even for those for which ${\rm{\Delta }}{f}_{\mathrm{red}}\ne 1$. At high shock temperatures, the small pressure buildup ahead of the shock slows the gas down slightly, making Equation (16) less accurate by at most 3.5% for the range of parameters shown. Given the relatively weak (logarithmic, with a small prefactor) dependence of the entropy s(P, T) on P, outside of dissociation or ionization regions, this will not be an important source of inaccuracy. The ram pressure varies from ${P}_{\mathrm{ram}}\approx {10}^{-4}$ bar to 0.2 bar for the range of parameters discussed here.

We compare the postshock entropies ${s}_{\mathrm{post}}$ in Figure 5 based on the actual T to ${s}_{\mathrm{post}}$ using as input $T={T}_{\mathrm{shock}}$ from Equation 6(b), i.e., taking ${\eta }^{\mathrm{kin}}={\rm{\Delta }}{f}_{\mathrm{red}}=1$ for all simulations. The match is very good, which reflects the overall good match of temperature, on which the entropy depends only logarithmically.

To derive the slope of the temperature throughout the infalling gas, let us begin with the general relationships

$$\begin{eqnarray}&&{{acT}}_{\mathrm{rad}}^{4}=\displaystyle \frac{L}{4\pi {r}^{2}}\displaystyle \frac{1}{{f}_{\mathrm{red}}}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\,{f}_{\mathrm{red}}=R\lambda (R),\end{eqnarray} \tag{ 18 }$$

with ${{aT}}_{\mathrm{rad}}^{4}\equiv {E}_{\mathrm{rad}}$. The first equation is nothing but a rewriting of the definition of ${f}_{\mathrm{red}}$, and the second follows from Equation (17) and the definition of the radiation quantity R, given by

$$\begin{eqnarray}&&R\equiv \displaystyle \frac{1}{{\kappa }_{{\rm{R}}}\rho }\left|\displaystyle \frac{d\mathrm{ln}{E}_{\mathrm{rad}}}{{dr}}\right|\end{eqnarray} \tag{ 19 }$$

in spherical symmetry. (In this section we use the symbol "$d$" instead of "∂" because of time independence.) It is equal to the ratio of the photon mean free path to the "${E}_{\mathrm{rad}}$ scale height" (Paper I). Taking the derivative of Equation (17) yields

$$\begin{eqnarray}&&\displaystyle \frac{d\ \mathrm{ln}\ {T}_{\mathrm{rad}}^{4}}{d\ \mathrm{ln}\ r}=-2+\displaystyle \frac{d\ \mathrm{ln}\ L}{d\ \mathrm{ln}\ r}-\displaystyle \frac{d\ \mathrm{ln}\ {f}_{\mathrm{red}}}{d\ \mathrm{ln}\ r}.\end{eqnarray} \tag{ 20 }$$

The last term, by Equation (18), is

$$\begin{eqnarray}&&\displaystyle \frac{d\ \mathrm{ln}\ {f}_{\mathrm{red}}}{d\ \mathrm{ln}\ r}=\displaystyle \frac{d\ \mathrm{ln}\ R}{d\ \mathrm{ln}\ r}\left(1+\displaystyle \frac{d\mathrm{ln}\lambda }{d\ \mathrm{ln}\ R}\right).\end{eqnarray} \tag{ 21 }$$

These expressions are exact and general. We now proceed to expand the last equation.

By the definition of R, the first factor in Equation (21) is

$$\begin{eqnarray}&&\displaystyle \frac{d\ \mathrm{ln}\ R}{d\ \mathrm{ln}\ r}=\displaystyle \frac{d}{d\ \mathrm{ln}\ r}\left(\mathrm{ln}\left[\displaystyle \frac{d\ \mathrm{ln}\ {T}_{\mathrm{rad}}^{4}}{d\ \mathrm{ln}\ r}\right]-\mathrm{ln}\left[{\kappa }_{{\rm{R}}}\rho r\right]\right)\end{eqnarray} \tag{ 22 }$$

because ${E}_{\mathrm{rad}}={{aT}}_{\mathrm{rad}}^{4}$. Inserting recursively Equation (20) into Equation (22) is likely not fruitful as it generates derivatives of ${f}_{\mathrm{red}}$ of ever-higher order. However, it will prove instructive to perform this once. In full generality, this yields

$$\begin{eqnarray}&&\displaystyle \frac{d\ \mathrm{ln}\ R}{d\ \mathrm{ln}\ r}=\displaystyle \frac{{d}_{{\left(\mathrm{ln}r\right)}^{2}}^{2}\mathrm{ln}L-{d}_{{\left(\mathrm{ln}r\right)}^{2}}^{2}\mathrm{ln}{f}_{\mathrm{red}}}{-2+{d}_{\mathrm{ln}r}\mathrm{ln}L-{d}_{\mathrm{ln}r}\mathrm{ln}{f}_{\mathrm{red}}}-\displaystyle \frac{d\mathrm{ln}{\kappa }_{{\rm{R}}}}{d\ \mathrm{ln}\ r}+\displaystyle \frac{1}{2},\end{eqnarray} \tag{ 23 }$$

where ${d}_{{x}^{n}}^{n}f(x)\equiv {d}^{n}f(x)/{{dx}}^{n}$. For this derivation we have made use of the fact that $\rho \propto {r}^{-3/2}$ in the accretion flow.

The second factor in Equation (21) depends only on the choice of the flux limiter and can be computed easily independently of a simulation. If one takes the flux limiter used in Ensman (1994),

$$\begin{eqnarray}&&\lambda (R)=\displaystyle \frac{1}{3+R},\end{eqnarray} \tag{ 24 }$$

one obtains for the derivative term the simple expression

$$\begin{eqnarray}&&1+\displaystyle \frac{d\mathrm{ln}\lambda }{d\ \mathrm{ln}\ R}=3\lambda .\end{eqnarray} \tag{ 25 }$$

Note that this flux limiter is actually not physical (Levermore 1984), but it recovers the correct limits of free-streaming and diffusion (see Paper I). When the rational approximation of the Levermore & Pomraning (1981) flux limiter

$$\begin{eqnarray}&&\lambda (R)=\displaystyle \frac{2+R}{6+3R+{R}^{2}},\end{eqnarray} \tag{ 26 }$$

which we take by default in this work, is used instead, the result is

$$\begin{eqnarray}&&1+\displaystyle \frac{d\mathrm{ln}\lambda }{d\ \mathrm{ln}\ R}=1-\displaystyle \frac{{R}^{2}\left(R+4\right)}{(2+R)(6+3R+{R}^{2})}.\end{eqnarray} \tag{ 27 }$$

For $R\ll 2$ (the diffusion limit), this function (the right-hand side) is equal to 3λ. In the other limit of $R\gg 2$, it converges to simply λ; in any case, it remains roughly proportional to the flux limiter λ, a fact we shall use in the discussion below.

Combining Equations (20), (21), (23), and (25), we obtain

$$\begin{eqnarray}&&\displaystyle \frac{d\ \mathrm{ln}\ T}{d\ \mathrm{ln}\ r}\approx -\displaystyle \frac{1}{2}+\displaystyle \frac{3\lambda (R)}{4}\left[\displaystyle \frac{{d}_{{\left(\mathrm{ln}r\right)}^{2}}^{2}\mathrm{ln}{f}_{\mathrm{red}}}{2+{d}_{\mathrm{ln}r}\mathrm{ln}{f}_{\mathrm{red}}}-\displaystyle \frac{d\mathrm{ln}{\kappa }_{{\rm{R}}}}{d\ \mathrm{ln}\ r}+\displaystyle \frac{1}{2}\right]\end{eqnarray} \tag{ 28 }$$

in the limit of a radially constant luminosity (in the sense that $\left|d\ \mathrm{ln}\ L/d\ \mathrm{ln}\ r\right|\ll 2$, which is the case here as shown below). The temperature here was written as T because we find that ${T}_{\mathrm{gas}}={T}_{\mathrm{rad}};$ if the gas and radiation are not coupled, T should be taken to refer to ${T}_{\mathrm{rad}}$ only. We used Equation (25) because it makes it explicit that the second term on the right-hand side of Equation (28) is (in general roughly) proportional to the flux limiter $\lambda (R)$. While at least in this form, Equation (28) is not predictive because ${f}_{\mathrm{red}}(r)$ is needed for its derivatives, it does exhibit the link between the temperature and the opacity slopes, as we now discuss.

Looking in detail at Equation (28), we see that the first term on the right-hand side corresponds to the result that $T\propto {r}^{-1/2}$ for a radially constant luminosity. It seems likely that this term will dominate in the free-streaming (often termed "optically thin") regime because then λ (and thus the second term) goes to zero.

For the temperature slope to differ significantly from −1/2, it is sufficient for either the opacity slope or the term with derivatives of ${f}_{\mathrm{red}}$ to be large, and necessary for the radiation to be sufficiently diffusive (λ not too small).

The second term of Equation (28) is interesting: it relates the local slope of the opacity to that of the temperature. In the free-streaming regime, the flux limiter λ goes to zero and thus also the second term in Equation (28) because it is multiplied by λ). Therefore, even strong variations of ${\kappa }_{{\rm{R}}}$ with radius will only minorly affect the temperature structure; this is as expected for the limit of infinite optical mean free path, in which the radiation does not interact with the opacity carrier. In the other limiting case of diffusion, $\lambda \approx 1/3$ and the opacity variations are important. The evaporation of dust as the material moves in leads to a large $d\mathrm{ln}{\kappa }_{{\rm{R}}}/d\ \mathrm{ln}\ r$, and if λ is not too small, this will be able to slow the decrease of the temperature outward.

Moreover, it does not seem necessary that the radiation be free streaming (${f}_{\mathrm{red}}\to 1$) in order to have T ∝ r^−1/2, which is usually associated with free streaming. If ${\kappa }_{{\rm{R}}}$ is constant with radius and (as perhaps a consequence) ${f}_{\mathrm{red}}$ is sufficiently constant, and λ has an intermediate value (e.g., λ ∼ 0.1), the −1/2 term in Equation (28) can dominate.

We show in the top panel of Figure 6 the luminosity, temperature, reduced flux, and opacity on logarithmic scales for the nominal case shown in Figure 9. Indeed, the luminosity is effectively radially constant.

The middle panel of Figure 6 shows the terms on the right-hand side of Equation (28) and parts thereof. In the free-streaming region at r = 2–8 R_J, the first and second derivatives of ${f}_{\mathrm{red}}$ are nearly zero, and the strong values of the opacity slope do not bring T away from an r^−1/2 scaling because the radiation is nearly free-streaming: ${f}_{\mathrm{red}}\approx 1$ (see top panel), so that λ ≈ 0. In the flat temperature part at r = 8–16 R_J (shaded region), the opacity slope term dominates, but the term with the derivatives of ${f}_{\mathrm{red}}$ is similar in magnitude, with their signed sum dominating the −1/2 term. The result, along with a lower ${f}_{\mathrm{red}}$ (i.e., higher λ), is a different (namely, almost zero) temperature slope. At $r\gt 16\,{R}_{{\rm{J}}}$, the radiation remains somewhat diffusive with ${f}_{\mathrm{red}}\sim 0.3$, but the opacity is nearly constant. The temperature slope is therefore set by the leading −1/2 term as well as the two other terms in the brackets and the factor of ≈3λ (depending on the flux limiter model). The net result is that the −1/2 dominates.

The slope from Equation (28) is shown in the bottom panel of Figure 6 and compared to the actual value. The agreement is excellent, which mainly confirms the approximation ${dL}/{dr}\approx 0$ because the equation is otherwise exact for a freefall density profile. The choice of the flux limiter barely makes a difference.

Note finally that Equations (19) and (20) imply that when the ratio L1/${f}_{\mathrm{red}}$ is radially constant, the radiation quantity R is given by

$$\begin{eqnarray}&&R=\displaystyle \frac{2}{{\kappa }_{{\rm{R}}}\rho r}\,\,\,\,\,\,({\rm{c}}{\rm{o}}{\rm{n}}{\rm{s}}{\rm{t}}{\rm{a}}{\rm{n}}{\rm{t}}\,L/{f}_{{\rm{r}}{\rm{e}}{\rm{d}}}).\end{eqnarray} \tag{ 29 }$$

This had been derived in Paper I, above their Equation (25).

This analysis provides an approximate analytical understanding of the link between the opacity and the temperature. What is less clear from this derivation is how much the geometry accurately reflects the realistic transport of radiation, even within the gray approximation, and to what extent it depends on the approximations (e.g., concerning the angular distribution of the specific intensity) inherent to moments-based methods such as FLD. Nevertheless, it might prove insightful to attempt a similar derivation using M1 radiation transport (e.g., González et al. 2007; Hanawa & Audit 2014).

Figure 6 suggests that there is a relationship between the reduced flux ${f}_{\mathrm{red}}$ and the Rosseland mean opacity ${\kappa }_{{\rm{R}}}$. Defining the logarithmic slope

$$\begin{eqnarray}&&\beta \equiv \displaystyle \frac{d\ \mathrm{ln}\ (L/{f}_{\mathrm{red}})}{d\ \mathrm{ln}\ r},\end{eqnarray} \tag{ 30 }$$

this relationship can be obtained by combining Equations (17)–(19) to yield

$$\begin{eqnarray}&&\,\displaystyle \frac{1}{{f}_{\mathrm{red}}}=1+\displaystyle \frac{3}{\left|2-\beta \right|}{\kappa }_{{\rm{R}}}\rho r\,\end{eqnarray} \tag{ 31a }$$

$$\begin{eqnarray}&&\,=\,1+\displaystyle \frac{3}{2}{\kappa }_{{\rm{R}}}\rho r\,\,\,\,\,\,({\rm{c}}{\rm{o}}{\rm{n}}{\rm{s}}{\rm{t}}{\rm{a}}{\rm{n}}{\rm{t}}L/{f}_{{\rm{r}}{\rm{e}}{\rm{d}}}),\end{eqnarray} \tag{ 31b }$$

which holds anywhere in the flow. This generalizes Equation (25) in Paper I. We used λ from Equation (24), but it is trivial to repeat the derivation for another flux limiter.

Zoom In Zoom Out Reset image size

Equation (31) explains why ${f}_{\mathrm{red}}$ drops in the dust destruction region in Figure 6 (gray band there): the opacity increases such that ${\kappa }_{{\rm{R}}}\rho r\gtrsim 1$, thereby leading to a drop in ${f}_{\mathrm{red}}$. Physically, this has the intuitive explanation that the radiation becomes less freely streaming and starts to diffuse. In the dust destruction region, $\beta \approx 1.87$ (and relatively constant) due to the outward decreasing ${f}_{\mathrm{red}}$. From Equation 31(a) this implies that $1/{f}_{\mathrm{red}}$ will be larger, i.e., ${f}_{\mathrm{red}}$ smaller, than predicted by Equation 31(b), but the trend is the same. Quantitiatively, this works well in this example, but to obtain the exact value of ${f}_{\mathrm{red}}$, one would have to use the Levermore & Pomraning (1981) flux limiter because this was chosen for the simulations. Thus, Equation (31) explains the sudden drop of ${f}_{\mathrm{red}}$ where the opacity suddenly increases to become important in the sense of ${\kappa }_{{\rm{R}}}\rho r\gtrsim 1$.

5.3.1. Calculation Setup

We now turn to explaining the shock temperatures found in Figure 3, in particular the points that deviate from Equation 6(b). We recall that this temperature more precisely refers to the radiation temperature ${T}_{\mathrm{rad}}$, should it ever be found to differ from ${T}_{\mathrm{gas}}$, which is not the case here, however.

Looking at the results of Section 4.1 again, there is a tight anticorrelation (not shown) between ${f}_{\mathrm{red}}$ and ${\kappa }_{{\rm{R}}}\rho r$ both evaluated immediately upstream of the shock. Empirically across our grid of models, $\beta ({{r}_{\mathrm{shock}}}^{+})\sim -0.1$ typically, which in absolute value is $\ll 2$. Therefore, the tight relationship comes from Equation 31(b). Thus while in general, $\beta ({{r}_{\mathrm{shock}}}^{+})$ is not known but could perhaps be estimated, in the following analysis we restrict ourselves to $\beta ({{r}_{\mathrm{shock}}}^{+})\ll 2$ and consequently use $R({{r}_{\mathrm{shock}}}^{+})=2/{\kappa }_{{\rm{R}}}{\rho }_{\mathrm{ff}}{R}_{{\rm{p}}}$ at the shock (Equation (29)).

With this, Equation 6(a) can be rewritten as

$$\begin{eqnarray}&&\sigma {T}_{\mathrm{shock}}^{4}=\displaystyle \frac{{L}_{\mathrm{dnstr}}+{\eta }^{\mathrm{kin}}({T}_{\mathrm{shock}})\tfrac{{{GM}}_{{\rm{p}}}\dot{M}}{{R}_{{\rm{p}}}}}{{{f}_{\mathrm{red}}}^{+}({T}_{\mathrm{shock}})16\pi {{R}_{{\rm{p}}}}^{2}}\,\end{eqnarray} \tag{ 32a }$$

$$\begin{eqnarray}&&\,\approx \,\left(1+\displaystyle \frac{3}{2}{\kappa }_{{\rm{R}}}({T}_{\mathrm{shock}}){\rho }_{\mathrm{ff}}{R}_{{\rm{p}}}\right)\displaystyle \frac{{L}_{\mathrm{dnstr}}+\tfrac{{{GM}}_{{\rm{p}}}\dot{M}}{{R}_{{\rm{p}}}}}{16\pi {{R}_{{\rm{p}}}}^{2}},\end{eqnarray} \tag{ 32b }$$

which is an implicit equation for ${T}_{\mathrm{shock}}$. The second line is valid when the resulting ${\mathscr{M}}\gtrsim 2.5$ (so that ${\eta }^{\mathrm{kin}}\approx 1$), which should be verified a posteriori, and was written for the Ensman (1994) flux limiter through Equation (31). The downstream luminosity ${L}_{\mathrm{dnstr}}$ is the sum of the luminosity coming from the deep interior and of the compression luminosity: ${L}_{\mathrm{dnstr}}={L}_{\mathrm{int}}\,+{L}_{\mathrm{compr}}$. While ${L}_{\mathrm{compr}}$ might depend on ${T}_{\mathrm{shock}}$, we effectively absorb this dependency into ${L}_{\mathrm{int}}$, treated as a free parameter. Still, the generalization ${L}_{\mathrm{dnstr}}\to {L}_{\mathrm{dnstr}}({T}_{\mathrm{shock}})$ could be readily made. With these assumptions, in Equation 32(b) the shock temperature enters on the right-hand side only through the opacity.

In the 100% efficiency limit, the shock temperature that solves Equation 32(b) thus has the limiting cases

$$\begin{eqnarray}&&{T}_{\mathrm{sh},\mathrm{fs}}\equiv \,{T}_{\mathrm{shock}}({\kappa }_{{\rm{R}}}\rho {R}_{{\rm{p}}}\ll 1)\,\end{eqnarray} \tag{ 33a }$$

$$\begin{eqnarray}&&=\,{{\ell }}^{1/4}{\left(\displaystyle \frac{G}{16\pi \sigma }\right)}^{1/4}\displaystyle \frac{{M}_{{\rm{p}}}^{1/4}{\dot{M}}^{1/4}}{{{R}_{{\rm{p}}}}^{3/4}}\,\end{eqnarray} \tag{ 33b }$$

$$\begin{eqnarray}\begin{array}{rcl} & = & 2315\,{\rm{K}}\,{{\ell }}^{1/4}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{2{R}_{{\rm{J}}}}\right)}^{-3/4}\\ & & \times {\left(\displaystyle \frac{\dot{M}}{{10}^{-2}{M}_{\oplus }{\mathrm{yr}}^{-1}}\right)}^{1/4}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{1{M}_{{\rm{J}}}}\right)}^{1/4},\end{array}\end{eqnarray} \tag{ 33c }$$

$$\begin{eqnarray}&&{T}_{\mathrm{sh},\mathrm{diff}}\equiv \,{T}_{\mathrm{shock}}({\kappa }_{{\rm{R}}}\rho {R}_{{\rm{p}}}\gg 1)\,\end{eqnarray} \tag{ 34a }$$

$$\begin{eqnarray}&&\,=\,{{\ell }}^{1/4}{\left(\displaystyle \frac{3\sqrt{G}}{128\sqrt{2}{\pi }^{2}\sigma }\right)}^{1/4}\displaystyle \frac{{\dot{M}}^{1/2}{\kappa }_{{\rm{R}}}^{1/4}{M}_{{\rm{p}}}^{1/8}}{{{R}_{{\rm{p}}}}^{7/8}}\,\end{eqnarray} \tag{ 34b }$$

$$\begin{eqnarray}\begin{array}{lll} & \,= & 3220\,{\rm{K}}\,{{\ell }}^{1/4}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{2{R}_{{\rm{J}}}}\right)}^{-7/8}{\left(\displaystyle \frac{{\kappa }_{{\rm{R}}}}{1{\mathrm{cm}}^{2}{{\rm{g}}}^{-1}}\right)}^{1/4}\\ & & \times {\left(\displaystyle \frac{\dot{M}}{{10}^{-2}{M}_{\oplus }{\mathrm{yr}}^{-1}}\right)}^{1/2}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{1{M}_{{\rm{J}}}}\right)}^{1/8},\end{array}\end{eqnarray} \tag{ 34c }$$

with ${\ell }=1+{L}_{\mathrm{dnstr}}/{L}_{\mathrm{acc},\max }$, without needing to assume that ${L}_{\mathrm{dnstr}}$/${L}_{\mathrm{acc},\max }$ is small. We defined ${T}_{\mathrm{sh},\mathrm{fs}}$ as the "free-streaming shock temperature" given by Equation 6(b), but now generalized to include a downstream luminosity ${L}_{\mathrm{dnstr}}$ (we consider ${L}_{\mathrm{dnstr}}=0$ in this work). Similarly, ${T}_{\mathrm{sh},\mathrm{diff}}$ is the "diffusive shock temperature," which obtains when the preshock material is diffusive. The maximum accretion luminosity is ${L}_{\mathrm{acc},\max }\,={{GM}}_{{\rm{p}}}\dot{M}/{R}_{{\rm{p}}}$, with the actual luminosity at the shock ${L}_{\mathrm{acc}}\,={\eta }^{\mathrm{kin}}{L}_{\mathrm{acc},\max }$ (Paper I). The opacity ${\kappa }_{{\rm{R}}}$ appearing in the expressions for ${T}_{\mathrm{sh},\mathrm{diff}}$ is either the constant opacity or less trivially, the value satisfying Equation 32(b), as detailed below. Equation 34(c) should replace Equation 28(b) of Paper I because the prefactor there was inexact and the exponent of the ${R}_{{\rm{p}}}$ factor was missing a minus sign. Equations 34(b) and (c) were written to leading order in the quantity ${\kappa }_{{\rm{R}}}{\rho }_{\mathrm{ff}}{r}_{\mathrm{shock}}$, in which case they are also independent of the flux limiter. Note that if ${L}_{\mathrm{dnstr}}$ were to dominate ${L}_{\mathrm{acc},\max }$ (${\ell }\gg 1$), the functional dependence of the shock temperature on the various parameters would be different from what Equations (33) and (34) naively suggest.

As an example, for the $(\gamma =1.1,\kappa =1$ cm² g ${}^{-1})$ simulation in Figure 2 of Paper I, Equation (34) predicts ${T}_{\mathrm{sh},\mathrm{fs}}\approx 3600$ K, which is close to the actual shock temperature 3500 K. This is satisfying, especially given that there ${\kappa }_{{\rm{R}}}\rho {R}_{{\rm{p}}}=2.3$, which is not entirely $\gg 1$, and that the Levermore & Pomraning (1981) flux limiter had been used and not the Ensman (1994) one as for Equations 34(b) and (c), which makes a difference in this transition regime between free streaming and diffusion.

In general, we solve Equation 32(b) by writing it as

$$\begin{eqnarray}{\left(\displaystyle \frac{{T}_{\mathrm{shock}}}{{T}_{\mathrm{sh},\mathrm{fs}}}\right)}^{4}=\left\{\begin{array}{ll}1+\tfrac{3}{2}x & ({\rm{E}}94)\\ \left(1+\tfrac{3}{2}x+\tfrac{3}{2}{x}^{2}\right)/\left(1+x\right) & (\mathrm{LP}81)\end{array},\right.\end{eqnarray} \tag{ 35 }$$

with

$$\begin{eqnarray}&&x({T}_{\mathrm{shock}})\equiv {\kappa }_{{\rm{R}}}({T}_{\mathrm{shock}}){\rho }_{\mathrm{ff}}{R}_{{\rm{p}}}.\end{eqnarray} \tag{ 36 }$$

The first and second versions of the right-hand side of Equation (35) hold for the flux limiters of Ensman (1994, E94) and Levermore & Pomraning (1981, LP81), respectively. The equation is solved by root finding, simply starting slightly below $T={T}_{\mathrm{sh},\mathrm{fs}}$ and stepping up in temperature to find an intersection of the two sides of Equation (35).

5.3.2. Semianalytical Shock Temperature Solutions

In Figure 7 we display the shock temperature as a function of planet radius according to Equation (35). We consider ${M}_{{\rm{p}}}=1.3$ and 5 M_J, take $\dot{M}={10}^{-4,-3,-2}\,{M}_{\oplus }\,{{\rm{yr}}}^{-1}$, and vary ${R}_{{\rm{p}}}$ from 1 to about 7 R_J, but truncating at lower ${R}_{{\rm{p}}}$ for the lowest accretion rate. The downstream luminosity is taken to be zero or smaller than but not entirely negligible compared to the accretion luminosity, ${\mathrm{log}}_{10}\left({L}_{\mathrm{dnstr}}/{L}_{\odot }\right)={\mathrm{log}}_{10}\left(\dot{M}/{M}_{\oplus }\,{{\rm{yr}}}^{-1}\right)-1$ for definiteness. For reference, the corresponding preshock (freefall) velocity is indicated on the top axis. The solution is also compared to the free-streaming shock temperature ${T}_{\mathrm{sh},\mathrm{fs}}$ (Equation (33)).

Zoom In Zoom Out Reset image size

Figure 7 reveals that opacity effects can alter the shock temperature only for a relatively small range of parameters. Specifically, for $\dot{M}\sim {10}^{-3}$–${10}^{-2}\,{M}_{\oplus }\,{{\rm{yr}}}^{-1}$, planets of low mass (${M}_{{\rm{p}}}\lesssim 3\,{M}_{{\rm{J}}}$) with moderate to large radii (${R}_{{\rm{p}}}\approx 3$–10 R_J) could have their ${T}_{\mathrm{shock}}$ increased by up to tens of percent, by a few 100 K. This might be important especially for a nonzero ${L}_{\mathrm{int}}$, leading to a higher shock temperature than naively expected.

5.3.3. Regimes of the Shock Temperature

To understand these results graphically, we show in the bottom row of Figure 7 the fourth root of the left- and right-hand sides of Equation (35). We see that the points that deviate in Figure 3 can do so for two reasons. First of all, a high constant (temperature- and density-independent) opacity will increase the shock temperature by reducing the effective speed of light ahead of the shock. The threshold for this is ${\kappa }_{{\rm{R}}}{\rho }_{\mathrm{ff}}{R}_{{\rm{p}}}\sim 1$, as mentioned above, and an example is for $(\dot{M}={10}^{-3}\,{M}_{\oplus }\,{{\rm{yr}}}^{-1},{M}_{{\rm{p}}}=1.3\,{M}_{{\rm{J}}},{R}_{{\rm{p}}}=3\,{R}_{{\rm{J}}})$, labeled "A," where the opacity of the most refractory dust component leads to a higher shock temperature than ${T}_{\mathrm{sh},\mathrm{fs}}\,=1000$ K: ${T}_{\mathrm{shock}}\approx 1080$ K for the nominal dust model and ${T}_{\mathrm{shock}}\approx 1200$ K for the curve sticking out most in Figure 2. Second of all, the case of nonconstant high (${\kappa }_{{\rm{R}}}{\rho }_{\mathrm{ff}}{R}_{{\rm{p}}}\gtrsim 1$) opacities opens up the possibility of multiple solutions.

If ${\kappa }_{{\rm{R}}}\rho {R}_{{\rm{p}}}\ll 1$ at the free-streaming shock temperature ${T}_{\mathrm{sh},\mathrm{fs}}$, it is necessary and sufficient for the opacity slope

$$\begin{eqnarray}&&{\alpha }_{\kappa }\equiv {\left(\displaystyle \frac{\partial \mathrm{ln}{\kappa }_{{\rm{R}}}}{\partial \mathrm{ln}T}\right)}_{\rho }\end{eqnarray} \tag{ 37 }$$

to satisfy ${\alpha }_{\kappa }\gt {\alpha }_{\kappa }^{\mathrm{crit}}=1$ at higher temperatures and for a sufficient T range in order to have one high-T solution. A drop of ${\alpha }_{\kappa }$ below ${\alpha }_{\kappa }^{\mathrm{crit}}$ at a higher temperature will lead to a third solution. The stronger the slopes, the closer these higher-T solutions will be to ${T}_{\mathrm{sh},\mathrm{fs}}$.

If on the other hand ${\kappa }_{{\rm{R}}}\rho {R}_{{\rm{p}}}\gtrsim 1$ at ${T}_{\mathrm{sh},\mathrm{fs}}$ (i.e., the preshock gas would be diffusive already if it were at that temperature), ${T}_{\mathrm{sh},\mathrm{fs}}$ will not be a solution. Provided the opacity slope ${\alpha }_{\kappa }\lt 1$ at some higher temperature, there will only be a single¹¹ solution at high temperature. Again, the stronger the slopes $| {\alpha }_{\kappa }|$ in the increasing and decreasing parts, the closer the solution will be to ${T}_{\mathrm{sh},\mathrm{fs}}$.

If ${\kappa }_{{\rm{R}}}\rho {R}_{{\rm{p}}}\gtrsim 1$ at ${T}_{\mathrm{sh},\mathrm{fs}}$ and ${\alpha }_{\kappa }$ remained $\gt 1$ at higher temperatures, however, there would formally be no solution: try it as it may by increasing its temperature, the shock would not be able to radiate away the kinetic energy. It is not clear, however, whether in this case we would still find a Mach number such that ${\eta }^{\mathrm{kin}}=1$, or if the model otherwise breaks down. Fortunately, for gas mixtures as considered here, ${\alpha }_{\kappa }$ does drop again, so that the situation does not arise.

Qualitatively, these results do not depend much on different model settings. Including a low downstream luminosity (coming from compression, an interior luminosity, or both) does not significantly change the shock temperature(s) of Figure 7. Changing the opacity (the dust model or, at higher temperature, the abundance of water; Malygin et al. 2014) or the flux limiter also has a clear but limited effect.

5.3.4. Mach Number

Within the restriction of negligible ${L}_{\mathrm{dnstr}}$, the Mach number ${\mathscr{M}}$ increases with increasing ${M}_{{\rm{p}}}$ or $1/{R}_{{\rm{p}}}$, but, perhaps counterintuitively, decreases with increasing accretion rate. Using Equation (6), the Mach number is

$$\begin{eqnarray}&&{\mathscr{M}}=2{\left(\displaystyle \frac{{G}^{3}\pi \sigma {\rm{\Delta }}{f}_{{\rm{r}}{\rm{e}}{\rm{d}}}}{{\eta }^{{\rm{k}}{\rm{i}}{\rm{n}}}}\right)}^{1/8}\sqrt{\displaystyle \frac{\mu {m}_{{\rm{H}}}}{\gamma {k}_{{\rm{B}}}}}\,{\left(\displaystyle \frac{{M}_{{\rm{p}}}^{3}}{\dot{M}{R}_{{\rm{p}}}}\right)}^{1/8}\,\end{eqnarray} \tag{ 38a }$$

$$\begin{eqnarray}\begin{array}{r}\begin{array}{c}\begin{array}{ccc} & = & 30{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{1.3{M}_{{\rm{J}}}}\right)}^{3/8}\,{\left(\displaystyle \frac{\mu }{2.29}\right)}^{1/2}\,{\left(\displaystyle \frac{\gamma }{1.44}\right)}^{-1/2}\\ & & \times \,{\left(\displaystyle \frac{\dot{M}}{{10}^{-2}{M}_{\oplus }{{\rm{y}}{\rm{r}}}^{-1}}\right)}^{-1/8}\,{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-1/8},\end{array}\end{array}\end{array}\end{eqnarray} \tag{ 38b }$$

where we took ${\eta }^{\mathrm{kin}}=1$ and ${\rm{\Delta }}{f}_{\mathrm{red}}=1$ for the second expression. Thus the Mach number depends only moderately on the mass (${\mathscr{M}}\propto {M}_{{\rm{p}}}^{3/8\approx 0.4}$) and accretion rate (given that the latter ranges over several orders of magnitude) and barely on the radius.

We show in Figure 8 the preshock Mach number obtained from Equation 32(b) with μ = 1.23, and also plot Equation (38) for reference. Mach numbers increase with decreasing accretion rate, reaching ${\mathscr{M}}\sim 30$ (in particular at higher masses, not shown) for $\dot{M}={10}^{-4}\,{M}_{\oplus }\,{\mathrm{yr}}^{-1}$. The Mach number at ${T}_{\mathrm{sh},\mathrm{fs}}$, given by Equation (38), provides an upper bound but ${\mathscr{M}}$ still remains securely above ${\mathscr{M}}\approx 2.5$. There, ${\eta }^{\mathrm{kin}}$ is near 100%, which justifies our approximation a posteriori. Also, taking ${L}_{\mathrm{dnstr}}$ not large (see Figure 8) but nonzero barely changes these results.

Zoom In Zoom Out Reset image size

5.3.5. Discussion of the Analytical Shock Temperature

While the use of an ideal but nonperfect EOS is deferred to a later article, we study in this section the importance of the constant μ and γ. For conciseness, we refer to a particular (μ, γ) combination as an EOS.

Figure 9 shows the profiles for the four extreme combinations (μ, γ), where μ = 1.23 (atomic) or μ = 2.29 (molecular), and γ = 1.1 (occurring for dissociation or ionization) or γ = 5/3 (monatomic, no dissociation or ionization) together with the nominal case (μ = 2.29, γ = 1.44), i.e., for a perfect-gas H₂–He mixture. The density in the accretion flow being set by mass conservation, it is independent of the EOS, but the jump across the shock scales as ${\rm{\Delta }}\rho \propto \mu$. The pressure in the accretion flow does scale as P ∝ 1/μ, but being a dynamical quantity, the ram pressure (i.e., the postshock pressure) is the same in all cases. While this might be counterintuitive, the shock temperature (see inset) is also independent of (μ, γ), which was to be expected from the analytical estimates of the shock temperature because the microphysical parameters do not enter anywhere. The opacities ${\kappa }_{{\rm{R}}}$ and ${\kappa }_{{\rm{P}}}$ are the same here because we use the same tables in all cases (computed using an nonperfect EOS) while the opacity is fundamentally a function of temperature and density (but not mean molecular weight).

Zoom In Zoom Out Reset image size

6.1.1. Entropy

One of the quantities that change with the EOS is the entropy calculated with the chosen (μ, γ). Note that it differs from the entropy that would be obtained with nonperfect EOS shock simulations. This (μ, γ)-dependent entropy is given in general (but for constant (μ, γ), i.e., outside of chemical reactions such as conversion of ortho- to parahydrogen, dissociation, or ionization) by

$$\begin{eqnarray}&&s={s}_{0}+\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{\mathrm{ln}(10)}{\mu }{\mathrm{log}}_{10}\displaystyle \frac{T}{{T}_{0}}-\displaystyle \frac{\mathrm{ln}(10)}{\mu }{\mathrm{log}}_{10}\displaystyle \frac{P}{{P}_{0}},\end{eqnarray} \tag{ 39 }$$

where s₀ is a constant, and T₀ and P₀ are an arbitrary reference temperature and pressure, respectively. The form is similar as a function of (P, ρ) (e.g., Rafikov 2016).

The radial profile of the entropy given by Equation (39) is shown in Figure 9 for various EOS. We plot in fact the entropy difference with respect to the preshock location because we are interested in the radial profiles of s for the individual (μ, γ) cases. It increases outward for γ = 5/3 because it is greater than the critical value of 4/3 (Section 3.3.2 of Paper I) but decreases for γ = 1.1 (because then γ < 4/3) outward in the accretion flow. Convection is not expected there because of the supersonic motion, however.

6.1.2. Luminosity

Another difference between the various EOS is in the luminosity profile, as Figure 9 shows. We now proceed to explain it. The luminosity L varies throughout the accretion flow by an amount that depends on the EOS, not directly on the optical depth. This perhaps surprising statement can be derived from considering the steady-state (time-independent) version of the evolution equation for the total energy (so that the ±Λ terms cancel even in nonequilibrium radiation transport; Equation (3)) and writing the gravity field term as a potential:

$$\begin{eqnarray}&&\displaystyle \frac{dL}{dr}=| \dot{M}| \,\left(\displaystyle \frac{dh}{dr}+\displaystyle \frac{d}{dr}\left[\displaystyle \frac{1}{2}{v}^{2}-\displaystyle \frac{{GM}_{{\rm{p}}}}{r}\right]\right)\end{eqnarray} \tag{ 40a }$$

$$\begin{eqnarray}&&=\,| \dot{M}| \,\left(\displaystyle \frac{dh}{dr}-\displaystyle \frac{1}{\rho }\displaystyle \frac{dP}{dr}\right)\,\end{eqnarray} \tag{ 40b }$$

$$\begin{eqnarray}&&=\,| \dot{M}| \,\left(\displaystyle \frac{{de}_{{\rm{i}}{\rm{n}}{\rm{t}}}}{dr}+P\displaystyle \frac{d(1/\rho )}{dr}\right)\,\end{eqnarray} \tag{ 40c }$$

$$\begin{eqnarray}&&=\,| \dot{M}| T\displaystyle \frac{{ds}}{{dr}},\,\end{eqnarray} \tag{ 40d }$$

where $h=H/\rho ={e}_{\mathrm{int}}+P/\rho$ is the specific enthalpy per mass and is given by $h=\gamma /(\gamma -1){k}_{{\rm{B}}}T/(\mu {m}_{{\rm{H}}})$ for a perfect EOS. Only Equation 40(a) was derived in Paper I. It did not require identifying the gas and radiation temperatures with each other, leaving the derivation valid also for 2-T radiation transport. Equation 40(b) is obtained from $\dot{M}=4\pi {r}^{2}\rho v$ (Equation 1(a)) and the time-independent momentum equation (Equation 1(b)). Equation 40(d) follows from the first law of thermodynamics,

$$\begin{eqnarray}&&{Tds}={{de}}_{\mathrm{int}}+{Pd}\left(\displaystyle \frac{1}{\rho }\right),\end{eqnarray} \tag{ 41 }$$

with the entropy s here not normalized by ${k}_{{\rm{B}}}\,{{m}_{{\rm{H}}}}^{-1}$. Note that Equation (40) both above and below the shock. In the pre-shock region, it is the small deviation of the velocity from a freefall profile that makes ${dL}/{dr}\propto {Tds}/{dr}$ and not $\propto \,dh/dr$.

In Paper I, we had argued that the radial non-constancy (over ∼ Hill-sphere scales) of the luminosity can be understood from the enthalpy profile when the second term in Equation 40(a) is negligible. While this is formally true, Equation (40) provides a more complete derivation that reveals that the relevant quantity is the entropy, without requiring restrictions on any term. Thus the radial luminosity profile is set not directly by the optical depth. Nevertheless, the analysis in Section 5.1 has shown that there is a link between the local opacity and its slope on the one hand and the temperature (and thus the entropy) profile on the other.

In the case of a perfect EOS, Equation 40(d) becomes

$$\begin{eqnarray}&&\displaystyle \frac{{dL}}{{dr}}=\displaystyle \frac{| \dot{M}| }{\gamma -1}\displaystyle \frac{{k}_{{\rm{B}}}T}{\mu {m}_{{\rm{H}}}}\left[\displaystyle \frac{d\ \mathrm{ln}\ T}{{dr}}-(\gamma -1)\displaystyle \frac{d\mathrm{ln}\rho }{{dr}}\right]\end{eqnarray} \tag{ 42a }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{{k}_{{\rm{B}}}| \dot{M}| }{\mu {m}_{{\rm{H}}}}\displaystyle \frac{T}{r}\left[\displaystyle \frac{1}{\gamma -1}\displaystyle \frac{d\ \mathrm{ln}\ T}{d\mathrm{ln}r}+\displaystyle \frac{3}{2}\right],\,\end{eqnarray} \tag{ 42b }$$

where the second line holds for a freefall density profile. Because the factor T/r in Equation 42(b) decreases outward, it is the layers closest to the shock that contribute most to the change in L between ${r}_{\mathrm{shock}}$ and ${r}_{\max }$, at least for a constant logarithmic temperature slope. Equation (42) reveals that, all other things being equal, the change in luminosity from the shock to the Hill sphere will be more important as γ tends to 1 (i.e., more isothermal), or for smaller mean molecular weight μ. The choice of γ can change the sign of ${dL}/{dr}$, whereas μ cannot.

In fact, an estimate for the change in luminosity between the shock and a given distance in the flow, in particular the accretion radius, can be derived by using the constant-$(L/{f}_{\mathrm{red}})$ temperature profile (see Equation (20))

$$\begin{eqnarray}&&T(r)={T}_{\mathrm{shock}}{\left(\displaystyle \frac{{r}_{\mathrm{shock}}}{r}\right)}^{-1/2},\end{eqnarray} \tag{ 43 }$$

where ${T}_{\mathrm{shock}}$ is the shock temperature, which by Equation (35) can be higher than the free-streaming temperature. (Note that $T\propto {r}^{-1/2}$ can hold not only for free-streaming.) In Figure 6(c) or Figure 9, the T(r) profiles are rarely steeper than Equation (43), so that it will provide an upper bound. Inserting Equation (43) in Equation (42) yields

$$\begin{eqnarray}&&\displaystyle \frac{{{dL}}_{(43)}}{{dr}}=-\displaystyle \frac{| \dot{M}| }{\gamma -1}\displaystyle \frac{{k}_{{\rm{B}}}{T}_{\mathrm{shock}}}{\mu {m}_{{\rm{H}}}}\displaystyle \frac{4-3\gamma }{2{R}_{{\rm{p}}}}{\left(\displaystyle \frac{r}{{R}_{{\rm{p}}}}\right)}^{3/2},\end{eqnarray} \tag{ 44 }$$

using the subscript "(43)" on L to recall that a temperature profile given by Equation (43) was assumed. Thus, the change (drop or increase) in ${L}_{(43)}$ between ${R}_{{\rm{p}}}$ and ${R}_{\mathrm{acc}}$, neglecting terms of order $\sqrt{{R}_{{\rm{p}}}/{R}_{\mathrm{acc}}}\ll 1$, is given by

$$\begin{eqnarray}&&{\rm{\Delta }}{L}_{(43)}={\int }_{{R}_{{\rm{p}}}}^{{R}_{\mathrm{acc}}}\displaystyle \frac{{dL}}{{dr}}=\displaystyle \frac{3\gamma -4}{\gamma -1}\displaystyle \frac{{k}_{{\rm{B}}}{T}_{\mathrm{shock}}| \dot{M}| }{\mu {m}_{{\rm{H}}}},\end{eqnarray} \tag{ 45 }$$

with the luminosity reaching the local circumstellar disk (the nebula) equal to

$$\begin{eqnarray}&&L({R}_{\mathrm{acc}})\approx {L}_{\mathrm{int}}+{L}_{\mathrm{acc}}+{\rm{\Delta }}{L}_{(43)}.\end{eqnarray} \tag{ 46 }$$

Equation (45) is more general than Equation (32) in Paper I. In the usual limit ${R}_{\mathrm{acc}}\gg {R}_{{\rm{p}}}$, the luminosity change ${\rm{\Delta }}{L}_{(43)}$ is thus independent of ${R}_{\mathrm{acc}}$, and thus also of ${k}_{\mathrm{Lissauer}}=1/3$ (as we assume here) for the accretion radius ${R}_{\mathrm{acc}}\approx {k}_{\mathrm{Lissauer}}{R}_{\mathrm{Hill}}$. For γ > 4/3, the luminosity increases outward (${\rm{\Delta }}{L}_{(43)}\gt 0$), whereas for γ < 4/3 there is a decrease in luminosity (${\rm{\Delta }}{L}_{(43)}\lt 0$).

Equation (45) is an upper bound in the case γ < 4/3 because any flattening of the luminosity profile as in Figure 9 will slow down the decrease of T(r). For γ > 4/3 the estimate is more like a lower bound, but because the critical γ = 4/3 ≈ 1.33 is not far from γ = 1.44 for molecular hydrogen with helium, it is also roughly equal to the actual change in that case. This can be seen graphically in Figure 9.

Thus, in Figure 9, L(r) increases outward for the high-γ cases (γ = 5/3) as well as for the nominal case (γ = 1.44), while it slightly decreases overall for the low-γ cases (γ = 1.1). The mean molecular weight also plays a role, and, in all, L at roughly the Hill radius varies by roughly 10% for this specific $(\dot{M},{M}_{{\rm{p}}},{R}_{{\rm{p}}})$ case over all $(\gamma ,\mu )$ combinations. In Section 7.3 we look at this more generally.

However, the size of the luminosity jump ΔL at the shock is very nearly the same across all simulations in Figure 9. Because in all cases the immediate upstream region is in the free-streaming regime, the small differences in $L({{r}_{\mathrm{shock}}}^{+})$ are related to the slightly different downstream luminosities or, equivalently, ${f}_{\mathrm{red}}({{r}_{\mathrm{shock}}}^{-})$.

6.1.3. Global Physical Efficiency

6.1.4. Other Quantities

6.1.5. Summary of the Effect of the EOS

Simulations such as those presented here do not spatially resolve the physical Zel'dovich spike (see Appendix B of Vaytet et al. 2013b), which is typically much thicker than the mean free path of the gas particles but still orders of magnitude smaller than a photon mean free path (Zel'dovich & Raizer 1967; Drake 2007). The peak temperatures should be very high but the cooling behind the peak also very quick, so that the fate of dust grains passing through this shock is not obvious a priori. Calculating their time-dependent sublimation and recondensation including nonequilibrium effects is beyond the scope of this work. Because our standard assumption is to use time-independent (equilibrium) abundances from the simple model of Isella & Natta (2005), in which the destruction temperature is simply ${T}_{\mathrm{dest}}=1220\times {\rho }_{-11}^{0.0195}$ K, where ${\rho }_{-11}\,\equiv \rho \times {10}^{11}\,{\mathrm{cm}}^{3}\,{{\rm{g}}}^{-1}$, we consider in this section the other extreme, namely that the grains are not destroyed. Note that the question of the dust (non-)destruction in the shock poses itself only for shocks at low enough temperature that solid grains are still present the incoming material (e.g., the $\approx 1\,{M}_{{\rm{J}}}$ cases at $\dot{M}={10}^{-3}\,{M}_{\oplus }\,{{\rm{yr}}}^{-1}$ in Figure 3).

Figure 10 shows the resulting shock structure for constant dust abundance in red and equilibrium abundances in blue. In all relevant quantities (ρ, v, T, L, ${f}_{\mathrm{red}}$), with the obvious exception of ${\kappa }_{{\rm{R}}}$ and ${\kappa }_{{\rm{P}}}$, the profiles are essentially identical for these two cases. This also holds throughout the simulation domain (not shown). In particular, the shock temperature is the same, and the Zel'dovich spike stays very optically thin (as it should), increasing from, very roughly, ${\rm{\Delta }}{\tau }_{{\rm{R}}}\sim {10}^{-5}$ to 10⁻². Note that the opacity is too high in the always-dust case (see blue solid line in opacity panel), but this does not affect the temperature structure of higher postshock regions, i.e., directly below the shock.

Zoom In Zoom Out Reset image size

As a more extreme case, we also switched off the dust contribution entirely, shown by gray profiles in Figure 10. This represents the limiting case of the reduction of Szulágyi et al. (2017) by a factor of 10 relative to the interstellar medium, on the grounds that the growth of dust grains into larger aggregates diminishes the opacity. As pointed out by Uyama et al. (2017), dust also tends to settle to the midplane while accretion comes from higher up in the circumstellar disk, so that the accretion flow onto a planet is actually possibly devoid of dust.

Leaving out the dust opacity lowered the Rosseland mean by four orders of magnitude near the shock but the Planck mean only by a factor of two due to the high Planck opacity of the gas. The consequence was a decrease in shock temperature by only 100 K, associated with a jump in the reduced flux increasing from ${\rm{\Delta }}{f}_{\mathrm{red}}\approx 0.7$ to ${\rm{\Delta }}{f}_{\mathrm{red}}\approx 1$.

Our conclusion from this test is that the upstream ${f}_{\mathrm{red}}$ (given ${f}_{\mathrm{red}}\approx 0$ downstream), and hence ultimately the opacity there, is important in setting the shock temperature, but that the dust destruction in the Zel'dovich spike is not. This seems to imply that one would need to follow the time-dependent evaporation of the dust approaching the shock in each simulation because the outer parts are always cold enough for dust to be present.¹² However, this will affect only a very small fraction of cases. Indeed, it requires (i) preshock temperatures above the (density-dependent) dust destruction temperature but (ii) not too high to not already have the dust destroyed at a large distance ahead from shock, in which case the dust would certainly be evaporated. The transition region is very narrow in temperature (ΔT ∼ 100 K; Semenov et al. 2003) compared to the range of shock temperatures we can expect (see, e.g., Figure 3). For these cases where it is relevant, a zeroth-order approach would be to compare the flow timescale ${t}_{\mathrm{flow}}=r/v$ to the evaporation time as given by the Polanyi-Wagner formula (see, e.g., Equation (20) of Grassi et al. 2017). However, this is an only marginally important consideration and will not be investigated further.

The main question driving this work is whether the processing of the accreting gas through the shock leads to hot starts or cold starts (Marley et al. 2007). While a detailed coupling of the shock results to formation calculations is needed to resolve this, there are several ways of estimating the outcome beforehand. We now consider them in turn.

7.1.1. Luminosity-based Argument: Shock Heating versus Internal Luminosity

The classical understanding of the extreme outcomes considers that the accreting gas carries only kinetic energy and that for cold starts, all of this incoming energy is radiated away at the shock, with no energy brought into the planet, while for hot starts, none of the incoming energy should leave the system but instead is added to the planet. However, this view neglects the thermal energy of the gas that is brought into the planet, which comes from preheating by the radiation from the shock. In Paper I we have shown that this is in fact not negligible, leading to the definition of ${\eta }^{\mathrm{phys}}$ instead of ${\eta }^{\mathrm{kin}}$ (see Section 4.4). Therefore, it is possible for most of the accretion energy to be radiated away at the shock, but to still have an important heating of the planet by the shock if the inward heating is high compared to the internal luminosity.

To quantify this, at least when following the global energetics alone, one can look at the effective heating by the shock ${Q}_{\mathrm{shock}}^{+}$ (to be defined shortly) relative to the internal luminosity ${L}_{\mathrm{int}}$,

$$\begin{eqnarray}&&{q}_{\mathrm{shock},\mathrm{rel}}^{+}\equiv \displaystyle \frac{{Q}_{\mathrm{shock}}^{+}}{{L}_{\mathrm{int}}}.\end{eqnarray} \tag{ 47 }$$

A high ${q}_{\mathrm{shock},\mathrm{rel}}^{+}$ would suggest that the shock can heat up the planet appreciably.

The effective heating of the planet by the shock, in turn, is the net rate of total energy that is not lost from the system and therefore added to the planet. By the definition of ${\eta }^{\mathrm{phys}}$ (Equation (10)), this is

$$\begin{eqnarray}&&{Q}_{\mathrm{shock}}^{+}=\left(1-{\eta }^{\mathrm{phys}}\right)| \dot{E}({r}_{\max })| =| \dot{E}({{r}_{\mathrm{shock}}}^{-})| ,\end{eqnarray} \tag{ 48 }$$

where ${{r}_{\mathrm{shock}}}^{-}$ is immediately downstream of the shock. From the definition of $\dot{E}$ (Equation (11)), we have in the isothermal limit that is applicable to our results

$$\begin{eqnarray}&&{Q}_{\mathrm{shock}}^{+}=\dot{M}\left[h({{r}_{\mathrm{shock}}}^{-})+{e}_{\mathrm{kin}}({{r}_{\mathrm{shock}}}^{-})\right]\,\end{eqnarray} \tag{ 49a }$$

$$\begin{eqnarray}&&=\,\left(1+\displaystyle \frac{\gamma -1}{2{\gamma }^{2}{{\mathscr{M}}}^{2}}\right)\dot{M}h({T}_{\mathrm{shock}})\,\end{eqnarray} \tag{ 49b }$$

$$\begin{eqnarray}&&=\,\left(1+\displaystyle \frac{\gamma -1}{2{\gamma }^{2}{{\mathscr{M}}}^{2}}\right)\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{{k}_{{\rm{B}}}{T}_{\mathrm{shock}}}{\mu {m}_{{\rm{H}}}}\dot{M}\,\end{eqnarray} \tag{ 49c }$$

$$\begin{eqnarray}&&=\,\left(\displaystyle \frac{2}{(\gamma -1){{\mathscr{M}}}^{2}}+\displaystyle \frac{1}{{\gamma }^{2}{{\mathscr{M}}}^{4}}\right){L}_{\mathrm{acc},\max },\end{eqnarray} \tag{ 49d }$$

using also that for an isothermal shock the post- and preshock densities are related by ${\rho }_{2}={\rho }_{1}\gamma {{\mathscr{M}}}^{2}$ and that ρ₂v₂ = ρ₁v₁. We recall that ${L}_{\mathrm{acc},\max }={{GM}}_{{\rm{p}}}\dot{M}/{R}_{{\rm{p}}}$. Equation 49(b) in particular shows that what is added to the planet is mostly the enthalpy of the gas, with a small additional term corresponding to the leftover kinetic energy (from the postshock settling velocity of the gas). The latter is vanishingly small in the limit of a high preshock Mach number.

The shock heating relative to ${L}_{\mathrm{acc},\max }$ is plotted against Mach number in Figure 11(a). It may seem surprising that for low Mach numbers ${\mathscr{M}}\lt 2$–4 (depending on γ), we have that ${Q}_{\mathrm{shock}}^{+}\gt {L}_{\mathrm{acc},\max }$. However, this simply comes from the fact that the gas at ${r}_{\max }$ does not bring only kinetic energy but also enthalpy with it.

Zoom In Zoom Out Reset image size

For the range ${\mathscr{M}}\approx 7$–35 (looking at Figures 4 and 8) and with γ = 1.44 as we used here, we see from Figure 11(a) that ${Q}_{\mathrm{shock}}^{+}\approx (10$%–0.4%$){L}_{\mathrm{acc}}$, respectively. In this high-Mach number limit (valid actually already for ${\mathscr{M}}\gtrsim 1$), the effective heating rate can be simplified to

$$\begin{eqnarray}&&\,{Q}_{\mathrm{shock}}^{+}=\displaystyle \frac{2}{\gamma -1}\displaystyle \frac{{L}_{\mathrm{acc}}}{{{\mathscr{M}}}^{2}}\,\end{eqnarray} \tag{ 50a }$$

$$\begin{eqnarray}&&\,=\,\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{{k}_{{\rm{B}}}}{\mu {m}_{{\rm{H}}}}{\left(\displaystyle \frac{{{GM}}_{{\rm{p}}}{\dot{M}}^{5}}{16\pi \sigma {{R}_{{\rm{p}}}}^{3}}\right)}^{1/4}{{\ell }}^{1/4}\,\end{eqnarray} \tag{ 50b }$$

$$\begin{eqnarray}\begin{array}{rcl} & = & 3\,\times {10}^{-4}\,{L}_{\odot }\,\displaystyle \frac{{\gamma }_{1.44}}{{\gamma }_{1.44}-1}{\left(\displaystyle \frac{\mu }{2.29}\right)}^{-1}{{\ell }}^{1/4}\\ & & \times {\left(\displaystyle \frac{\dot{M}}{{10}^{-2}{M}_{\oplus }{\mathrm{yr}}^{-1}}\right)}^{5/4}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-3/4}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{1/4},\end{array}\end{eqnarray} \tag{ 50c }$$

where again ${\ell }=1+{L}_{\mathrm{dnstr}}/{L}_{\mathrm{acc},\max }$, the downstream luminosity ${L}_{\mathrm{dnstr}}={L}_{\mathrm{int}}+{L}_{\mathrm{compr}}$ being the sum of the luminosity coming from the deep interior and the compression luminosity, and writing γ_1.44 ≡ γ/1.44. For Equations 50(b) and (c), we took the free-streaming limit for the shock temperature (Equation (33)), which was seen to hold over most of parameter space. As in Equation (33), ℓ formally depends on $(\dot{M},{M}_{{\rm{p}}},{R}_{{\rm{p}}})$, but this dependence is negligible in the limit of ${L}_{\mathrm{dnstr}}\ll {L}_{\mathrm{acc},\max }$.

Mordasini et al. (2017) have demonstrated that not only a "hot accretion", but even a "cold nominal" formation scenario leads to warm starts. Therefore, we estimate ${q}_{\mathrm{shock},\mathrm{rel}}^{+}$ with the cold starts of Marley et al. (2007), which represent a conservative lower limit to the entropies and luminosities of planets during their formation. Because Marley et al. (2007) do not report the internal luminosities of their planets during formation, we need to estimate them. Their Figure 4 shows that right after formation, the cold starts stay within 90% of the initial luminosity for a timescale ${t}_{90}\lesssim 1\,\mathrm{Myr}$ for 1 M_J, 1–2 Myr for 2 M_J, and ∼6 Myr and increasing for 4 M_J and up.¹³ The time spent in runaway accretion in their models is around ${t}_{\mathrm{acc}}=0.1\,\mathrm{Myr}$ to 0.5 Myr for 1 to 10 M_J, respectively. Thus even for ${M}_{{\rm{p}}}=1\,{M}_{{\rm{J}}}$, the cooling timescale ${t}_{90}$ is much longer than ${t}_{\mathrm{acc}}$, and it should be a reasonable approximation to assume that the internal luminosity and the radius do not change considerably between the last stages of the runaway accretion phase—when, for example, more than half the planet mass has been assembled—and right after.

Figure 11(b) displays the approximate heating by the shock (Equation 50(c)). Considering the extreme combinations of μ and γ as in Section 6, we find ${Q}_{\mathrm{shock}}^{+}\approx {10}^{-4}$–${10}^{-3}\,{L}_{\odot }$ for $\dot{M}\,\approx {10}^{-2}\,{M}_{\oplus }\,{\mathrm{yr}}^{-1}$ and ${Q}_{\mathrm{shock}}^{+}\approx {10}^{-5}$–${10}^{-4}\,{L}_{\odot }$ for ${10}^{-3}\,{M}_{\oplus }\,{\mathrm{yr}}^{-1}$. The radius is fixed at ${R}_{{\rm{p}}}=1.5\,{R}_{{\rm{J}}}$, a typical value for forming planets (Mordasini et al. 2012a), but varying it in a reasonable range ${R}_{{\rm{p}}}\approx 1$–3 R_J does not change the outcome qualitatively.

We compare in Figure 11(b) the shock heating to the internal luminosity of the Marley et al. (2007) cold starts. In their model, most of the mass is accreted with $\dot{M}\approx {10}^{-2}\,{M}_{\oplus }\,{\mathrm{yr}}^{-1}$ with a linear decrease of $\dot{M}$ toward the end (Hubickyj et al. 2005; Marley et al. 2007). Thus, the relevant heating is around ${Q}_{\mathrm{shock}}^{+}\gt {10}^{-4}\,{L}_{\odot }$. This is one to two orders of magnitude higher (!) than the (post-)formation luminosities ${L}_{\mathrm{int}}\,\approx {10}^{-6}\,{L}_{\odot }$ of Marley et al. (2007). At $1\,{M}_{{\rm{J}}}$, the heating could be only moderate, but for ${M}_{{\rm{p}}}\gtrsim 2\,{M}_{{\rm{J}}}$, the conclusion becomes more secure. Moreover, taking ${L}_{\mathrm{dnstr}}\ne 0$ into account for the estimate would only lead to a lower Mach number and thus a higher ${Q}_{\mathrm{shock}}^{+}$.

Thus, based on this a posteriori comparison of the internal luminosity and of the energy input rate, the shock is expected to heat up planets in the "cold classical" approach (Marley et al. 2007; Mordasini et al. 2017). The importance of the shock should increase with planet mass.

7.1.2. Shock-temperature-based Argument

7.1.3. Entropy-based Argument

Next we discuss the opacities. We have seen in detail in this work that the Rosseland mean controls the extent to which radiation is diffusing as opposed to freely streaming. As for the Planck opacity, it is an important factor in determining the extent to which the opacity carrier (gas or dust) and the radiation are in equilibrium. Indeed, the inflowing matter will equilibrate with the outgoing radiation, leading to ${T}_{\mathrm{gas}}\approx {T}_{\mathrm{rad}}$, if the energy exchange time (controlled by the absorption coefficient $\rho {\kappa }_{{\rm{P}}};$ see Malygin et al. 2017) is shorter than the time needed for the gas to reach the shock (controlled by ${v}_{\mathrm{ff}}$). Therefore we take a critical look at uncertainties concerning the opacities.

• 1.  
  We find that everywhere except in the shock (in the Zel'dovich spike), the preshock temperatures of the gas and radiation are equal. Simulations that use nonequilibrium radiation transport with tables with unrealistically low Planck mean values (see Figure 2(b)) might not find that the shock is able to preheat the gas. Thus it is an important result that the Planck values are high enough for the radiation and gas to be coupled.
• 2.  
  A quite uncertain aspect of dust opacity computations is the distribution of dust grain properties (size, porosity, composition, etc.) and also their sublimation. However, we found that the exact opacity in the accretion flow does not matter for the shock temperature. This effectively removes a source of uncertainty and makes the derived shock temperatures more robust.
• 3.  
  Nevertheless, the presence of dust in the accretion flow was seen to affect the temperature at and beyond the dust destruction front (see Stahler et al. 1980; Vaytet et al. 2013a) and thus also the luminosity at the Hill sphere. With a dust destruction front, the temperature there and beyond remains higher by up to a factor of several compared to the expression for a constant-$(L/{f}_{\mathrm{red}})$ profile, Equation (43) (although with the same power-law dependence). The decrease or increase in luminosity between the shock and the Hill sphere is also different from the case without a dust destruction front (see Section 6.1). However, if the dust in the incoming flow is strongly depleted relative to the interstellar medium abundance assumed in Semenov et al. (2003), the flow will tend more to be in the free-streaming regime, and its temperature is thus given by Equation (43).
• 4.  
  We conducted tests as in Paper I with constant low opacity. Typically, as we verified separately (not shown), at lower values ${\kappa }_{{\rm{R}}}={\kappa }_{{\rm{P}}}\sim 0.01$ cm² g⁻¹, the radiation and gas temperatures stay decoupled even in the high-density postshock region. This is entirely unrealistic, however, given that at these densities ($\rho \approx {10}^{-10}$–${10}^{-8}\,{\text{g cm}}^{-3}$) and temperatures (${T}_{\mathrm{gas}}\approx 500$–5000 K; see Figure 4 of Paper I) the Planck opacity is rather of order ${\kappa }_{{\rm{P}}}\sim 10$ cm² g⁻¹ (Figure 2(b)).
• 5.  
  Finally, we related the behavior of the opacity to that of the temperature and thus also of the luminosity (Section 5.1). This makes it now possible to understand the "bursts" in L seen, e.g., at 0.3 au in Figure 8 of Vaytet et al. (2013a), who also use FLD, but kept the frequency dependency. These bursts are associated with sharp opacity transitions in the respective wavelength band (in particular at the dust destruction front) and with slight changes of slope in the temperature.

We note that it is in the Zel'dovich spike that nonequilibrium (2-T) effects lead to the formation of observable spectral tracers of accretion onto protostars and brown dwarfs (Hartmann et al. 2016; Santamaría-Miranda et al. 2018). This emission is discussed for the shock onto the circumplanetary disk in Aoyama et al. (2018) and for the accretion shock on the planet surface in a forthcoming publication (Y. Aoyama et al. 2019, in preparation).

In these first two papers (Paper I, this work) we have restricted ourselves to a perfect EOS (constant μ and γ). To first order, this should not affect our main results. However, (i) the luminosity in the accretion flow (and thus at the Hill sphere), (ii) the postshock compression luminosity, and thus also (iii) the more precise value of ${T}_{\mathrm{shock}}$ (through the ${\rm{\Delta }}{f}_{\mathrm{red}}$ factor) should all be affected to some extent by the EOS, at least for some combinations of $(\dot{M},{M}_{{\rm{p}}},{R}_{{\rm{p}}})$.

Concerning item (i), we studied in Section 6.1.2 an estimate ${\rm{\Delta }}{L}_{(43)}$ of the drop or increase in the luminosity between the shock and the Hill radius. In the limit of a perfect EOS and of temperature profile $T\propto {r}^{-1/2}$, Equation (45) shows that the ratio $| {\rm{\Delta }}{L}_{(43)}| /{L}_{\mathrm{acc}}$ is highest for high accretion rates, low masses, and high radii. Over $\dot{M}={10}^{-3}$–${10}^{-2}\,{M}_{\oplus }\,{\mathrm{yr}}^{-1}$, ${M}_{{\rm{p}}}\,=1$–10 M_J, and ${R}_{{\rm{p}}}=1$–5 R_J (an even wider parameter space than what we consider for our simulations), the relative drop or increase is never larger than $| {\rm{\Delta }}{L}_{(43)}| /{L}_{\mathrm{acc}}\approx 10 \%$–15%, taking the extreme case of $(\mu =1.23,\gamma =1.1)$. Taking the actual temperature profile into account (as opposed to assuming T ∝ r^−1/2 everywhere) will change this somewhat, but usually only to make it smaller. In any case, the variation in L across the Hill sphere is unimportant compared to the effect of other simplifications of our model. We note that for molecular hydrogen and neutral helium (our default case), the actual change in L (i.e., as measured from the simulation and not using a simple T ∝ r^−1/2 profile) is less than 2% for any $(\dot{M},{M}_{{\rm{p}}},{R}_{{\rm{p}}})$ considered here.

We report already here that preliminary estimates suggest that for simulations with a full EOS, the change in luminosity across the Hill radius is at most approximately 20%, which is thus noticeable but also not large. Details will be presented in a forthcoming publication.

Finally, the relative smallness of the luminosity change ${\rm{\Delta }}{L}_{(43)}$ justifies a posteriori the assumption of constant $(L/{f}_{\mathrm{red}})$ made to derive it. Indeed, a relative change $| {\rm{\Delta }}{L}_{(43)}| /{L}_{\mathrm{acc}}\,\approx 10 \%$–15% over 2–3 dex in radius (from ${R}_{{\rm{p}}}\approx 1$–3 R_J to ${R}_{\mathrm{acc}}\,\approx {k}_{\mathrm{Lissauer}}{R}_{\mathrm{Hill}}({M}_{{\rm{p}}})\approx 250$–500 R_J for ${M}_{{\rm{p}}}\approx 1$–10 M_J at 5.2 au) corresponds to an approximate average slope (see Equation (30)) of at most $| \beta | \approx {\mathrm{log}}_{10}(1.15)/2\,=\,0.04$ if ${f}_{\mathrm{red}}$ is constant or $| \beta | \approx 0.4$ at most if ${f}_{\mathrm{red}}$ also changes by a factor of approximately 3, as in Figure 9. Thus even for these conservative estimates, we find $| \beta | \ll 2$, justifying a posteriori the assumption of a constant L/${f}_{\mathrm{red}}$ used to derive the change in L.