The Planetary Accretion Shock. II. Grid of Post-Shock Entropies and Radiative Shock Efficiencies for Non-Equilibrium Radiation Transport

TL;DR

This paper uses advanced radiation-hydrodynamical simulations to analyze the properties of accretion shocks in gas giant formation, revealing high post-shock entropies and efficiencies that favor warm or hot planetary starts.

Contribution

It provides a comprehensive grid of post-shock entropies and efficiencies using non-equilibrium radiation transport, advancing understanding of planetary formation shocks.

Findings

Shock temperature is mainly determined by the free-streaming limit.

Post-shock entropies are high, making cooling unlikely.

Radiation efficiencies are very high (>97%), supporting warm/hot start scenarios.

Abstract

In the core-accretion formation scenario of gas giants, most of the gas accreting onto a planet is processed through an accretion shock. In this series of papers we study this shock since it is key in setting the forming planet's structure and thus its post-formation luminosity, with dramatic observational consequences. We perform one-dimensional grey radiation-hydrodynamical simulations with non-equilibrium (two-temperature) radiation transport and up-to-date opacities. We survey the parameter space of accretion rate, planet mass, and planet radius and obtain post-shock temperatures, pressures, and entropies, as well as global radiation efficiencies. We find that usually, the shock temperature T_shock is given by the "free-streaming" limit. At low temperatures the dust opacity can make the shock hotter but not significantly. We corroborate this with an original semi-analytical derivation of T_shock . We also estimate the change in luminosity between the shock and the nebula. Neither T_shock nor the luminosity profile depend directly on the optical depth between the shock and the nebula. Rather, T_shock depends on the immediate pre-shock opacity, and the luminosity change on the equation of state (EOS). We find quite high immediate post-shock entropies (S ~ 13-20 kB/mH), which makes it seem unlikely that the shock can cool the planet. The global radiation efficiencies are high (eta^phys > 97%) but the remainder of the total incoming energy, which is brought into the planet, exceeds the internal luminosity of classical cold starts by orders of magnitude. Overall, these findings suggest that warm or hot starts are more plausible.

Figures

d of (Mordasini et al., [2017](#bib.bib64)). Most points are below ${L_{\textnormal{int}}}/{L_{\rm acc}}\sim 1$. Overall, the higher the accretion rate or the mass, the less important the interior luminosity. Using the hot-start population should yield somewhat weaker shock since the radii are larger (leading to the “core mass effect”; Mordasini, [2013](#bib.bib61)) but overall the results should be similar (see figure 13 of Mordasini et al., [2017](#bib.bib64)).

Equations

$$\frac{\partial \rho }{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\rho v\right)$$

$$\frac{\partial \left(\rho v\right)}{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\varrho v^2\right)+\frac{\partial P}{\partial r}$$

$$\rho \frac{\partial v}{\partial t}+\rho v\frac{\partial v}{\partial r}+\frac{\partial P}{\partial r}=-\varrho g,$$

$$\frac{\partial \left(E_{\text{kin}}+E_{\text{int}}\right)}{\partial t}+\nabla \left(\left[E_{\text{kin}}+E_{\text{int}}+P\right]v\right)$$

$$\frac{\partial E_{\text{rad}}}{\partial t}+\nabla F_{\mathrm{rad}}$$

$$F_{\mathrm{rad}}$$

$$D_{\text{F}}$$

$$R$$

$$\Lambda$$

$$=c{\kappa }_{\text{P}}\rho \left(a{T_{\text{gas}}}^4-E_{\text{rad}}\right).$$

$$10^{-3}⩽$$

$$1.3⩽$$

$$1.6\lesssim$$

$$\sigma {T_{\mathrm{shock}}}^4$$

$$\approx \frac{G{M}_{\text{p}}\dot{M}}{16\pi {R_{\text{p}}}^3},$$

$$f_{\text{red}}\equiv F_{\mathrm{rad}}/c{E}_{\text{rad}}$$

$$L_{\mathrm{acc}}=\frac{G{M}_{\text{p}}\dot{M}}{R_{\text{p}}}.$$

$$\displaystyle\Delta{\tau_{\mathrm{R}}}\sim 3\left(\frac{{\kappa_{\textnormal{R}}}}{1~{}\mbox{cm${}^{2}$\,g${}^{-1}$}}\right)$$

$$\times \sqrt{\left(\frac{1M_{\text{J}}}{M_{\text{p}}}\right)\left(\frac{2R_{\text{J}}}{r_{\mathrm{shock}}}\right)},$$

$${\eta }^{\mathrm{phys}}\equiv \frac{\dot{E}(r_{\max })-\dot{E}({r_{\mathrm{shock}}}^{-})}{\dot{E}(r_{\max })},$$

$$\dot{E}(r)\equiv -|\dot{M}|\left[e_{\text{kin}}(r)+h(r)+\Delta \Phi (r,r_{\mathrm{shock}})\right],$$

$$\Delta \Phi (r,r_0)=-G{M}_{\text{p}}\left(\frac{1}{r}-\frac{1}{r_0}\right).$$

$${\eta }_{\text{isoth}}^{\mathrm{phys}}$$

$$=\frac{\left({\gamma }^2{\mathcal{M}}^4-1\right)(\gamma -1)}{{\gamma }^2{\mathcal{M}}^2\left[(\gamma -1){\mathcal{M}}^2+2\right]},$$

$${\eta }_{\text{isoth}}^{\mathrm{kin}}=1-\frac{1}{{\gamma }^2{\mathcal{M}}^4}$$

$${\eta }^{\mathrm{kin}}\equiv \frac{\Delta F_{\mathrm{rad}}}{\frac{1}{2}{\rho }_{+}{v_{+}}^3},$$

$$P_{\mathrm{ram}}=\rho {v_{\text{ff}}}^2=\frac{\dot{M}}{4\pi {R_{\text{p}}}^2}{v}_{\text{ff}}=\frac{\sqrt{2G}}{4\pi }\frac{{M_{\text{p}}}^{1/2}\dot{M}}{{R_{\text{p}}}^{5/2}}.$$

$$ac{T_{\text{rad}}}^4$$

$$f_{\text{red}}$$

$$R\equiv \frac{1}{{\kappa }_{\text{R}}\rho }\left|\frac{\mathrm{d}\ln E_{\text{rad}}}{\mathrm{d}r}\right|$$

$$\frac{\mathrm{d}\ln {T_{\text{rad}}}^4}{\mathrm{d}\ln r}=-2+\frac{\mathrm{d}\ln L}{\mathrm{d}\ln r}-\frac{\mathrm{d}\ln f_{\text{red}}}{\mathrm{d}\ln r}.$$

$$\frac{\mathrm{d}\ln f_{\text{red}}}{\mathrm{d}\ln r}$$

$$\frac{\mathrm{d}\ln R}{\mathrm{d}\ln r}=\frac{\mathrm{d}}{\mathrm{d}\ln r}\left(\ln \left[\frac{\mathrm{d}\ln {T_{\text{rad}}}^4}{\mathrm{d}\ln r}\right]-\ln \left[{\kappa }_{\text{R}}\rho r\right]\right)$$

$$\frac{\mathrm{d}\ln R}{\mathrm{d}\ln r}=$$

$$\lambda (R)=\frac{1}{3+R},$$

$$1+\frac{\mathrm{d}\ln \lambda }{\mathrm{d}\ln R}=3\lambda .$$

Full text

accretion — planets and satellites: formation — planets and satellites: gaseous planets — methods: numerical — methods: analytical — radiative transfer

††software: PLUTO (Mignone et al., 2007, 2012), Makemake (Kuiper et al., 2010; Kuiper, Yorke, & Mignone, in prep.)

1 Introduction

With its first direct detections already some ten to fifteen years ago (Chauvin et al., 2004; Marois et al., 2008), the technique of direct imaging has started to reveal a scarce but interesting population of planets or very-low-mass substellar objects at large separations from their host stars (Bowler, 2016; Bowler & Nielsen, 2018; Wagner et al., 2019). The formation mechanism of individual detections is often not obvious but gravitational instability as well as core accretion (with the inclusion of $N$-body interactions during the formation phase and in the first few million years afterwards) are likely candidates to explain the origin of at least some of these systems (e.g., Marleau et al., 2019). These different formation pathways may imprint into the observed brightness of the planets (Baruteau et al., 2016; Mordasini et al., 2017).

To interpret the brightness measurements however requires knowing the post-formation luminosity of planets of different masses. Formation models, principally the ones of the California group (Pollack et al., 1996; Bodenheimer et al., 2000; Marley et al., 2007; Lissauer et al., 2009; Bodenheimer et al., 2013) and of the Bern group (Alibert et al., 2005; Mordasini et al., 2012b, a, 2017) seek to predict this luminosity within the approximation of spherical accretion. They need to assume something about the efficiency of the gas accretion shock at the surface of the planet during runaway gas accretion. This efficiency is defined as the fraction of the total energy influx which is re-radiated into the local disc and thus does not end up being added to the planet. The extremes are known as “cold starts” and “hot starts” (Marley et al., 2007) and their post-formation luminosities can differ by orders of magnitude.

In a recent series of papers, Berardo et al. (2017), Berardo & Cumming (2017), and Cumming et al. (2018) have begun calculating the structure of accreting planets following Stahler (1988). Crucially, they take into account that in the settling zone below the shock, the continuing compression of the post-shock layers leads to a non-constant luminosity. They find that the thermal influence of the shock on the evolution of the planet during accretion depends on the contrast between the entropy of the (outer convective zone of the) planet and that of the post-shock gas. This approach promises eventually to lead to more realistic predictions of the post-formation luminosity111 This statement holds given an accretion history, which, admittedly, is however uncertain since it depends in part on the migration behaviour of the planet, itself fraught with uncertainty, but does require, as a boundary condition, knowledge of the temperature of the shock.

While global three-dimensional (radiation-)hydrodynamical simulations of the protoplanetary and of the circumplanetary discs (Klahr & Kley, 2006; Machida et al., 2008; Tanigawa et al., 2012; D’Angelo & Bodenheimer, 2013; Szulágyi et al., 2016, 2017) have the potential of providing a realistic answer as to the post-shock conditions and its efficiency, they still have a limited spatial resolution of $\Delta x\sim 1~{}{R_{\textnormal{J}}}$ at best at the position of the planet, despite their high dynamical range of spatial scales. Also due to the computational cost, they are (currently) unable to survey the large input parameter space, which covers a factor of a few in planet radius, an order of magnitude in mass, and several orders or magnitude in accretion rate and internal luminosity (e.g., Mordasini et al., 2012b, 2017).

In the present series of papers, we use one-dimensional models of the gas accretion to take a careful look at shock microphysics. In Marleau et al. (2017, hereafter Paper I), we introduced our approach, presented a detailed analysis of results for one combination of formation parameters, and discussed the shock efficiency for a certain range of parameters. However, we restricted ourselves to equilibrium radiation transport, in which the gas and radiation temperatures are assumed to be equal everywhere, and assumed a perfect222 This is also sometimes termed a “constant EOS” but should not be refered to only as an “ideal gas”, as is unfortunately often done. Indeed, the latter only needs fulfill $P=\rho/(\mu{m_{\textnormal{H}}}){k_{\textnormal{B}}}T$ with $\mu$ not necessarily constant. A non-ideal EOS also includes for example quantum degeneracy effects.

equation of state (EOS), i.e., a constant mean molecular weight $\mu$ and heat capacity ratio333 For a perfect EOS, the various adiabatic indices $\Gamma_{1,2,3}$, the heat capacity ratio $\gamma$, as well as ${\gamma_{\rm eff}}\equiv P/{e_{\textnormal{int}}}+1$, where ${e_{\textnormal{int}}}$ is the internal energy, are all equal. Therefore we always write $\gamma$ in this paper. $\gamma$. We found that the shock is isothermal and supercritical and that the efficiencies can be as low as 40 %.

Here, we relax the assumption of equilibrium radiation transport, look more carefully at the role of opacity, and consider a wider parameter space than in Paper I, exploring systematically the dependence on accretion rate, mass, and planet radius. We also extend considerably our analytical derivations. We again assume a perfect EOS, varying $\mu$ and $\gamma$, and focus on the cases in which the internal luminosity ${L_{\textnormal{int}}}$ is much smaller than the shock luminosity ${L_{\rm acc}}$.

This paper is structured as follows: In Section 2 we estimate in what regions of $({\dot{M}},{M_{\textnormal{p}}},{R_{\textnormal{p}}},{L_{\textnormal{int}}})$, where ${\dot{M}}$ is the accretion rate while ${M_{\textnormal{p}}}$ and ${R_{\textnormal{p}}}$ are the planet mass and radius, the assumption ${L_{\textnormal{int}}}\ll{L_{\rm acc}}$ holds, which guides our choice of parameters. Section 3 reviews briefly our set-up and details the relevant microphysics, including the updated opacities. The main thrust of this paper is in Section 4, which presents and analyses results, including shock temperatures, global radiative shock efficiencies, and post-shock entropies, for a grid of simulations. In Section 5 we present semi-analytical derivations of the shock temperature and of the temperature profile in the accretion flow and compare these to our results. In Section 6 we investigate carefully the effect of different perfect EOS in non-equilibrium radiation transport and of dust destruction in the Zel’dovich spike. This motivates us to derive analytically the drop or increase of the luminosity across the Hill sphere. While we do not calculate the structure of forming planets using our shock results yet, Section 7 explores whether hot or cold starts are expected, and presents a further discussion. Finally, Section 8 summarizes this work and presents our conclusions.

2 Estimate of negligible internal luminosity

The main formation parameters are the accretion rate onto the planet, the planet mass, the planet radius, and the internal luminosity of the planet, denoted respectively by ${\dot{M}}$, ${M_{\textnormal{p}}}$, ${R_{\textnormal{p}}}={r_{\rm shock}}$ (the position of the shock ${r_{\rm shock}}$ defining the radius of the very nearly hydrostatic protoplanet), and ${L_{\textnormal{int}}}$. In Marleau et al. (2017) and this work, we focus on the case of negligible ${L_{\textnormal{int}}}$. This luminosity comes from the contraction and cooling of the planet interior and is generated almost entirely within a small fraction of the planet’s volume, where most of the mass resides. Specifically, this means ${L_{\textnormal{int}}}\ll\Delta L$, where the luminosity jump at the shock $\Delta L\sim{L_{\rm acc}}$. This reduces the number of free parameters in our study but, mainly, is also expected to be the limit in which the shock simulations we perform are the most relevant.

We can estimate in what cases neglecting the interior luminosity should be best justified. In the Marley et al. (2007) formation calculations, which represent the extreme case of cold accretion (Mordasini et al., 2017), it is blatantly obvious that ${L_{\textnormal{int}}}$ is negligible compared to ${L_{\rm acc}}$ (see their figure 3 and the discussion in Section 7.1). For more “moderate” (i.e., warmer) versions of cold starts, Figure 1 shows the ratio ${L_{\textnormal{int}}}/{L_{\rm acc}}$ for cold-start population synthesis planets as in figure 7d of (Mordasini et al., 2017). Most points are below ${L_{\textnormal{int}}}/{L_{\rm acc}}\sim 1$. Overall, the higher the accretion rate or the mass, the less important the interior luminosity. Using the hot-start population should yield somewhat weaker shock since the radii are larger (leading to the “core mass effect”; Mordasini, 2013) but overall the results should be similar (see figure 13 of Mordasini et al., 2017).

If we focus on ${\dot{M}}\gtrsim 10^{-3}~{}{M_{\oplus}}\,\textrm{yr}^{-1}$ and masses above a few ${M_{\textnormal{J}}}$, simulations with ${R_{\textnormal{p}}}\approx 1.5$–3 ${R_{\textnormal{J}}}$ should be the ones in which the accretion shock is the most relevant. At ${\dot{M}}\lesssim 10^{-4}~{}{M_{\oplus}}\,\textrm{yr}^{-1}$, the interior luminosity frequently represents an appreciable fraction of the accretion luminosity, even up to high planets masses. Therefore, we do not consider these lower rates in this paper. Of course, this estimate is not self-consistent since these formation and evolution calculations assume a fixed shock efficiency $\eta=100\,\%$, whereas we find clearly smaller values (in the sense that the heating of the planet by the post-shock material is important relative to the planet’s interior luminosity; see Section 7.1 and Paper I). Nevertheless, the estimate should provide reasonable guidance.

3 Simulation approach

As in Paper I, we use the PLUTO code (Mignone et al., 2007, 2012) to solve in a time-explicit fashion for the hydrodynamics. The one-dimensional, spherically-symmetric mass and momentum conservation equations are respectively

[TABLE]

where $P$ is the pressure and the local gravitational acceleration is $g=|-G{M_{\textnormal{p}}}/r^{2}|$, the self-gravity of the gas being negligible. Note that, using mass conservation, the momentum equation can also be written in the form

[TABLE]

i.e., without the geometry factors $r^{\pm 2}$ in the flux term.

The energy equation is similar to in Paper I but we make use of the newest version of the flux-limited diffusion (FLD) solver Makemake. We return to this in Section 3.2.

3.1 Set-up

The simulation domain extends from close to the planet’s surface almost out to the “accretion radius” ${R_{\textnormal{acc}}}$ (Bodenheimer et al., 2000; Paper I), at one third (${k_{\textnormal{Lissauer}}}=1/3$ Lissauer et al., 2009) of the Hill radius. A grid with 1000 uniformly-spaced cells between the inner edge of the domain, ${r_{\rm min}}$, and ${r_{\rm min}}+0.5~{}{R_{\textnormal{J}}}$ and 1000 geometrically stretched cells from ${r_{\rm min}}+0.5~{}{R_{\textnormal{J}}}$ to the outer edge of the domain, ${r_{\rm 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 post-shock structure.

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

We use a CFL number of 0.8 (we verified that $\textrm{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 free-fall 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 which defines the surface of the planet at ${r_{\rm shock}}={R_{\textnormal{p}}}$ (both variables are used interchangeably throughout, with ${r_{\rm shock}}^{-}$ (${r_{\rm shock}}^{+}$) denoting the downstream (upstream) location). This shock surface often, but not always, moves outwards over time, although at a rate which is slow even compared to the strongly sub-sonic post-shock velocity. The deepest pressure in the atmosphere we simulate is greater than the post-shock (ram) pressure by one to several orders of magnitude, depending on the amount of accumulated mass and the local gravitational acceleration $g=G{M_{\textnormal{p}}}/r^{2}$, where $r$ is the radial coordinate in spherical geometry and $G$ 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\times 10^{6}$ 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^{4}$–$10^{6}$ yr. (Even after ${t_{\textnormal{stop}}}=2\times 10^{7}$ s, our usual stopping time, some starting structures were still visible at a position $r$ deep below the shock for which ${t_{\textnormal{stop}}}\lesssim-\int_{r}^{{r_{\rm shock}}}(1/v){\rm d}x$, in at least some simulations; however, these regions are not of interest here.)

3.2 Radiation transport

A significant improvement since Paper I is the change from equilibrium to non-equilibrium (or two-temperature, 2-$T$) flux-limited diffusion (FLD) in the radiation transport routine Makemake (Kuiper et al., 2010; Kuiper, Yorke, & Mignone, in prep.) 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; Kuiper et al., 2010; Klassen et al., 2014; Commerçon et al., 2011b)

[TABLE]

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

[TABLE]

where ${D_{\textnormal{F}}}$ is the diffusion coefficient, $\lambda$ the flux limiter (see Section 5.1 below), ${\kappa_{\textnormal{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 density444In angularly-resolved radiation transport methods, it is in fact possible for the radiation to flow up the ${E_{\textnormal{rad}}}$ gradient (McClarren & Drake, 2010; Jiang et al., 2014). However, this can occur only for subcritical shocks and does not represent a large effect, thus not affecting our work.. In Equation (3), the combined cooling and heating or exchange term is

[TABLE]

We use the second line, since we take the source function $S$ to be the Planck function $B={\sigma}{T_{\textnormal{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 (LTE). Our grey approximation suffices to make ${\kappa_{\textnormal{P}}}$ and ${\kappa_{\textnormal{E}}}$, the Planck (blackbody-weighted) and ${E_{\textnormal{rad}}}$-weighted mean opacities respectively, be the same. Equations (3a) and (3b) are followed separately to prevent the $\Lambda$ terms from cancelling. This system states that the material (gas or dust, assumed here to have the same temperature) is losing energy at the rate $c{\kappa_{\textnormal{P}}}\rho a{T_{\textnormal{gas}}}^{4}$ but absorbing energy (photons) at the rate $c{\kappa_{\textnormal{E}}}\rho{E_{\textnormal{rad}}}$, and conversely for the radiation.

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

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

3.3 Microphysics: opacities

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 modified slightly the opacity routine555Original version available under http://www2.mpia-hd.mpg.de/home/henning/Dust_opacities/Opacities/opacities.html. to leave the evaporation of the most refractory component to be handled in the main code and to not add the contribution from the Helling et al. (2000) gas opacities.

Instead, 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 towards higher constant $\mu$ and $\gamma$, to evaluate the single-temperature ${\kappa_{\textnormal{P}}}$ at the radiation temperature ${T_{\textnormal{rad}}}\equiv({E_{\textnormal{rad}}}/a)^{1/4}$ and not at the gas temperature ${T_{\textnormal{gas}}}$. Fortunately, it is also justified. Indeed, looking at ${\kappa_{\textnormal{P}}}({T_{\textnormal{rad}}},{T_{\textnormal{gas}}})$ for fixed densities, one generally incurs a smaller mistake when using ${\kappa_{\textnormal{P}}}\approx{\kappa_{\textnormal{P}}}({T_{\textnormal{rad}}},{T_{\textnormal{rad}}})$ than with ${\kappa_{\textnormal{P}}}\approx{\kappa_{\textnormal{P}}}({T_{\textnormal{gas}}},{T_{\textnormal{gas}}})$. We also evaluate the dust opacities at ${T_{\textnormal{rad}}}$ but note that this makes barely a difference since we will find that the radiation and gas are always well coupled (${T_{\textnormal{gas}}}={T_{\textnormal{rad}}}$) in the region where the dust opacity matters (${T_{\textnormal{gas}}}\lesssim 1500$ K).

We remind 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_{\textnormal{P}}}$ is shown in Figure 2 and compared to the Rosseland mean ${\kappa_{\textnormal{R}}}$. In the low-temperature, dust-dominated regime ${\kappa_{\textnormal{P}}}$ is similar to ${\kappa_{\textnormal{R}}}$ with ${\kappa_{\textnormal{P}}}\sim{\kappa_{\textnormal{R}}}\sim 1$–10 cm2 g*-1*. For comparison, the different models of Semenov et al. (2003) are shown, with the exception of the porous models. Indeed, fluffy dust aggregates are not expected in protoplanetary discs due to compactification from collisions, as seems to be borne out by observations (Cuzzi et al., 2014; Kataoka et al., 2014, and references therein).

While the Rosseland mean opacity past dust evaporation (near $T\approx 1500$ K) drops to ${\kappa_{\textnormal{R}}}\approx 10^{-2}$ cm2 g*-1* (see figure 1 of Paper I), the Planck mean does not drop much below ${\kappa_{\textnormal{P}}}\approx 1$ cm2 g*-1* and even increases (between a few thousand to $10^{4}$ K) to ${\kappa_{\textnormal{P}}}=10^{2}$–104 cm2 g*-1*, depending on the density.

Notice 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 non-equilibrium radiation transport with the Bell & Lin (1994) Rosseland mean as their Plank opacity are effectively assuming much less coupling between the opacity carrier (dust or gas) and the 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 however depend also on the local density and velocity, as briefly discussed in Section 7.2.

4 Grid of simulations: results and analysis

We now present and discuss results for a grid of simulations for the macrophysical parameters

[TABLE]

We consider a mixture of molecular hydrogen and helium with helium mass fraction $Y=0.25$, yielding $\mu=2.29$ and $\gamma=1.44$. At high temperatures, the hydrogen should dissociate (Szulágyi & Mordasini, 2017) but we defer simulations taking this into account to another paper in this series (Marleau et al., in prep.).

The results are summarised in Figures 3, 4, and 5 and discussed in the following: Figure 3 shows, as a function of the macrophysical parameters, the resulting shock temperature (Section 4.1), the luminosity at roughly the Hill radius (Section 4.2), and the optical depth to the Hill radius (Section 4.3), whereas Figure 4 shows the global physical efficiency as a function of the pre-shock Mach number (Section 4.4). Finally, in Figure 5 we display, as a function of the macrophysical parameters, the efficiency (Section 4.5) and the post-shock entropy (Section 4.6).

Note that the grid is irregular in shock position because we considered the same ${r_{\rm min}}$ values for all masses and accretion rates but the shock moves at different rates ${\rm d}{r_{\rm shock}}/{\rm d}t$; over the course of $2\times 10^{7}$ s, which we use as the maximal simulation time since it is more than enough for the profiles to reach a quasi-steady state, the ranges of radial positions which the shock covers often do not overlap between different $({\dot{M}},{M_{\textnormal{p}}},{r_{\rm min}})$. We added several simulations at higher ${r_{\rm min}}\leqslant 2.9~{}{R_{\textnormal{J}}}$ for $({\dot{M}}=10^{-3}~{}{M_{\oplus}}\,\mathrm{yr}^{-1},{M_{\textnormal{p}}}=1.3~{}{M_{\textnormal{J}}})$.

4.1 Shock temperature

For our choice of macrophysical parameters, the shock temperature is always above 1000 K. In Figure 3, simulations with a given $({\dot{M}},{M_{\textnormal{p}}})$ but different ${r_{\rm min}}$ (symbol size) are seen to lead to the same shock temperature at a given radius ${r_{\rm shock}}={R_{\textnormal{p}}}$. This confirms that, even though the post-shock 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 (Kuiper et al., in prep.) does not modify the simulation outcomes.

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

[TABLE]

where $\eta^{\rm kin}\equiv\Delta{F_{\rm rad}}/(0.5\rho v^{3})\approx{L_{\rm acc}}/(G{M_{\textnormal{p}}}{\dot{M}}/{R_{\textnormal{p}}})$ is the normalised jump in luminosity and $\Delta{f_{\textnormal{red}}}\equiv{f_{\textnormal{red}}}({r_{\rm shock}}^{+})-{f_{\textnormal{red}}}({r_{\rm shock}}^{-})$ the jump in ${f_{\textnormal{red}}}$ at the shock. This equation will be revisited in Section 5.3. The free-fall density ${\rho_{\textnormal{ff}}}={\dot{M}}/(4\pi r^{2}{v_{\textnormal{ff}}})$ is evaluated ahead of the shock, with ${v_{\textnormal{ff}}}=\sqrt{2G{M_{\textnormal{p}}}/r}$ the approximate free-fall velocity (Paper I). The reduced flux

[TABLE]

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

As a consequence of this, Equation (6b) holds very accurately for almost all simulations; we are almost always in the limit of equation (28a) of Paper I, discussed as Equation (33) below. There is an exception to this, namely for the lowest mass (1.3 ${M_{\textnormal{J}}}$) at a low accretion rate ($10^{-3}~{}{M_{\oplus}}\,\mathrm{yr}^{-1}$) and towards larger radii (${R_{\textnormal{p}}}\gtrsim 2.5~{}{R_{\textnormal{J}}}$). In this case, the post-shock temperature is higher than predicted from Equation (6b) by $\Delta T\approx 50$ K. This is because ${f_{\textnormal{red}}}$ is lower ahead of the shock, with e.g. ${f_{\textnormal{red}}}^{+}\equiv{f_{\textnormal{red}}}({r_{\rm shock}}^{+})=0.65$ for the ${T_{\rm shock}}\approx 1100$ K case. Indeed, since ${F_{\rm rad}}\propto{f_{\textnormal{red}}}{T_{\textnormal{rad}}}^{4}$ and we find ${T_{\textnormal{gas}}}={T_{\textnormal{rad}}}$, a lower ${f_{\textnormal{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_{\textnormal{red}}}$ can also be pictured as a smaller effective speed of light, so that ${E_{\textnormal{rad}}}$ must increase so as to have the same ${F_{\rm rad}}={c_{\textnormal{eff}}}{E_{\textnormal{rad}}}$. We return to these points in Section 5.3 and show that also another effect is at play. Note however that this is the gas temperature but not the effective temperature, with only the latter setting the spectral shape.

4.2 Luminosity at the Hill sphere

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

[TABLE]

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_{\textnormal{acc}}}$. This radius is however much larger than the planet radius since we are considering the detached phase during giant planet formation (Mordasini et al., 2012b).

We also produced a grid with $(\mu=1.1,\gamma=1.1)$ (not shown), which increases the contrast with the present situation. One example will be discussed in Section 6.1 below. Over the grid, though, the shock temperatures were the same, as one would expect from Equation (6) since they do not depend on the EOS explicitly. The situation could be different for an ideal but non-perfect EOS, however, 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 $\approx 10\,\%$ compared to the immediate shock upstream luminosity. The luminosity profile is discussed in Section 6.1.

4.3 Optical depth of the infalling gas

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

[TABLE]

shown as dotted lines in Figure 3c. The agreement with the nominal models is not too rough (to at worse 1 dex) considering that we took a constant ${\kappa_{\textnormal{R}}}=1$ cm2 g*-1* in Equation (9) in this estimate for all simulations. Looking at the radial profiles, the maximum value of the actual ${\kappa_{\textnormal{R}}}(r)$ is near 7 (3) cm2 g*-1* for ${\dot{M}}=10^{-3}$ ($10^{-2}$) ${M_{\oplus}}$ yr*-1*. However, contrary to what was explained in Paper I, what is relevant is the immediate pre-shock opacity, as we shall see in Section 5.

The significance of Figure 3c can be appreciated only in conjunction with panels (a) and (b). It shows that the total optical depth, at least up to moderate depths of $\Delta{\tau_{\mathrm{R}}}\sim 10$, does not set the shock temperature nor the shock luminosity (which is nearly identical to the Hill-radius luminosity in the present case). In fact, the deviations of ${T_{\rm shock}}$ from Equation (6) do not occur at the highest optical depths. (What causes these deviations is explained in Section 5.) Also, when, as is often the case, the layers between the $\Delta{\tau_{\mathrm{R}}}\sim 1$ surface and the shock are not sufficiently diffusive (${f_{\textnormal{red}}}\lesssim 0.1$), it does not hold that ${T_{\rm shock}}^{4}\approx\frac{3}{4}(\Delta{\tau_{\mathrm{R}}}(r)+\frac{2}{3}){T_{\textnormal{eff}}}^{4}$, where ${T_{\textnormal{eff}}}$ is the temperature at the radius where $\Delta{\tau_{\mathrm{R}}}=2/3$, the photosphere, and where $\Delta{\tau_{\mathrm{R}}}(r)=-\int_{{r_{\rm max}}}^{r}\rho{\kappa_{\textnormal{R}}}\,{\rm d}r^{\prime}$ is the optical depth measured from the outer radius inwards. 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 ${\rm d}\tau/{\rm d}r={\kappa_{\textnormal{R}}}\rho$ is small.

Finally, we note that there is nothing special about the $\Delta{\tau_{\mathrm{R}}}\sim 1$ surface. As Figure 3c shows, most simulations at ${\dot{M}}=10^{-3}~{}{M_{\oplus}}$ yr*-1* remain optically thin or barely optically thick above the shock, yet their temperature structure (not shown) is qualitatively identical to cases with an optically thick accretion flow. In any case, $\Delta{\tau_{\mathrm{R}}}$ is not well defined due to its dependence on the outer integration radius. We come back to considerations of shock temperature in Section 5.3.

4.4 Shock efficiency against Mach number

When calculating a planet structure for given $({\dot{M}},{M_{\textnormal{p}}},{R_{\textnormal{p}}})$, only the post-shock $({P_{\textnormal{post}}},{T_{\rm 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 $\rho$; from this and ${\dot{M}}$, the boundary condition on the velocity $v$ and thus the luminosity profile $L(r)$ in the post-shock region follow (Berardo et al., 2017). When however following only the global energetics (see the review in sect. 2.1 of Berardo et al., 2017), we argued in Paper I that one should use $\eta^{\rm phys}$ and not $\eta^{\rm kin}$.

We show in Figure 4 the “global physical efficiency” $\eta^{\rm phys}$ of the shock against the pre-shock Mach number. The quantity $\eta^{\rm phys}$ is measured as ( Paper I, their equation (18))

[TABLE]

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

[TABLE]

where ${e_{\textnormal{kin}}}=\frac{1}{2}v^{2}$, ${e_{\textnormal{int}}}$, and $h={e_{\textnormal{int}}}+P/\rho$ are respectively the kinetic energy, internal energy density, and the enthalpy per unit mass, and $\Phi$ is the external potential. The $\Delta\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_{0}$ to $r$ given by

[TABLE]

We found that for all simulations here, using the tenth cell below the shock as the post-shock location was a robust prescription. The shock itself was identified by the ${\rm d}v/{\rm d}r<0$ and $\Delta P/\min(P)>5$ criterion in appendix B of Mignone et al. (2012).

The efficiency $\eta^{\rm 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)),

[TABLE]

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

[TABLE]

This definition was derived from energy conservation for the 1-$T$ case but it is also meaningful here since 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

[TABLE]

where $\Delta{F_{\rm rad}}$ is the jump in radiative flux at the shock and $\rho_{+}$ and $v_{+}$ the density and velocity upstream of the shock. This measured $\eta^{\rm kin}$ was found to match Equation (14). As mentioned above, since ${\mathcal{M}}\gtrsim 2.5$ and the shock is isothermal, we have $\eta^{\rm kin}\approx 100$ %. Thus, locally at the shock, the whole incoming kinetic energy is converted to radiation.

4.5 Shock efficiency against formation parameters

Figure 5 shows the global physical efficiency $\eta^{\rm phys}$ but now as an explicit function of the (possibly observable, macrophysical) formation parameters $({\dot{M}},{M_{\textnormal{p}}},{R_{\textnormal{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 phase666 This can easily be verified for instance with the results of the Bern planet formation code in the “Evolution” section of the Data Analysis Centre for Exoplanets (DACE) platform under https://dace.unige.ch by plotting ${R_{\textnormal{p}}}(t)$ against ${M_{\textnormal{p}}}(t)$. . Since the gas accretion onto a planet usually slows down with time, at least in the single-embryo-per-disc simulations of Mordasini et al. (2012b), our simulations clearly suggest that the efficiency $\eta^{\rm 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^{\rm phys}$ that is small early after detachment should lead to the accretion of “hot” (high-entropy) material into the planet. If the planet’s accretion timescale is much shorter than its global cooling timescale—i.e., ${t_{\rm acc}}\equiv{M_{\textnormal{p}}}/{\dot{M}}\ll{t_{\rm KH}}\equiv G{M_{\textnormal{p}}}^{2}/({R_{\textnormal{p}}}{L_{\textnormal{p}}})$, where ${t_{\rm acc}}$ and ${t_{\rm KH}}$ are respectively the accretion and Kelvin–Helmholtz timescales, 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^{\rm 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_{\rm acc}}$ and ${t_{\rm KH}}$ (i.e., for the planet as a whole).

4.6 Post-shock entropy

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); Cumming et al. (2018). Since these planet structure calculations always use a full EOS including chemical reactions (dissociation and ionisation), 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, since ${T_{\rm 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 (Marleau et al., in prep.).

To calculate the entropy, we can use the ideal-gas form $P=\rho/(\mu{m_{\textnormal{H}}}){k_{\textnormal{B}}}T$ (but with variable $\mu$) since degeneracy starts being relevant only at a conservative limit of $\rho\sim 10^{-2}$ g cm*-3* (cf. figure 1 of SCvH), several orders of magnitude above even the highest post-shock densities ${\rho_{\textnormal{post}}}\propto{\rho_{\textnormal{pre}}}{\mathcal{M}}^{2}$, with ${\rho_{\textnormal{pre}}}\sim 10^{-13}$–$10^{-10}$ g cm*-3* (see figure 4 of Paper I) and ${\mathcal{M}}\ll 100$, so that ${\rho_{\textnormal{post}}}\ll 10^{-6}$ g cm*-3*. We use the Saha equation and Sackur–Tetrode formula as implemented777See https://github.com/andrewcumming/gasgiant.

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)\ln 2=0.52~{}{k_{\textnormal{B}}}\,{m_{\textnormal{H}}}^{-1}$, where $Y=0.25$ is the helium mass fraction, than the ones 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 post-shock entropies ${s_{\rm post}}$ are shown in Figure 5. For the range of shock positions ${r_{\rm shock}}\approx 1.5$–3 ${R_{\textnormal{J}}}$ and masses ${M_{\textnormal{p}}}=1.3$–10 ${M_{\textnormal{J}}}$ shown, the entropies are, for ${\dot{M}}=10^{-3}~{}{M_{\oplus}}\,\textrm{yr}^{-1}$, mostly around ${s_{\rm post}}\approx 12$–14 but go up to ${s_{\rm post}}\approx 19$ (dropping, also in the following, the usual units of ${k_{\textnormal{B}}}\,{m_{\textnormal{H}}}^{-1}$ for clarity). For ${\dot{M}}=10^{-2}~{}{M_{\oplus}}\,\textrm{yr}^{-1}$, the whole range ${s_{\rm post}}\approx 13$–20 is covered. These values are high compared to the post-formation entropy of planets, which is at most around 10–14 ${k_{\textnormal{B}}}\,{m_{\textnormal{H}}}^{-1}$ according to current, though not definitive, predictions (Berardo & Cumming, 2017; Berardo et al., 2017; Mordasini, 2013; Mordasini et al., 2017). However, we caution and emphasize that this post-shock entropy is not the same as the entropy below the post-shock 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) post-shock entropy analytically, we briefly discuss the ram pressure

[TABLE]

We inserted the expressions for free-fall from infinity and find that Equation (16) holds very well for all simulations, even for those for which $\Delta{f_{\textnormal{red}}}\neq 1$. At high shock temperatures, the small pressure build-up ahead of the shock slows down the gas 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 pre-factor) dependence of the entropy $s(P,T)$ on $P$, outside of dissociation or ionisation regions, this will not be an important source of inaccuracy. The ram pressure varies from $P_{\rm ram}\approx 10^{-4}$ bar to $0.2$ bar for the range of parameters discussed here.

We compare the post-shock entropies ${s_{\rm post}}$ in Figure 5 based on the actual $T$ to ${s_{\rm post}}$ using as input $T={T_{\rm shock}}$ from Equation (6b), i.e., taking $\eta^{\rm kin}=\Delta{f_{\textnormal{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.

5 Analytics of the temperature

ahead of and at the shock

As will be seen in Figure 9, the temperature profile upstream of the shock shows variations which appear related to variations in opacity. This modifies the temperature at the outer edge of the accretion flow (i.e., the local nebula temperature) compared to what a naive extrapolation $T(r)={T_{\rm shock}}\left({r_{\rm shock}}/r\right)^{-1/2}$ would predict. Also, and even more importantly, the temperature at the shock is obviously a key outcome of our simulations. The usual expression for the shock temperature, with our modification of a factor $(1/4)^{1/4}$ (Equation (6)), provides in general a very good estimate. However, we have seen in Figure 3 that there can be small deviations. We now turn to the task of understanding both the temperature profile in the accretion flow and the shock temperature by analytical means. We also discuss the link between the reduced flux and the Rosseland opacity.

5.1 Temperature profile in the accretion flow

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

[TABLE]

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

[TABLE]

in spherical symmetry. (In this section we use the symbol “${\rm d}$” instead of “$\partial$” because of time independence.) It is equal to the ratio of the photon mean free path to the “${E_{\textnormal{rad}}}$ scale height” (Paper I). Taking the derivative of Equation (17) yields

[TABLE]

The last term, by Equation (18), is

[TABLE]

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

[TABLE]

since ${E_{\textnormal{rad}}}=a{T_{\textnormal{rad}}}^{4}$. Inserting recursively Equation (20) into Equation (22) is likely not fruitful as it generates derivatives of ${f_{\textnormal{red}}}$ of ever-higher order. However, it will prove instructive to perform this once. In full generality, this yields

[TABLE]

where ${\rm d}^{n}_{x^{n}}f(x)\equiv{\rm d}^{n}f(x)/{\rm d}x^{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),

[TABLE]

one obtains for the derivative term the simple expression

[TABLE]

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

[TABLE]

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

[TABLE]

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

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

[TABLE]

in the limit of a radially constant luminosity (in the sense that $\left|{\rm d}\ln L/{\rm d}\ln r\right|\ll 2$, which is the case here as shown below). The temperature here was written as $T$ since we find that ${T_{\textnormal{gas}}}={T_{\textnormal{rad}}}$; if the gas and radiation are not coupled, $T$ should be taken to refer to ${T_{\textnormal{rad}}}$ only. We used Equation (25) since it makes it explicit that the second term on the righthand 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 since ${f_{\textnormal{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 since then $\lambda$ (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_{\textnormal{red}}}$ to be large, and necessary for the radiation to be sufficiently diffusive ($\lambda$ 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 $\lambda$ goes to zero and thus also the second term in Equation (28 since it is multiplied by $\lambda$). Therefore, even strong variations of ${\kappa_{\textnormal{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 ${\rm d}\ln{\kappa_{\textnormal{R}}}/{\rm d}\ln r$, and if $\lambda$ is not too small, this will be able to slow the decrease of the temperature outwards.

Also, it does not seem necessary that the radiation be free streaming (${f_{\textnormal{red}}}\rightarrow 1$) in order to have $T\propto r^{-1/2}$, which is usually associated with free streaming. If ${\kappa_{\textnormal{R}}}$ is constant with radius and (as perhaps a consequence) ${f_{\textnormal{red}}}$ is sufficiently constant, and $\lambda$ has an intermediate value (e.g., $\lambda\sim 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 constant radially.

The middle panel of Figure 6 shows the terms on the righthand side of Equation (28) and parts thereof. In the free-streaming region at $r=2$–8 ${R_{\textnormal{J}}}$, the first and second derivatives of ${f_{\textnormal{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_{\textnormal{red}}}\approx 1$ (see top panel), so that $\lambda\approx 0$. In the flat temperature part at $r=8$–16 ${R_{\textnormal{J}}}$ (shaded region), the opacity slope term dominates but the term with the derivatives of ${f_{\textnormal{red}}}$ is similar in magnitude, with their signed sum dominating over the $-1/2$ term. The result, along with a lower ${f_{\textnormal{red}}}$ (i.e., higher $\lambda$), is a different (namely, almost zero) temperature slope. At $r>16~{}{R_{\textnormal{J}}}$, the radiation remains somewhat diffusive with ${f_{\textnormal{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 $\approx 3\lambda$ (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 ${\rm d}L/{\rm d}r\approx 0$ since the equation is otherwise exact for a free-fall density profile. The choice of the flux limiter barely makes a difference.

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

[TABLE]

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 grey approximation, and to what extent it depends on the approximations (for instance 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).

5.2 Link between reduced flux and opacity

Figure 6 suggests that there is a relationship between the reduced flux ${f_{\textnormal{red}}}$ and the Rosseland mean opacity ${\kappa_{\textnormal{R}}}$. Defining the logarithmic slope

[TABLE]

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

[TABLE]

which holds anywhere in the flow. This generalises equation (25) in Paper I. We used $\lambda$ from Equation (24) but it is trivial to repeat the derivation for another flux limiter.

Equation (31) explains why ${f_{\textnormal{red}}}$ drops in the dust destruction region in Figure 6 (grey band there): the opacity increases such that ${\kappa_{\textnormal{R}}}\rho r\gtrsim 1$, thereby leading to a drop in ${f_{\textnormal{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 outwards decreasing ${f_{\textnormal{red}}}$. From Equation (31a) this implies that $1/{f_{\textnormal{red}}}$ will be larger, i.e., ${f_{\textnormal{red}}}$ smaller, than predicted by Equation (31b), but the trend is the same. Quantitiatively, this works well in this example but to obtain the exact value of ${f_{\textnormal{red}}}$ one would have to use the Levermore & Pomraning (1981) flux limiter since this was chosen for the simulations. Thus, Equation (31) explains the sudden drop of ${f_{\textnormal{red}}}$ where the opacity suddenly increases to become important in the sense of ${\kappa_{\textnormal{R}}}\rho r\gtrsim 1$.

5.3 Schock temperature

5.3.1 Calculation set-up

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

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

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

[TABLE]

which is an implicit equation for ${T_{\rm shock}}$. The second line is valid when the resulting ${\mathcal{M}}\gtrsim 2.5$ (so that $\eta^{\rm 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_{\rm dnstr}}$ is the sum of the luminosity coming from the deep interior and of the compression luminosity: ${L_{\rm dnstr}}={L_{\textnormal{int}}}+{L_{\rm compr}}$. While ${L_{\rm compr}}$ might depend on ${T_{\rm shock}}$ we effectively absorb this dependency into ${L_{\textnormal{int}}}$, treated as a free parameter. Still, the generalisation ${L_{\rm dnstr}}\rightarrow{L_{\rm dnstr}}({T_{\rm shock}})$ could be readily made. With these assumptions, in Equation (32b) the shock temperature enters on the righthand side only through the opacity.

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

[TABLE]

[TABLE]

with $\ell=1+{L_{\rm dnstr}}/{L_{\rm acc,\,max}}$, without needing to assume that ${L_{\rm dnstr}}/{L_{\rm acc,\,max}}$ is small. We defined ${T_{\rm sh,\,fs}}$ as the “free-streaming shock temperature” given by Equation (6b) but now generalised to include a downstream luminosity ${L_{\rm dnstr}}$ (we consider ${L_{\rm dnstr}}=0$ in this work). Similarly, ${T_{\rm sh,\,diff}}$ is the “diffusive shock temperature”, which obtains when the pre-shock material is diffusive. The maximum accretion luminosity is ${L_{\rm acc,\,max}}=G{M_{\textnormal{p}}}{\dot{M}}/{R_{\textnormal{p}}}$, with the actual luminosity at the shock ${L_{\rm acc}}=\eta^{\rm kin}{L_{\rm acc,\,max}}$ (Paper I). The opacity ${\kappa_{\textnormal{R}}}$ appearing in the expressions for ${T_{\rm sh,\,diff}}$ is either the constant opacity or less trivially the value satisfying Equation (32b), as detailed below. Equation (34c) should replace equation (28b) of Paper I since the prefactor there was inexact and the exponent of the ${R_{\textnormal{p}}}$ factor was missing a minus sign. Equations (34b) and (34c) were written to leading order in the quantity ${\kappa_{\textnormal{R}}}{\rho_{\textnormal{ff}}}{r_{\rm shock}}$, in which case they are also independent of the flux limiter. Note that should ${L_{\rm dnstr}}$ dominate over ${L_{\rm 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$ cm2 g${}^{-1})$ simulation in figure 2 of Paper I, Equation (34) predicts ${T_{\rm sh,\,fs}}\approx 3600$ K, which is close to the actual shock temperature 3500 K. This is satisfying, especially given that there ${\kappa_{\textnormal{R}}}\rho{R_{\textnormal{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 (34b) and (34c), which makes a difference in this transition regime between free streaming and diffusion.

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

[TABLE]

with

[TABLE]

The first and second versions of the righthand 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_{\rm sh,\,fs}}$ and stepping up in temperature to find an intersection of the two sides of Equation (35).

5.3.2 Semi-analytical shock temperature solutions

In Figure 7 we display the shock temperature as a function of planet radius according to Equation (35). We consider ${M_{\textnormal{p}}}=1.3$ and 5 ${M_{\textnormal{J}}}$, take ${\dot{M}}=10^{-4,-3,-2}~{}{M_{\oplus}}\,\mbox{yr}^{-1}$, and vary ${R_{\textnormal{p}}}$ from 1 to about 7 ${R_{\textnormal{J}}}$, truncating however at lower ${R_{\textnormal{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, $\log_{10}\left({L_{\rm dnstr}}/{L_{\odot}}\right)=\log_{10}\left({\dot{M}}/{M_{\oplus}}\,\mbox{yr}^{-1}\right)-1$ for definiteness. For reference, the corresponding pre-shock (free-fall) velocity is indicated on the top axis. The solution is also compared to the free-streaming shock temperature ${T_{\rm sh,\,fs}}$ (Equation (33)).

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}}\,\mbox{yr}^{-1}$, planets of low mass (${M_{\textnormal{p}}}\lesssim 3~{}{M_{\textnormal{J}}}$) with moderate to large radii (${R_{\textnormal{p}}}\approx 3$–10 ${R_{\textnormal{J}}}$) could have their ${T_{\rm shock}}$ increased by up to tens of percent, by a few 100 K. This might be important especially for a non-zero ${L_{\textnormal{int}}}$, leading to a higher shock temperature than expected naively.

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 deviating 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_{\textnormal{R}}}{\rho_{\textnormal{ff}}}{R_{\textnormal{p}}}\sim 1$, as mentioned above, and an example is for $({\dot{M}}=10^{-3}~{}{M_{\oplus}}\,\mbox{yr}^{-1},{M_{\textnormal{p}}}=1.3~{}{M_{\textnormal{J}}},{R_{\textnormal{p}}}=3~{}{R_{\textnormal{J}}})$, labelled “A”, where the opacity of the most refractory dust component leads to a higher shock temperature than ${T_{\rm sh,\,fs}}=1000$ K: ${T_{\rm shock}}\approx 1080$ K for the nominal dust model and ${T_{\rm shock}}\approx 1200$ K for the curve sticking out most in Figure 2. Second of all, the case of non-constant high (${\kappa_{\textnormal{R}}}{\rho_{\textnormal{ff}}}{R_{\textnormal{p}}}\gtrsim 1$) opacities opens up the possibility of multiple solutions.

If ${\kappa_{\textnormal{R}}}\rho{R_{\textnormal{p}}}\ll 1$ at the free-streaming shock temperature ${T_{\rm sh,\,fs}}$, it is necessary and sufficient for the opacity slope

[TABLE]

to satisfy $\alpha_{\kappa}>\alpha_{\kappa}^{\rm 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}^{\rm 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_{\rm sh,\,fs}}$.

If on the other hand ${\kappa_{\textnormal{R}}}\rho{R_{\textnormal{p}}}\gtrsim 1$ at ${T_{\rm sh,\,fs}}$ (i.e., the pre-shock gas would be diffusive already if it were at that temperature), ${T_{\rm sh,\,fs}}$ will not be a solution. Provided the opacity slope $\alpha_{\kappa}<1$ at some higher temperature, there will only be a single888 In the case that the opacity shows several non-monotonicities at higher temperatures the number of roots is of course different and the discussion would need to be adapted. This is presumably the case at the iron opacity “bump” at $\log T\approx 5.3$ (Iglesias et al., 1992; Jiang et al., 2015), and in principle also for a very metal-rich gas. 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_{\rm sh,\,fs}}$.

If however ${\kappa_{\textnormal{R}}}\rho{R_{\textnormal{p}}}\gtrsim 1$ at ${T_{\rm sh,\,fs}}$ and $\alpha_{\kappa}$ remained $>1$ at higher temperatures there would be formally 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 though whether in this case we would still find a Mach number such that $\eta^{\rm kin}=1$, or whether 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 small downstream luminosity (coming from compression, an interior luminosity, or both) does not change significantly 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_{\rm dnstr}}$, the Mach number ${\mathcal{M}}$ increases with increasing ${M_{\textnormal{p}}}$ or $1/{R_{\textnormal{p}}}$ but, perhaps counterintuitively, decreases with increasing accretion rate. Using Equation (6),

[TABLE]

where we took $\eta^{\rm kin}=1$ and $\Delta{f_{\textnormal{red}}}=1$ for the second expression. Thus the Mach number depends only moderately on the mass (${\mathcal{M}}\propto{M_{\textnormal{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 pre-shock Mach number obtained from Equation (32b) with $\mu=1.23$, and also plot Equation (38) for reference. Mach numbers increase with decreasing accretion rate, reaching ${\mathcal{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_{\rm sh,\,fs}}$, given by Equation (38), provides an upper bound but ${\mathcal{M}}$ still remains securely above ${\mathcal{M}}\approx 2.5$. There, $\eta^{\rm kin}$ is near 100 %, which justifies a posteriori our approximation. Also, taking ${L_{\rm dnstr}}$ not large (see Figure 8) but non-zero barely changes these results.

5.3.5 Discussion of the analytical shock temperature

The preceding analysis has shown that, for an isothermal radiative shock with negligible downstream luminosity, the shock temperature is determined by the conditions immediately upstream. There can be small to very important deviations from the free-streaming result ${T_{\rm sh,\,fs}}$. These occur either because of the dust contribution to the opacity (at low temperature) or because of the high opacity and high pre-shock density at high temperature and accretion rates. Nevertheless, the ${T_{\rm sh,\,fs}}$ solution remains valid for a large part of parameter space.

Interestingly, our derivations do not depend explicitly on the equation of state (EOS), specifically on the choice of and (non-)constancy of $\mu$ and $\gamma$ across or ahead of the shock. Thus the expressions (32–34) for the shock temperature could apply also to the case of a general EOS. This is should certainly be the case for $({\dot{M}},{M_{\textnormal{p}}},{R_{\textnormal{p}}})$ combinations for which the pre- and post-shock points have the same $\gamma$ and $\mu$ values. It might also hold more generally but this will be investigated in a subsequent article.

Finally, note that the importance of a shock temperature increased with respect to ${T_{\rm sh,\,fs}}$ for lower planet masses and larger radii than shown here, as might be relevant during the early stages of detachment (runaway accretion), will have to be assessed with separate formation calculations as in Berardo et al. (2017).

6 Importance of the Equation of State and opacity

Next we verify the robustness of our results by varying different parts of the microphysics that go into our simulations. This also gives us occasion to show global profiles of the accretion flow; we had restricted ourselves in figure 2 of Paper I to the vicinity of the shock.

6.1 Dependence of pre-shock-region quantities on the EOS

While the use of an ideal but non-perfect EOS is deferred to a later article, we study in this section the importance of the constant $\mu$ and $\gamma$. For conciseness, we will refer to a particular $(\mu,\gamma)$ combination as an EOS.

Figure 9 shows the profiles for the four extreme combinations $(\mu,\gamma)$, where $\mu=1.23$ (atomic) or $\mu=2.29$ (molecular), and $\gamma=1.1$ (occurring for dissociation or ionization) or $\gamma=5/3$ (monatomic, no dissociation or ionization) together with the nominal case $(\mu=2.29,\gamma=1.44)$, i.e., for a perfect-gas H2–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 $\Delta\rho\propto\mu$. The pressure in the accretion flow does scale as $P\propto 1/\mu$ but, being a dynamical quantity, the ram pressure (i.e., the post-shock pressure) is the same in all cases. While this might be counterintuitive, also the shock temperature (see inset) is independent of $(\mu,\gamma$), which was to be expected from the analytical estimates of the shock temperature since the microphysical parameters do not enter anywhere. The opacities ${\kappa_{\textnormal{R}}}$ and ${\kappa_{\textnormal{P}}}$ are the same here because we use in all cases the same tables (computed using a non-perfect EOS) while the opacity is fundamentally a function of temperature and density (but not mean molecular weight).

6.1.1 Entropy

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

[TABLE]

where $s_{0}$ is a constant, and $T_{0}$ and $P_{0}$ are an arbitrary reference temperature and pressure, respectively. The form is similar as a function of $(P,\rho)$ (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 pre-shock location, since we are interested in the radial profiles of $s$ for the individual $(\mu,\gamma)$ cases. It increases outwards for $\gamma=5/3$ since it is greater than the critical value of 4/3 (section 3.3.2 of Marleau et al., 2017) but decreases for $\gamma=1.1$ (since then $\gamma<4/3$) outwards in the accretion flow. Convection is however not expected there because of the supersonic motion.

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 which 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 $\pm\Lambda$ terms cancel even in non-equilibrium radiation transport; Equation (3)) and writing the gravity field term as a potential:

[TABLE]

where $h=H/\rho={e_{\textnormal{int}}}+P/\rho$ is the specific enthalpy per mass and is given by $h=\gamma/(\gamma-1){k_{\textnormal{B}}}T/(\mu{m_{\textnormal{H}}})$ for a perfect EOS. Only Equation (40a) 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 (40b) is obtained from ${\dot{M}}=4\pi r^{2}\rho v$ (Equation (1a)) and the time-independent momentum equation (Equation (1b)). Equation (40d) follows from the first law of thermodynamics

[TABLE]

with the entropy $s$ here not normalised by ${k_{\textnormal{B}}}\,{m_{\textnormal{H}}}^{-1}$. Note that, in the pre-shock region, it is the small deviation of the velocity from a free-fall profile which makes ${\rm d}L/{\rm d}r\propto T{\rm d}s/{\rm d}r$ and not $\propto{\rm d}h/{\rm d}r$.

In Paper I, we had argued that the radial non-constancy (over $\sim$ Hill-sphere scales) of the luminosity can be understood from the enthalpy profile when the second term in Equation (40a) is negligible. While this is formally true, Equation (40) provides a more complete derivation which 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 (40d) becomes

[TABLE]

where the second line holds for a free-fall density profile. Since the factor $T/r$ in Equation (42b) decreases outwards, it is the layers closest to the shock which contribute most to the change in $L$ between ${r_{\rm shock}}$ and ${r_{\rm 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 $\gamma$ tends to 1 (i.e., more isothermal), or for smaller mean molecular weight $\mu$. The choice of $\gamma$ can change the sign of ${\rm d}L/{\rm d}r$ whereas $\mu$ 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_{\textnormal{red}}})$ temperature profile (see Equation (20)),

[TABLE]

where ${T_{\rm shock}}$ is the shock temperature, which by Equation (35) can be larger than the free-streaming temperature. (Note that $T\propto r^{-1/2}$ can hold not only for free-streaming.) In Figure 6c or 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

[TABLE]

using the subscript “(43)” on $L$ to remind that a temperature profile given by Equation (43) was assumed. Thus, the change (drop or increase) in ${L_{(\ref{Gl:T(r) konst L/fred})}}$ between ${R_{\textnormal{p}}}$ and ${R_{\textnormal{acc}}}$, neglecting terms of order $\sqrt{{R_{\textnormal{p}}}/{R_{\textnormal{acc}}}}\ll 1$, is given by

[TABLE]

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

[TABLE]

Equation (45) is more general than equation (32) in Paper I. In the usual limit ${R_{\textnormal{acc}}}\gg{R_{\textnormal{p}}}$, the luminosity change $\Delta{L_{(\ref{Gl:T(r) konst L/fred})}}$ is thus independent of ${R_{\textnormal{acc}}}$, and thus also of ${k_{\textnormal{Lissauer}}}=1/3$ (as we assume here) for the accretion radius ${R_{\textnormal{acc}}}\approx{k_{\textnormal{Lissauer}}}{R_{\textnormal{Hill}}}$. For $\gamma>4/3$, the luminosity increases outwards ($\Delta{L_{(\ref{Gl:T(r) konst L/fred})}}>0$), whereas for $\gamma<4/3$ there is a decrease in luminosity ($\Delta{L_{(\ref{Gl:T(r) konst L/fred})}}<0$).

Equation (45) is an upper bound in the case $\gamma<4/3$ because any flattening of the luminosity profile as in Figure 9 will slow down the decrease of $T(r)$. For $\gamma>4/3$ the estimate is more like a lower bound but since the critical $\gamma=4/3\approx 1.33$ is not far from $\gamma=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 outwards for the high-$\gamma$ cases ($\gamma=5/3$) as well as for the nominal case ($\gamma=1.44$), while it overall decreases slightly for the low-$\gamma$ cases ($\gamma=1.1$). The mean molecular weight also plays a role, and, in all, $L$ at roughly the Hill radius varies by roughly 10 percent for this specific $({\dot{M}},{M_{\textnormal{p}}},{R_{\textnormal{p}}})$ case over all $(\gamma,\mu)$ combinations. In Section 7.3 we look at this more generally.

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

6.1.3 Global physical efficiency

Finally, depending on both $\gamma$ and $\mu$, the global physical efficiency ranges from $\eta^{\rm phys}=84\,\%$ to 98 %, increasing with $\mu$ but decreasing with $\gamma$. This large difference comes from different $\gamma$ values but also the changing Mach number (through $\mu$), as can be seen from Equation (38). One can wonder whether a low $\mu$ caused by ionization could lead to much lower Mach numbers and thus efficiencies. However, since the Mach numbers we find are all high ($10\lesssim{\mathcal{M}}\propto 1/\sqrt{\mu}$), even a change by a factor of at the very most $\sqrt{0.6/2.3}=0.5$ (from ionized to molecular gas) would leave ${\mathcal{M}}\gtrsim 5>2.5$ and thus the shock supercritical for a given ${T_{\rm shock}}$ and ${v_{\textnormal{ff}}}$. This statement should still be revisited with simulations using the full EOS since ${\mathcal{M}}$ is of course a function of ${T_{\rm shock}}$.

6.1.4 Other quantities

The other panels of Figure 9 are shown for completeness. The profiles are qualitatively similar between all simulations (see the detailed description in Paper I) and we can mention here that we find the gas and radiation are always well coupled, which is presumably related to the sufficiently high opacity. Also, as in Paper I, the radiative precursor to the shock is larger than the simulation box, as even a cursory comparison to standard supercritical and subcritical shock structures (e.g., Ensman, 1994) reveals.

6.1.5 Summary of the effect of the EOS

In summary, depending on the EOS, a different luminosity reaches the accretion radius but the variation is moderate (tens of percent) across the relevant input parameter range. However, both the shock temperature and the ram (post-shock) pressure are independent of the perfect EOS, at least in the limit of an isothermal shock. This is an important result. We conjecture that the post-shock temperature and pressure will not be different when using instead an ideal but non-perfect EOS. This will need to be verified but if if holds, and if the shock is still isothermal, it implies that the post-shock entropy ${s_{\rm post}}$, which depends only on $T$ and $P$, will be the same as found here. As for the shock efficiency, it clearly depends on the EOS (see Figure 9).

6.2 Influence of the dust opacity

Simulations such as the ones 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 yet 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 non-equilibrium effects is beyond the scope of this work. Since 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_{\textnormal{dest}}}=1220\times{\rho_{-11}}^{0.0195}$ K, where $\rho_{-11}\equiv\rho\times 10^{11}~{}\mathrm{cm}^{3}\,\mathrm{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_{\textnormal{J}}}$ cases at ${\dot{M}}=10^{-3}~{}{M_{\oplus}}\,\mbox{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 ($\rho$, $v$, $T$, $L$, ${f_{\textnormal{red}}}$), with the obvious exception of ${\kappa_{\textnormal{R}}}$ and ${\kappa_{\textnormal{P}}}$, the profiles are essentially identical for these two cases. This holds also 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, $\Delta{\tau_{\mathrm{R}}}\sim 10^{-5}$ to $10^{-2}$. Note that the opacity is too high in the always-dust case (see blue solid line in opacity panel) but that this does not affect the temperature structure of higher post-shock regions, i.e., directly below the shock.

As a more extreme case, we also switched off entirely the dust contribution shown by grey profiles in Figure 10. This represents the limiting case of the reduction of Szulágyi et al. (2017) by a factor of ten relative to the ISM, on the grounds that the growth of dust grains into larger aggregates diminishes the opacity. Also, as pointed out by Uyama et al. (2017), dust tends to settle to the midplane while accretion comes from higher up in the circumstellar disc, 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 of the shock temperature by only 100 K, associated with a jump in the reduced flux increasing from $\Delta{f_{\textnormal{red}}}\approx 0.7$ to $\Delta{f_{\textnormal{red}}}\approx 1$.

Our conclusion from this test is that the upstream ${f_{\textnormal{red}}}$ (given ${f_{\textnormal{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, since the outer parts are always cold enough that dust should be present999It is easy to estimate, assuming $T\propto r^{-1/2}$, that only unrealistic parameter combinations could lead to evaporated dust (i.e., $T\gtrsim 1000$ K, taking the density dependence into account) at $r\sim{k_{\textnormal{Lissauer}}}{R_{\textnormal{Hill}}}$ or $r\sim{R_{\textnormal{Bondi}}}$.. However, this will affect only a very small fraction of cases. Indeed, it requires (i) pre-shock temperatures above the (density-dependent) dust destruction temperature but (ii) not too high to not already have the dust destroyed a large distance ahead from shock, in which case the dust would certainly be evaporated. The transition region is very narrow in temperature ($\Delta T\sim 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_{\textnormal{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 it will not be investigated further.

7 Discussion

We discuss some of the results presented here.

7.1 Hot starts or cold starts?

The main question driving this work is whether shock processing of the accreting gas leads to hot starts or to 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 vs. 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 be added to the planet. However, this view neglects the thermal energy of the gas that is brought into the planet, which comes from pre-heating 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^{\rm phys}$ instead of $\eta^{\rm 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 large compared to the internal luminosity.

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

[TABLE]

A large ${q^{+}_{\rm shock,\,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 which is not lost from the system and therefore added to the planet. By the definition of $\eta^{\rm phys}$ (Equation (10)), this is

[TABLE]

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

[TABLE]

using also that for an isothermal shock the post- and pre-shock densities are related by $\rho_{2}=\rho_{1}\gamma{\mathcal{M}}^{2}$ and that $\rho_{2}v_{2}=\rho_{1}v_{1}$. We recall that ${L_{\rm acc,\,max}}=G{M_{\textnormal{p}}}{\dot{M}}/{R_{\textnormal{p}}}$. Equation (49b) 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 post-shock settling velocity of the gas). The latter is vanishingly small in the limit of a large pre-shock Mach number.

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

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

[TABLE]

where again $\ell=1+{L_{\rm dnstr}}/{L_{\rm acc,\,max}}$, the downstream luminosity ${L_{\rm dnstr}}={L_{\textnormal{int}}}+{L_{\rm compr}}$ being the sum of the luminosity coming from the deep interior and the compression luminosity, and writing $\gamma_{1.44}\equiv\gamma/1.44$. For Equations (50b) and (50c), 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), $\ell$ formally depends on $({\dot{M}},{M_{\textnormal{p}}},{R_{\textnormal{p}}})$ but this dependence is negligible in the limit of ${L_{\rm dnstr}}\ll{L_{\rm 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^{+}_{\rm shock,\,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. Since 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$ Myr for 1 ${M_{\textnormal{J}}}$, 1–2 Myr for 2 ${M_{\textnormal{J}}}$, and $\sim 6$ Myr and increasing for 4 ${M_{\textnormal{J}}}$ and up101010 Alternatively, one could look at the $e$-folding times of the luminosity, which are 4 Myr for ${t_{\textnormal{exp}}}=1~{}{M_{\textnormal{J}}}$, 10 Myr for 2 ${M_{\textnormal{J}}}$, 60 Myr for 4 ${M_{\textnormal{J}}}$, and increasing. Surprisingly, the Kelvin–Helmholtz times ${t_{\rm KH}}=G{M_{\textnormal{p}}}^{2}/{R_{\textnormal{p}}}{L_{\textnormal{p}}}=1~{}\mathrm{Gyr}~{}({M_{\textnormal{p}}}/4~{}{M_{\textnormal{J}}})^{2}/[({R_{\textnormal{p}}}/1.3~{}{R_{\textnormal{J}}})({L_{\textnormal{p}}}/3\times 10^{-6}~{}{L_{\odot}})]$ are longer than this $e$-folding time by roughly a large factor of 30 (i.e., ca. 1.5 dex). The reason for this difference is not clear. Note that ${t_{\rm KH}}$ is within ca. $\pm 0.5$ dex of the cooling time $t_{S}={M_{\textnormal{p}}}\mkern 1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5muS/{L_{\textnormal{p}}}$ defined in Marleau & Cumming (2014), where $S$ is the entropy of the planet and $\mkern 1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern 1.5mu$ its mass-averaged temperature. Thus both $t_{S}$ and ${t_{\rm KH}}$ are hardly adequate proxies for ${t_{90}}$ or ${t_{\textnormal{exp}}}$, and one must exercise caution when using ${t_{\rm KH}}$ in timescale-based arguments. . The time spent in runaway accretion in their models is around ${t_{\rm acc}}=0.1$ Myr to 0.5 Myr for 1 to 10 ${M_{\textnormal{J}}}$, respectively. Thus even for ${M_{\textnormal{p}}}=1~{}{M_{\textnormal{J}}}$ the cooling timescale ${t_{90}}$ is much larger than ${t_{\rm 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 11b displays the approximate heating by the shock (Equation (50c)). Considering the extreme combinations of $\mu$ and $\gamma$ as in Section 6, we find ${Q^{+}_{\rm shock}}\approx 10^{-4}$–$10^{-3}~{}{L_{\odot}}$ for ${\dot{M}}\approx 10^{-2}~{}{M_{\oplus}}\,\mathrm{yr}^{-1}$ and ${Q^{+}_{\rm shock}}\approx 10^{-5}$–$10^{-4}~{}{L_{\odot}}$ for $10^{-3}~{}{M_{\oplus}}\,\mathrm{yr}^{-1}$. The radius is fixed at ${R_{\textnormal{p}}}=1.5~{}{R_{\textnormal{J}}}$, a typical value for forming planets (Mordasini et al., 2012a), but varying it in a reasonable range ${R_{\textnormal{p}}}\approx 1$–3 ${R_{\textnormal{J}}}$ does not change the outcome qualitatively.

We compare in Figure 11b 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}}$ towards the end (Hubickyj et al., 2005; Marley et al., 2007). Thus, the relevant heating is around ${Q^{+}_{\rm shock}}>10^{-4}~{}{L_{\odot}}$. This is one to two orders of magnitude higher (!) than the (post-)formation luminosities ${L_{\textnormal{int}}}\approx 10^{-6}~{}{L_{\odot}}$ of Marley et al. (2007). At $1~{}{M_{\textnormal{J}}}$ the heating could be only moderate but for ${M_{\textnormal{p}}}\gtrsim 2~{}{M_{\textnormal{J}}}$ the conclusion becomes more secure. Also, taking ${L_{\rm dnstr}}\neq 0$ into account for the estimate would only lead to a lower Mach number and thus a higher ${Q^{+}_{\rm 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

Berardo & Cumming (2017) and Cumming et al. (2018) followed the time-dependent internal structure of accreting planets with constant accretion rates. They specified their outer boundary conditions for the planet structure as a temperature $T_{0}$ at a pressure equal to the ram pressure, $P_{0}=P_{\rm ram}$. Berardo & Cumming (2017) report that setting $T_{0}$ as a fraction111111These authors write ${T_{0}}^{4}=\chi T_{\mathrm{acc}}^{4}+{T_{\textnormal{eff}}}^{4}$ but their $T_{\mathrm{acc}}$ differs from our equivalent quantity, ${T_{\rm sh,\,fs}}$, with $T_{\mathrm{acc}}=4^{1/4}{T_{\rm sh,\,fs}}\approx 1.4{T_{\rm sh,\,fs}}$ in the case of $\ell=0$ in Equation (33). We therefore use here $f$ and not $\chi$. Note also that $\chi$ in Cumming et al. (2018) was written as $\eta$ in Berardo & Cumming (2017) but that it should not be confused with the efficiency. Our $f$ here can be larger than 1, when the pre-shock ${\kappa_{\textnormal{R}}}\rho{R_{\textnormal{p}}}$ is large (see Section 5.3).

$f$ of the free-streaming temperature ${T_{\rm sh,\,fs}}$ (plus a relative contribution from the internal luminosity) led to fully radiative interiors at the end of formation for $f$ above a certain $f_{\rm min}$. Smaller values of $f$ resulted in convective interiors. The minimum fraction $f_{\rm min}$ was lower for larger accretion rates, with $f_{\rm min}\approx 1$ for ${\dot{M}}=10^{-3}$ and $f_{\rm min}\approx 0.1^{1/4}\approx 0.6$ for ${\dot{M}}=10^{-2}$ ${M_{\oplus}}\,\mathrm{yr}^{-1}$. Since we find $f\geqslant 1$—i.e., the temperature at the ram pressure matches the temperature in the $\eta^{\rm kin}=100\,\%$ limit (Equation (33))—we expect formation calculations using our results to lead to radiative planets. Note that even though we considered only ${L_{\textnormal{int}}}=0$ here, we expect the same result ($\eta^{\rm kin}=1$ and thus $f=1$) to hold for non-zero interior luminosity. It should however still be explored systematically, especially since the result of Berardo & Cumming (2017) was for a specific choice of pre-runaway entropy of the planet $S_{i}$, and it is the contrast between this entropy $S_{i}$ and the immediate post-shock entropy $S_{0}=S(T_{0},P_{0})$ that matters (Berardo et al., 2017).

7.1.3 Entropy-based argument

Since our post-shock entropies are larger or much larger than post-formation entropies (and thus, neglecting cooling, entropies during formation) of Mordasini et al. (2017), we certainly do not expect the shock to be able to cool the planet as it accretes. Berardo et al. (2017) found that they needed extremely low shock entropies (with temperatures of order ${T_{\rm shock}}\approx 150$ K) to reproduce the cold starts of Marley et al. (2007). We however find ${T_{\rm shock}}>1000$ K, down to ${\dot{M}}=10^{-3}~{}{M_{\oplus}}\,\mathrm{yr}^{-1}$. Thus the high post-shock entropy will at least slow down the cooling of the planet during its formation (the “stalling” regime of Berardo et al., 2017), if not heat the planet (“heating regime”), but should not allow for any decrease of the entropy (“cooling regime”). From this argument too, then, cold starts seem unlikely.

7.2 Opacities and gas–radiation coupling

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_{\textnormal{gas}}}\approx{T_{\textnormal{rad}}}$, if the energy exchange time (controlled by the absorption coefficient $\rho{\kappa_{\textnormal{P}}}$; cf. Malygin et al., 2017) is shorter than the time needed for the gas to reach the shock (controlled by ${v_{\textnormal{ff}}}$). Therefore we take a critical look at uncertainties concerning the opacities.

• •

We find that everywhere except in the shock (in the Zel’dovich spike), the pre-shock temperatures of the gas and radiation are equal. Simulations which use non-equilibrium radiation transport with tables with unrealistically low Planck mean values (cf. Figure 2b) might not find that the shock is able to pre-heat the gas. Thus it is important result that the Planck values are high enough for the radiation and gas to be coupled.

• •

A quite uncertain aspect of dust opacity computations is the distribution of dust grains 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.

• •

Nevertheless, the presence of dust in the accretion flow was seen to affect the temperature at and beyond the dust destruction front (cf. 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_{\textnormal{red}}})$ profile, Equation (43) (albeit with the same powerlaw 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 thus given by Equation (43).

• •

We conducted tests as in Paper I with constant low opacity. Typically, as we verified separately (not show), at lower values ${\kappa_{\textnormal{R}}}={\kappa_{\textnormal{P}}}\sim 0.01$ cm2 g*-1* the radiation and gas temperatures stay decoupled even in the high-density post-shock region. This is however entirely unrealistic given that at those densities ($\rho\approx 10^{-10}$–$10^{-8}~{}\textrm{g\,cm}^{-3}$) and temperatures (${T_{\textnormal{gas}}}\approx 500$–5000 K; cf. figure 4 of Paper I) the Planck opacity is rather of order ${\kappa_{\textnormal{P}}}\sim 10$ cm2 g*-1* (Figure 2b).

• •

Finally, we related the behaviour 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, however keeping 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 non-equilibrium (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 disc in Aoyama et al. (2018) and for the accretion shock on the planet surface in a forthcoming publication (Aoyama et al., in prep.).

7.3 Equation of State

In these first two papers (Marleau et al. 2017, this work) we have restricted ourselves to a perfect equation of state (EOS; constant $\mu$ and $\gamma$). 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 post-shock compression luminosity, and thus also (iii) the more precise value of ${T_{\rm shock}}$ (through the $\Delta{f_{\textnormal{red}}}$ factor) should all be affected to some extent by the EOS, at least for some combinations of $({\dot{M}},{M_{\textnormal{p}}},{R_{\textnormal{p}}})$.

Concerning item (i), we studied in Section 6.1.2 an estimate $\Delta{L_{(\ref{Gl:T(r) konst L/fred})}}$ 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 $|\Delta{L_{(\ref{Gl:T(r) konst L/fred})}}|/{L_{\rm 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_{\textnormal{p}}}=1$–10 ${M_{\textnormal{J}}}$, and ${R_{\textnormal{p}}}=1$–5 ${R_{\textnormal{J}}}$ (an even wider parameter space than what we consider for our simulations), the relative drop or increase is never larger than $|\Delta{L_{(\ref{Gl:T(r) konst L/fred})}}|/{L_{\rm 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\propto 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\propto r^{-1/2}$ profile) is less than 2 % for any $({\dot{M}},{M_{\textnormal{p}}},{R_{\textnormal{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 $\Delta{L_{(\ref{Gl:T(r) konst L/fred})}}$ justifies a posteriori the assumption of constant $(L/{f_{\textnormal{red}}})$ made to derive it. Indeed, a relative change $|\Delta{L_{(\ref{Gl:T(r) konst L/fred})}}|/{L_{\rm acc}}\approx 10$–15 % over 2–3 dex in radius (from ${R_{\textnormal{p}}}\approx 1$–3 ${R_{\textnormal{J}}}$ to ${R_{\textnormal{acc}}}\approx{k_{\textnormal{Lissauer}}}{R_{\textnormal{Hill}}}({M_{\textnormal{p}}})\approx 250$–500 ${R_{\textnormal{J}}}$ for ${M_{\textnormal{p}}}\approx 1$–10 ${M_{\textnormal{J}}}$ at 5.2 au) corresponds to an approximate average slope (see Equation (30)) of at most $|\beta|\approx\log_{10}(1.15)/2=0.04$ if ${f_{\textnormal{red}}}$ is constant or $|\beta|\approx 0.4$ at most if ${f_{\textnormal{red}}}$ also changes by a factor of ca. 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_{\textnormal{red}}}$ used to derive the change in $L$.

8 Summary and conclusions

In this series of papers (Marleau et al., 2017; this work; Marleau et al., in prep.) we take a detailed look at the physics of the accretion shock in planet formation. In this second paper, we have updated to disequilibrium (2-$T$) radiation transport (i.e., following the gas and radiation energy densities separately) and modern opacities, especially for the gas: we use the dust opacities of Semenov et al. (2003) but the gas opacities of Malygin et al. (2014), avoiding the too-low Planck mean opacities normally included in Semenov et al. (2003).We have also now surveyed a range of values for the formation parameters $({\dot{M}},{M_{\textnormal{p}}},{R_{\textnormal{p}}})$, assuming negligible ${L_{\textnormal{int}}}$ (see Figure 1), namely ${\dot{M}}=10^{-3}$–$10^{-2}~{}{M_{\oplus}}\,\mathrm{yr}^{-1}$, ${M_{\textnormal{p}}}\approx 1$–10 ${M_{\textnormal{J}}}$, ${R_{\textnormal{p}}}\approx 1.6$–3 ${R_{\textnormal{J}}}$. This has motivated us to several semi-analytical derivations along with comparisons to simulation outputs. We have kept the simplification of a perfect equation of state and focused on the case of molecular hydrogen with a cosmic admixture of helium ($\mu=2.29$, $\gamma=1.44$).

We now summarize our primary findings on different aspects. Concerning the thermal and radiative properties of the accretion flow:

Both our simulations and analytical theory show that the behaviour of the luminosity in the accretion flow is not the direct result of radiative transfer effects but rather depends on the equation of state (Section 6.1.2). The luminosity turns out to be radially constant to $\approx 2$ % for values of $\mu$ and $\gamma$ appropriate for H2 + He. Taking other values of $\mu$ and $\gamma$ increases or decreases the change in $L$ between the shock and roughly the Hill sphere. However, the maximum change is relatively small with $|\Delta{L_{(\ref{Gl:T(r) konst L/fred})}}|/{L_{\rm acc}}\lesssim 15$ % across the relevant parameter space and for any ($\gamma\geqslant 1.1,\mu\geqslant 1.23$) combination (Section 7.3). We highlight that the $\Delta{\tau_{\mathrm{R}}}\sim 1$ surface is not of any particular significance (Section 4.3). 2. 2.

Thanks to the sufficiently high Planck mean opacities (for which a contribution from the dust is not needed), the matter and radiation are very well coupled both ahead of and behind the shock, i.e., everywhere except in the Zel’dovich spike (Figures 2 and 9). In fact, non-equilibrium (2-$T$) radiation transport could be neglected when studying only the post-shock temperature and pressure or the global energetics. 3. 3.

As found in Paper I and confirmed here, the radiative precursor to the shock (Mihalas & Mihalas, 1984; Commerçon et al., 2011a) is larger than the simulation domain, which is roughly the Hill sphere, even in the case of somewhat high Rosseland optical depth ($\Delta{\tau_{\mathrm{R}}}\sim 10$). 4. 4.

The pre-shock region close to the planet, out to some Rosseland optical depth, is usually in the free-streaming regime and not, as one would expect for a supercritical shock, in the diffusion limit (e.g., figure 8 of Vaytet et al. 2013b). Thus the shock is a thick–thin shock in the classification of Drake (2006). At low shock temperatures ($T\lesssim 1500$ K) the dust is still present, making the pre-shock region somewhat diffusive and raising the shock temperature.

The shock properties were a focus of this study and we found the following:

As in Paper I, all shocks are isothermal and supercritical, and the Mach numbers are high enough for $\eta^{\rm kin}\approx 100\,\%$ of the incoming kinetic energy flux to be converted to radiation locally at the shock (see grey lines in Figure 4 and Figure 8). The post-shock pressure is equal to the ram pressure $P_{\rm ram}$. 2. 2.

The free-streaming analytical estimates of the shock temperature (Equation 6b) and of the upstream luminosity (Equation 8) were seen to hold very well over a large portion of parameter space. Importantly, we found out that this holds also for high optical depths between the shock and the nebula. Deviations of $\sim 5\,\%$ in $T$ occur at low shock temperatures (Figure 3a). 3. 3.

An important analytical development was the derivation of an implicit equation (Equation (35)) for the shock temperature ${T_{\rm shock}}$ given a Rosseland mean opacity function ${\kappa_{\textnormal{R}}}(\rho,{T_{\rm shock}})$. We solved this numerically (Figure 7). 4. 4.

Based on our analysis, ${T_{\rm shock}}$ should not be affected to first order by the use of a non-perfect ideal EOS (i.e., considering dissociation and ionization), since $\gamma$ and $\mu$ do not enter in the derivation of ${T_{\rm shock}}$ (see Equation (32)). However, (i) the luminosity in the accretion flow (and thus at the Hill sphere), (ii) the post-shock compression luminosity, and thus also (iii) the more precise value of ${T_{\rm shock}}$ (through the $\Delta{f_{\textnormal{red}}}$ factor; Equation (6)) should all be depend somewhat on the EOS. This will be assessed in a subsequent paper (Marleau et al., in prep.). 5. 5.

We calculated the post-shock entropies immediately below the shock using an EOS taking dissociation and ionisation into account (appendix A of Berardo et al., 2017). While this is formally not consistent with our (perfect-EOS) simulations, we argued this is likely accurate since ${T_{\rm shock}}$ and ${P_{\textnormal{post}}}=P_{\rm ram}$ are probably independent of the EOS. The immediate post-shock entropies were found to be between approximately 13 and 20 ${k_{\textnormal{B}}}\,{m_{\textnormal{H}}}^{-1}$ for our range of parameters (Section 4.6, Figure 5). These values are high compared to the post-formation entropy of planets, which is at most around 10–14 ${k_{\textnormal{B}}}\,{m_{\textnormal{H}}}^{-1}$ according to current, though not definitive, predictions (Berardo & Cumming, 2017; Berardo et al., 2017; Mordasini, 2013; Mordasini et al., 2017). However, we caution and emphasize that the actual entropy of the gas added to the planet, below the post-shock settling layer, is different from this immediate post-shock entropy. This is explored in Berardo et al. (2017) and Berardo & Cumming (2017).

Finally, a key output of our simulations was the efficiency of the shock:

We have measured the physical efficiency $\eta^{\rm phys}$ as a function of accretion rate ${\dot{M}}$, planet mass ${M_{\textnormal{p}}}$, and planet radius ${R_{\textnormal{p}}}$ (Figure 5) and derived it semi-analytically (Equations (13) and (38)). This efficiency captures the global energy recycling occurring in the pre-shock region (Paper I). We saw (Figure 5) that the efficiencies are always greater than roughly 97 % for the range of parameters considered here. 2. 2.

Naively, a high $\eta^{\rm phys}$ could suggest that the gas is added “cold” but the part not escaping (i.e., the heating of the planet heating by the shock) turns out to be much larger than the internal luminosity in the Marley et al. (2007) extreme cold starts (Section 7.1.1).

The semi-analytical work presented here revealed the reduced flux ${f_{\textnormal{red}}}={F_{\rm rad}}/(c{E_{\textnormal{rad}}})$ (Equation (7)) to be a powerful quantity for understanding the behaviour of the radiation field (free streaming or diffusing, often termed approximately “optically thin” and “optically thick”). This holds at least for the grey treatment of radiation transport used in this work in a spherically symmetric geometry. When $L/{f_{\textnormal{red}}}$ is sufficiently constant radially, we showed in Equation (31) that there is a simple relation between the reduced flux ${f_{\textnormal{red}}}$ and the opacity, which provides an intuitive understanding of ${f_{\textnormal{red}}}$.

The main results of our simulations are post-shock $(P,T)$ values and global efficiencies $\eta^{\rm phys}$. This is useful respectively for detailed modeling of the structure of accreting planets as in Berardo et al. (2017), Berardo & Cumming (2017), and Cumming et al. (2018), as well as for the one-zone, global approach of, e.g., Hartmann et al. (1997). The Bern model (Alibert et al., 2005; Mordasini et al., 2012b, 2017) is currently in between, calculating detailed planet structures but with the assumption of a radially constant luminosity. Note that the modelling of the energy transfer at the accretion shock is also relevant in the context of star formation (e.g., Baraffe et al., 2012; Geroux et al., 2016; Baraffe et al., 2017; Jensen & Haugbølle, 2018). Researchers interested in using our simulation results can take the semi-analytical formulæ presented above, including the opacity effects for the temperature, under the assumption of a perfect EOS.

The other main outcome is the amount of radiation reaching the Hill sphere. This should be useful input for studies of the thermo-chemical feedback of planets on the local protoplanetary disc, for instance as in a number of recent papers (Cleeves et al., 2015; Cridland et al., 2017; Stamatellos & Inutsuka, 2018; Rab et al., 2019). Within the simplification of a spherical accretion geometry, our results show that essentially all of the accretion shock luminosity is expected to reach the local nebula, and that a high Rosseland optical depth, at least up to $\Delta{\tau_{\mathrm{R}}}\sim 10$, does not lead to significant extinction of the bolometric shock luminosity in the accretion flow.

Finally, we have explored by different means whether our results point towards hot starts or cold starts (Section 7.1). As discussed above, the heating of the planet by the shock ${Q^{+}_{\rm shock}}$ (i.e., the flow rate of inward-going energy; Equation 48) was estimated to be much larger than the internal luminosity for the Marley et al. (2007) classical cold starts. This suggests that they are not entirely realistic. As for the “nominal cold start” or the “hot start” assumption during accretion, Mordasini et al. (2017) showed that both lead to warm or even hot starts. Taken together, all of this might explain why direct imaging observations, with the sole exception of 51 Eri b (Macintosh et al., 2015; Nielsen et al., 2019), have not been finding evidence for planets even consistent with cold starts.

We acknowledge the valuable support of Th. Henning and W. Benz for this project. It is a pleasure to thank A. Mignone and collaborators for their excellent open-source code PLUTO. We also thank W. Kley and H. Klahr for several helpful discussions; N. Malygin for alacritous and competent help as well as for providing routines to read and use his opacity tables in Makemake; D. Semenov for publicly-available fortran code for his opacity data; N. Turner, A. Cumming, C. P. Dullemond, R. P. Nelson, T. Tanigawa, M. Ikoma, D. N. C. Lin, J. Bouwman, and V. Elbakyan for interesting questions, comments, and discussions; and A. Emsenhuber (Bern), V. Lutz and J. Krüger (Tübingen), and U. Hiller (MPIA) for rapid and patient help with the respective computing clusters. We thank the referee for a useful report which lead to several clarifications and motivated a deeper analysis of some aspects. G-DM and RK acknowledge the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (KU 2849/7-1). G-DM and CM acknowledge support from the Swiss National Science Foundation under grant BSSGI0_155816 “PlanetsInTime”. Parts of this work have been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. RK acknowledges financial support within the Emmy Noether research group on “Accretion Flows and Feedback in Realistic Models of Massive Star Formation” funded by the German Research Foundation (DFG) under grants no. KU 2849/3-1 and KU 2849/3-2. Part of the simulations presented here were performed on the ba(t)chelor cluster at the MPIA. The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the Universität Tübingen, the state of Baden-Württemberg through bwHPC, and the German Research Foundation (DFG) through grant no. INST 37/935-1 FUGG. All figures were produced using gnuplot’s terminal epslatex with the font package fouriernc.