Open Access
Issue
A&A
Volume 667, November 2022
Article Number A95
Number of page(s) 19
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202244092
Published online 11 November 2022

© N. Oberg et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1 Introduction

A general feature of regular satellite formation theory is that the circumplanetary disk (CPD) consists of circumstellar material accreted from within the vicinity of the planet (Lubow et al. 1999; Canup & Ward 2002; Alibert et al. 2005; Shibaike et al. 2019; Ronnet & Johansen 2020). If the planet is massive enough to open a gap in the circumstellar disk, material continues to flow into the gap (Kley 1999; Teague et al. 2019) and falls nearly vertically onto the CPD (Tanigawa et al. 2012; Morbidelli et al. 2014). The CPD achieves a steady-state mass when this inflow is balanced by outflow where gas either spirals into the planet or is decreted beyond the Hill sphere (Canup & Ward 2002; Batygin & Morbidelli 2020). Independently of disk viscosity, stellar tides induce spiral waves in the CPD which transport angular momentum and promote accretion onto the planet at a rate on the order 10−7 MJ yr−1 (Rivier et al. 2012), suggesting that for a CPD of mass ~10−4 MJ (D’Angelo et al. 2002; Gressel et al. 2013; Szulágyi et al. 2014) infalling gas spends only a limited time inside the CPD before being lost (Canup & Ward 2002).

The timescale of radial dust drift in small disks is also predicted to be short (Pinilla et al. 2013; Shibaike et al. 2017; Rab et al. 2019). A CPD could lose mm-size grains within 102–103 yr to aerodynamic drag against highly sub-Keplerian gas due to a very steep radial pressure gradient (Zhu et al. 2018). The CPD is thus a very dynamical system, potentially with both inwards and outwards radial transport of both gas and dust on very short timescales.

The amount of time which gas and dust spend within the CPD becomes highly relevant if a chemical “reset” occurs. The infalling circumstellar material may pass through one or more accretion shocks (Lubow et al. 1999; Tanigawa et al. 2012; Zhu 2015; Szulágyi et al. 2016; Schulik et al. 2020) and can be heated ≥1000 K during accretion onto the CPD (Szulágyi 2017; Szulágyi & Mordasini 2017; Aoyama et al. 2018). For pre-shock velocities in excess of 90 km s−1 the shock can be sufficiently hot to leave most of the infalling gas atomic or ionized (Aoyama et al. 2018). The shock may also desorb the icy mantles of grains via sputtering and thermal desorption, which for pre-shock velocities in excess of 8–10 km s−1 can effectively strip a grain of H2O ice (Woitke et al. 1993; Aota et al. 2015; Tielens 2021). Small icy grains passing through the gap may also lose their volatile contents to photodesorption prior to shock-heating (Turner et al. 2012). We refer to a “reset” scenario if shock-heating or photodesorption effectively destroys molecules in the accretion flow to the CPD. In a reset scenario, the re-formation of ices in the CPD must compete with viscous accretion and decretion of gas and radial drift of dust. Alternatively, if circumstellar disk ices survive incorporation into the CPD, we refer to an “inheritance” scenario.

The Galilean satellites characteristically exhibit a radial compositional gradient of decreasing density with increasing distance from Jupiter. The inner moon Io is ice-free, while Europa has an ice mass fraction of ~6–9% (Schubert et al. 2004; Kuskov & Kronrod 2005) and Ganymede and Callisto have ice mass fractions of 40–55% (McKinnon 1997; Schubert et al. 2004; Sohl et al. 2002). Theoretically it appears challenging to reproduce the compositional gradient by tidal heating (Bierson & Steinbrügge 2021) or impact-driven escape of volatiles (Dwyer et al. 2013). Previously it has been proposed that the gradient was imprinted during the formation of the moons by a radial temperature gradient in the CPD (Lunine & Stevenson 1982), but the relevant chemical timescales have rarely been taken into account (Mousis & Alibert 2006). It is an open question whether the gradient can be produced primordially if the chemistry of infalling gas and dust is reset. By analyzing the composition and abundance of ices that are able to form within the relevant timescales we can place a lower bound on how efficiently angular momentum is transported within the CPD.

In this work, we investigated the balance of the competing timescales of ice formation, dust grain drift, and gas viscous flow to seek constraints on properties of the CPD such as viscosity, ice mass fraction, and dust-to-gas ratio. We considered the two opposing extreme cases of full chemical inheritance and reset in chemically evolving disk models utilizing a rate-based modeling approach. In Sect. 2, we describe our modeling setup and the assumptions that we make. In Sect. 3, we analyze the CPD time-dependent ice abundances for both the reset and inheritance cases, and place novel constraints on the properties of CPDs. In Sect. 4, we discuss the implications of these constraints and place them in the context of solar system formation. We consider also the role that radial grain drift has in competing with ice adsorption and desorption. In Sect. 5, we summarize our key findings.

2 Methods

We considered two opposing scenarios; one in which the molecular gas-phase and ice chemistry of the circumstellar disk is preserved during accretion onto the CPD (inherit) and one in which it is lost (reset). In the former case, the CPD is initially populated with gas and ices extracted from the circumstellar disk. In the latter case the disk is initially populated by a fully atomic gas and dust is free of ice. We followed the build-up of ices in the CPD over time in a thermochemical disk model (see Sect. 2.1) and extracted the molecular ice abundance and composition as a function of time.

The net inflow and outflow rate of gas and solids from and to the CPD is assumed to be zero, and that the disk mass is in steady-state. The rate of gas outflow then sets an upper limit on the applicable chemical evolutionary timescales. Hereafter we refer to this as the viscous timescale, defined as tvisc=MCPDM˙,${t_{{\rm{visc}}}} = {{{M_{{\rm{CPD}}}}} \over {\dot M}},$(1)

where MCPD is the mass of the CPD and M is the infall rate of circumstellar material onto the CPD. The relatively short ívisc considered in this work (see Sect. 2.1.2) may cause reactions with high activation energies to be kinetically inhibited, although it has been noted that the relatively high densities characteristic of CPDs may allow these reactions to proceed to equilibrium (de Pater & Lissauer 2010). Nevertheless we contrast the time-dependent results with the assumption of steady-state chemistry. In steady-state chemistry the disk is allowed to evolve for an indefinite time period until the rate of formation for every gas and ice species is balanced by the corresponding rate of destruction.

2.1 Disk modeling code

We used the radiation thermochemical disk modeling code PRODIMO1 (Protoplanetary Disk Model; Woitke et al. 2009, 2016; Kamp et al. 2010, 2017; Thi et al. 2011, 2020a) to explore the formation rates and resulting abundances of various ices in CPDs. ProDiMo uses a rate equation based approach to compute the gas chemistry using either a time-dependent or steady-state solver. The model represents a 2D slice through an axisymmetric disk, extending radially in distance from the planet r and vertically in distance from the disk midplane z. Our chemical network contains 13 elements and 235 atomic and molecular species. Where not explicitly specified we used the “large DIANA chemical standard” network as described in Kamp et al. (2017). Reaction rates are mainly selected from the UMIST2012 database (McElroy et al. 2013). Important three-body collider reactions are adopted from the UMIST2006 rate file (Woodall et al. 2007). Gas-phase reactions within the CPD produce molecular species which can then freeze-out on grain surfaces. The rate of ice formation is determined by the available grain surface area, dust temperature, and the rates of thermal, photo-, and cosmic-ray desorption (see Sect. 2.1.1 for a detailed description of ice formation).

We made the simplifying assumption that the chemical composition of the infalling material is distributed instantaneously and homogeneously throughout the disk (see Appendix E for a discussion of the potential impact of the rate of vertical gas mixing). We assumed that the CPD inherits micrometer to mm-sized dust grains directly from the circumstellar disk. The vertical dust stratification is calculated according to the method of Dubrulle et al. (1995) and kept fixed prior to the calculation of chemical evolution. The timescales of radial dust drift are calculated in a post-processing step described in Sect. 2.2. The temperature and radiation structure of the CPD is solved in steady-state and then kept fixed during chemical evolution of inherited or reset infalling material.

In Sect. 3.2, we considered the implications of including grain-surface chemistry reactions. With surface chemistry ProDiMo models explicitly the formation of H2 in the CPD for which the inclusion of additional chemical species such as hydrogenated PAH is necessitated (Thi et al. 2020b,a). The selection of the additional species and reactions in the surface chemistry network and their role in the eventual composition of the ices will be discussed in an accompanying work focused on the composition of the ices (Oberg et al., in prep.).

2.1.1 Ice formation

Where conditions in the disk are appropriate, ices can condense onto the grains in successive layers by physical van der Waals bonding (physisorption). The adsorption rate of species i onto a physisorption sites is Riads=4πa2vithndSis1,$R_i^{{\rm{ads}}} = 4\pi {a^2}v_i^{{\rm{th}}}{n_{\rm{d}}}{S_i}\,{{\rm{s}}^{ - 1}},$(2)

where 4πa2 is a dust grain surface area, a is the grain radius, vith${\rm{v}}_i^{{\rm{th}}}$ is the thermal speed (kBTg/2πmi)1/2, kB is the Boltzmann constant, Tg is the gas temperature, mi is the mass of the gas-phase species, nd is the dust number density, and Si is the sticking coefficient. The ice adsorption rate coefficients are then dn#,idt=Riadsnicm3s1,${{{\rm{d}}{n_{\# ,i}}} \over {{\rm{d}}t}} = R_i^{{\rm{ads}}}{n_i}\,{\rm{c}}{{\rm{m}}^{ - 3}}{{\rm{s}}^{ - 1}},$(3)

where ni is the number density of the gas-phase species.

Physisorbed species can desorb thermally, photodesorb, or desorb after a cosmic ray impact deposits energy in the grain. For physisorption there is no desorption barrier and the desorption energy is equal to the binding energy Eib$E_i^{\rm{b}}$. The Arrhenius formulation for the rate of thermal desorption is then Rides,th=v0,ieEib/kBTd in s1,$R_i^{{\rm{des,th}}} = {v_{0,i}}{{\rm{e}}^{ - E_i^{\rm{b}}/{k_{\rm{B}}}{T_{\rm{d}}}}}\,{\rm{in}}\,{{\rm{s}}^{ - 1}},$(4)

where thepre-exponential (frequency) factor v0,i is v0,i=2NsurfEibπ2mi.${v_{0,i}} = \sqrt {{{2{N_{{\rm{surf}}}}E_i^{\rm{b}}} \over {{\pi ^2}{m_i}}}} .$(5)

Nsurf is the density of surface binding sites 1.5 × 10 cm−2 (Hasegawa et al. 1992), and Td is the dust temperature. Adsorption energies of major volatile species are listed in Appendix B.

The photodesorption rate is computed using the UV field strength calculated by the radiative transfer and photodissociation cross-sections. The photodesorption rate for species is Rides,ph=πa2ndniactYiχFDraines1,$R_i^{{\rm{des,ph}}} = \pi {a^2}{{{n_{\rm{d}}}} \over {n_i^{{\rm{act}}}}}{Y_i}\chi {F_{{\rm{Draine}}}}{{\rm{s}}^{ - 1}},$(6)

where Yi is the photodesorption yield, niact$n_i^{{\rm{act}}}$ is the concentration of the species in the active layers, niact=niif  n#,tot4πNsurfa2nd         =ni(Nact/Nlayer)if  n#,tot>4πNsurfa2nd.$\matrix{{n_i^{{\rm{act}}} = {n_i}} \hfill & {{\rm{if}}\,\,{n_{\# ,tot}}4\pi {N_{{\rm{surf}}}}{a^2}{n_{\rm{d}}}} \hfill \cr {\,\,\,\,\,\,\,\,\, = {n_i}\left( {{N_{{\rm{act}}}}/{N_{{\rm{layer}}}}} \right)} \hfill & {{\rm{if}}\,\,{n_{\# ,tot}} > 4\pi {N_{{\rm{surf}}}}{a^2}{n_{\rm{d}}}.} \hfill \cr} $(7)

The number of physisorbed layers Nlayer = n#,tot /(4πNsurfa2nd) and n#,tot = ∑in#,i∙ is the total number of the physisorbed species, 4πNsurfa2 is the number of binding sites per layer, and Nact is the number of chemically active layers. χFDraine is a measure of the local UV energy density from the 2D continuum radiative transfer (Woitke et al. 2009). The rate of cosmic-ray induced desorption Rides,CR$R_i^{{\rm{des,CR}}}$ is calculated according to the method of Hasegawa & Herbst (1993). The total desorption rate Rides$R_i^{{\rm{des}}}$ is the sum of the thermal, photo- and cosmic ray induced desorption rates Rides,th+Rides,ph+Rides,CR$R_i^{{\rm{des,th}}} + R_i^{{\rm{des,ph}}} + R_i^{{\rm{des,CR}}}$. The desorption rate of the physisorbed species i is then dnidt=Ridesniactcm3s1,${{{\rm{d}}{n_i}} \over {{\rm{d}}t}} = R_i^{{\rm{des}}}n_i^{{\rm{act}}}{\rm{c}}{{\rm{m}}^{ - 3}}{{\rm{s}}^{ - 1}},$(8)

2.1.2 Properties of the disk models

The circumstellar disk. To generate the chemical abundances for our “inheritance” scenario we considered the properties of the surrounding circumstellar disk from which material is accreted onto the CPD. We used a two-step approach to model the chemistry in our circumstellar disk. As a first step, the initial conditions are derived from a zero-dimensional “molecular cloud” model, the parameters of which are listed in Table G.1. This stage represents 1.7 × 105 yr (the estimated lifetime of the Taurus Molecular Cloud TMC-1) of chemical evolution in a pre-collapse molecular cloud state (McElroy et al. 2013). The resulting chemical abundances of the majority of the most common species agree within a factor 10 with the observed abundances in TMC-1(see Fig. G.1). These abundances are then used as initial conditions for the 2D grid of cells in the circumstellar disk model in the second step.

In the second step, the circumstellar disk model is evolved for an additional 4 Myr to be consistent with the formation timeline of Jupiter proposed to account for the distinct isotopic population of meteorites found in the solar system wherein Jupiter undergoes runaway accretion >3.46Myr after the formation of calcium-aluminium rich refractory inclusions (Kruijer et al. 2017; Weiss & Bottke 2021). The surface density power law and physical extent of the circumstellar disk is based on a modified “Minimum Mass Solar Nebula” (MMSN; Hayashi 1981). A parameteric gap has been introduced which reduces the dust and gas density at 5.2 au centered on the location of Jupiter. The gap dimensions are parameterized by an analytical gap scaling relation derived from hydrodynamical simulations and are consistent with a circumstellar disk viscosity of α ~ 10−4 and disk age of 4 Myr (Kanagawa et al. 2016). A detailed description of the circumstellar disk model and gap structure methodology can be found in Oberg et al. (2020). The relevant circumstellar disk model parameters can be found in Table 1. The dust temperature, surface density profile, and UV field strength in and around the circumstellar disk gap can be seen in Fig. 1.

Finally we extract the chemical abundances from the circumstellar disk model at the outer edge of the gap at a radius of 8.2 au. The gap edge is defined as the radius at which the perturbed surface density profile reaches 50% of the unperturbed profile. As material flows into the gap from above one pressure scale height (Morbidelli et al. 2014), we extract our initial conditions for the CPD model at z = H (z = 0.5 au at r = 8.2 au). The ambient conditions at this point are listed in Table 2. The extracted abundances are then used as the initial chemical composition for our “inheritance” scenario CPD.

Throughout this work we quantify the iciness of solids with the ice mass fraction, fice=micemrock+mice,${f_{{\rm{ice}}}} = {{{m_{{\rm{ice}}}}} \over {{m_{{\rm{rock}}}} + {m_{{\rm{ice}}}}}},$(9)

where mice is the ice mass and mrock is the total rock (in this case, dust) mass. The fice of solids in the inherited circumstellar material is 0.48. At a single pressure scale height (z = H) there is a factor ~20 reduction of the initial, global dust-to-gas ratio due to settling which is calculated according to the method of Dubrulle et al. (1995) with αsetue = 10−2 such that the dust-togas ratio d/g at one scale height d/gz=H = 10−3.2. Nevertheless we test a range of different d/g values for the CPDs both above and below this value (10−4–10−2).

Survival of icy grains passing through the gap. Icy grains must orbit within the optically thin gap for an unknown amount of time prior to being accreted onto the CPD. We considered whether ices on grains can survive against thermal- or photodes-orption while crossing through the circumstellar disk gap. The dust and gas temperatures in the gap are not closely coupled as a consequence of the low densities. While the gas temperature extracted from the model at the midplane is ~200 K, the corresponding dust temperature is 48 ± 2 K. Given that pressures in the gap range from 10−12–10−10 bar, water ice is stable on the relevant timescales in the absence of irradiation (Lodders 2003). However, the actual ice abundance at the gap midplane is negligible in steady-state. Despite the presence of a shadowing inner disk, ice within the gap are sublimated as a result of the significant stellar background radiation scattered into the gap.

To assess the longevity of ices crossing through the gap we populated the gap region with the “inheritance” chemical abundances found exterior to the gap and produce snapshots at regular intervals. The resulting decline in ice abundance as a function of time is shown in Fig. 2 for various stellar luminosities. The differing luminosities correspond to expected properties of the sun for a stellar age of 0.1, 0.5, 1 and 4Myr (Siess et al. 2000). Turner et al. (2012) suggests that grains entering the gap are accreted generally within a single orbital period or 10 yr. We find that for a moon formation time >1 Myr after CAI formation (L* ≤ 2.34 L), grains retain >99% of their volatile content during gap-crossing if they reach the vicinity of the planet within 10 yr. A “full” inheritance scenario is thus not excluded by conditions within the gap, but would instead rely on shock-heating at the CPD surface.

The circumplanetary disks. The CPD is an actively-fed accretion disk with a steady-state mass proportional to its viscosity and mass infall rate. We considered an optically thick CPD of mass 10−7 Mo as well as a lower mass CPD (10−8 M) which is optically thin everywhere outside the orbit of Callisto (r > 0.03 RH where RH = 0.34 au is the Hill radius). For a Jovian-mass planet these represent planet-disk mass ratios of ~10−4 and 10−5, respectively. The CPDs are thus of the “gas-starved” type, and do not instantaneously contain the mass required to form a moon system as massive as the Galilean one.

The outer radius of the CPD is set as the planetary Hill radius RH, however an exponential decline in the surface density profile is parameterized to begin at RH/3 corresponding to the outermost stable orbit set by tidal interaction and angular momentum considerations (tidal truncation radius; Quillen & Trilling 1998; Ayliffe & Bate 2009; Martin & Lubow 2011). An empty inner magnetospheric cavity is assumed to exist as the result of magnetic interaction with the planet (Takata & Stevenson 1996; Shibaike et al. 2019; Batygin 2018). The parameterized radial surface density profiles of the high- and low-mass CPDs can be found in Fig. 3 together with the resulting optical extinction profiles derived from the continuum radiative transfer.

The small parameter variation grid of models exploring possible CPD mass and viscosity can be found in Table 3 along with the model id’s. The format of the id is (xy) where x is related to the CPD mass by MCPD = 10x M and y to the mass infall rate (and thus viscosity) by MCPD = 10−y M yr−1. The list of parameters which are common between the reference CPDs can be found in Table 4.

CPD viscosity. While the mechanism that produces angular momentum transport in accretion disks is not well understood, it is known that molecular viscosity alone is far too weak to explain observations (Shakura & Sunyaev 1973; Pringle 1981). The efficiency of angular momentum transport is parameterized by the dimensionless a-viscosity (Shakura & Sunyaev 1973) which for a circumstellar disk may have a broad distribution ranging from ~10−5–10−1 (Rafikov 2017; Ansdell et al. 2018; Villenave et al. 2022).

Disk gas with a sufficiently high ionization fraction couples to the stellar or planetary magnetic field such that the magnetorotational instability (MRI) induced turbulence may provide the source of the effective viscosity and produce α ≥ 10−3 (Balbus & Hawley 1991; Hawley et al. 1995; Balbus 2003). In the case of a CPD, MRI induced-turbulence may be inhibited by the short orbital periods and presence of magnetic dead-zones, limiting gas transport to a thin surface layer (Fujii et al. 2014). Even if the CPD were effectively inviscid, tidal interaction with the star may promote a minimum rate of angular momentum transport through the excitation of spiral waves (Rivier et al. 2012). In our ProDiMo model, the α-viscosity of the disk is not explicitly specified. Instead, a mass accretion rate is specified which controls the heating rate of the disk through viscous dissipation.

The disk mass, accretion rate, and viscosity are highly degenerate properties. 3D hydrodynamical modeling of gas delivery into the vicinity of a Jupiter mass planet at 5 au suggest M = 10−9.3 M yr−1 of gas (Szulágyi et al. 2022). Stellar tidal perturbation may produce a minimum accretion rate 10−9.7 M yr−1 (Rivier et al. 2012). In the PDS 70 system, two massive planets are observed to be accreting gas in a large dust cavity (Wagner et al. 2018; Keppler et al. 2018; Haffert et al. 2019). K-band observations of PDS 70 b with the VLT are consistent with M = 10−10.8–10−10.3 Myr−1 (Christiaens et al. 2019) with similar values estimated for PDS 70 c (Thanathibodee et al. 2019). HST UV and Hα imaging of the protoplanet PDS 70 b suggest M = 10−10.9–10-10.8 M yr−1 (Zhou et al. 2021). Based on these observational and theoretical constraints we adopt M = 10−10 M yr−1 (with a heating rate corresponding to α ≈ 10−2.7) and M = 10−11 M yr−1 (α ≈ 10−3.6) for the high-mass CPD, representing viscous timescales of 103 and 104 yr, respectively, over which the majority of the CPD mass is replaced by freshly accreted material (see Eq. (1)). For the low-mass CPD we adopt M = 10−11 M yr−1 and M = 10−12 M yr−1 to represent the same α-viscosities and tvisc.

The viscous heating is determined according to the method of D’Alessio et al. (1998). The half-column heating rate is Fvis(r)=3GMpM˙8πr3(1Rp/r)erg cm2s1,${F_{{\rm{vis}}}}\left( r \right) = {{3G{M_{\rm{p}}}\dot M} \over {8\pi {r^3}}}\left( {1 - \sqrt {{R_{\rm{p}}}/r} } \right){\rm{erg}}\,{\rm{c}}{{\rm{m}}^{ - 2}}{{\rm{s}}^{ - 1}},$(10)

where G is the gravitational constant, Mp is the planet mass, r is the distance to the planet, and Rp is the planetary radius. We convert the surface-heating to a heating rate per volume as Γvis(r,z)=Fvis(r)ρd(r,z)ρd(r,z)dzerg cm3s1,${{\rm{\Gamma }}_{{\rm{vis}}}}\left( {r,z} \right) = {F_{{\rm{vis}}}}\left( r \right){{{\rho _{\rm{d}}}\left( {r,z} \right)} \over {\int {{\rho _{\rm{d}}}\left( {r,z'} \right){\rm{d}}z'} }}{\rm{erg}}\,{\rm{c}}{{\rm{m}}^{ - 3}}{{\rm{s}}^{ - 1}},$(11)

where ρd is the mass density of the dust at radius r and height z. The heating is applied directly to the dust. The resulting mid-plane dust temperature profile of the CPDs can be found in Fig. 4. The temperature profiles have been calculated using a new diffusion-approximation radiative transfer solver, which is described in Appendix A.

Sources of CPD external irradiation. The rate of photodesorption plays an important role in determining the ice abundance in the outer optically-thin region of the CPDs. Potential sources of radiation include the planet, the star, and nearby massive cluster stars. For the planet we assume the runaway gas accretion phase has ended and that the luminosity has correspondingly declined to 10−5 L with a surface temperature of 1000 K where it remains relatively constant for 10 Myr (Hubickyj et al. 2005; Marley et al. 2007; Spiegel & Burrows 2012).

We parameterize the isotropic background radiation intensity with the dimensionless parameter χ. The background intensity is then the sum of a diluted 20 000 K blackbody field and the cosmic microwave background, Ivbg=χ·1.71·WdilBv(20000 K)+Bv(2.7 K),$I_v^{{\rm{bg}}} = \chi \,\cdot\,1.71\,\cdot\,{W_{{\rm{dil}}}}{B_v}\left( {20\,000\,{\rm{K}}} \right) + {B_v}\left( {2.7\,{\rm{K}}} \right),$(12)

where the dilution factor Wdil = 9.85357 × 10−17 such that a value of χ = 1 corresponds to the quiescent interstellar background, or “unit Draine field” (Draine & Bertoldi 1996; Röllig et al. 2007; Woitke et al. 2009). Irradiation by the young sun provides a minimum χ ~ 3000 at 5au (see Fig. 1) despite partial shadowing by an inner disk (Oberg et al. 2020). We adopted the same value for the strength of the FUV irradiation in the midplane at Jupiter’s location although it is contingent on our assumptions regarding the stellar UV luminosity and geometry of the inner circumstellar disk. Independently of these factors, Oberg et al. (2020) find that interstellar radiation results in a mean χ of O(3) in the gap, as the Sun is believed to have formed in a relatively dense stellar cluster (Adams 2010; Portegies Zwart 2019) containing massive OB stars (Pfalzner 2013). External irradiation heats dust and gas on the surface and in the outer regions of the optically-thin CPD midplane. 3D dust radiative transfer modelling of gap-embedded CPDs suggests scattered stellar radiation can be the dominant heating source in the outer regions of a CPD (Portilla-Revelo et al. 2022). We assumed that the outer edge of the CPD is in thermal equilibrium with the surroundings and set a CPD background temperature lower limit of 50 K equal to that of dust in the circumstellar disk gap (see Fig. 1).

CPD dust mass and grain size population. The dust-togas ratio is a key parameter that regulates the eventual ratio of ice to rock in the CPD. Planet-induced pressure bumps at the gap edges of the circumstellar disk may filter out dust particles above ~100 µm in size (Rice et al. 2006; Zhu et al. 2012). For a dust grain size population power law exponent of −3.5, minimum grain size amin 0.05 µm, maximum grain size amax 3000 µm, a filtering of grains larger than 100 µm would be reduced the mass of dust entering the gap by a factor ~30. As an opposing view Szulágyi et al. (2022) find that a planet within the gap can stir dust at the circumstellar disk midplane and produce a very high rate of dust delivery to the CPD, such that the dust-to-gas ratio can even be enhanced in the CPD to ~0.1 for a Jupiter mass planet at 5 au.

The dust population within the CPD may also rapidly evolve. It is possible that mm-sized grains are quickly lost to radial drift in 102–103 yr (Zhu et al. 2018). Similarly Rab et al. (2019) used the dust evolution code two-pop-py (Birnstiel et al. 2012) to show that for an isolated CPD onto which new material does not accrete, an initial dust-to-gas ratio of 10−2 can in only 104 yr be reduced to <10−3 and in 105 yr to <10−4. However, larger dust grains may become trapped in CPDs (Drążkowska & Szulágyi 2018; Batygin & Morbidelli 2020). Additionally, we expected that in an embedded, actively-accreting CPD the dust will continually be replenished and that a higher steady-state dust-to-gas ratio will be achieved. Given these considerations we tested various possible dust to gas ratios ranging from 10−4–10−2 in Appendix H.

Table 1

Parameters for the solar circumstellar disk.

thumbnail Fig. 1

2D dust temperature distribution around the circumstellar disk gap in the range 40–120 K (top). Solid, dashed, and dotted black contours indicate the relative UV field strength χ (see Eq. (12)). The white circle with black border represents the position of Jupiter and the horizontal bar indicates the physical extent of the CPD. The white cross with black border represents the location from which the inherited chemistry is extracted at one scale height z = H. Hydrogen nuclei column density (blue line) with unperturbed circumstellar disk density profile for comparison (light blue dashed line; bottom).

Table 2

Conditions at the gap edge “inheritance” point of the circumstellar disk at heliocentric distance 8.2 au, altitude 0.5 au (1 pressure scale height) above the midplane.

thumbnail Fig. 2

Ice abundance in the circumstellar disk gap normalized to the initial abundance of ice for various stellar luminosities as function of time after the onset of UV exposure. The width of each track represents the range of ice sublimation rates corresponding to variable conditions across the gap region (z = 0–0.5 au, r = 5.2–8.2 au).

Table 3

Variation of parameters for the circumplanetary disk model grid.

Table 4

Parameters common to the CPD models.

thumbnail Fig. 3

Surface density (red) and midplane FUV field strength in units of the Draine field χ (blue) of the high (10−7 M) and low (10−8 M) mass CPDs. The four circles indicate the present-day radial position of the Galilean satellites and their position on the ordinate is arbitrary. The striped gray region on the left indicates an empty inner magnetic cavity and the light gray region on the right indicates the gravitationally unstable zone outside of RH/3 (where RH = 0.35 au).

thumbnail Fig. 4

Radial midplane dust temperature profiles of the four reference CPDs (d/g = 10−33), labelled by the model id’s in Table 3. The four circles indicate only the present-day radial position of the Galilean satellites. The striped gray area on the left indicates the inner cavity. The light gray area on the right indicates the gravitationally unstable zone outside of RH/3 (where RH = 0.35 au).

2.2 Dust grain drift within the CPD

In ProDiMo, the chemistry is solved for each grid cell without accounting for radial dust or gas transport. Instead we used properties of the ProDiMo model output to inform the radial dust drift calculations in a post-processing step to compare timescales of chemical evolution vs. dust drift. Disk gas which is pressure-supported orbits at sub-Keplerian velocities such that larger grains feel a headwind and rapidly drift inwards (Weidenschilling 1977). In a circumstellar disk the degree of sub-Keplerianity can be <1% while in a CPD it could be as high as 20–80% due to significant gas pressure support (Armitage 2007; Drążkowska & Szulágyi 2018). We considered whether the timescale of grain drift can compete with the processes that shape grain composition In the Epstein regime, the radial grain drift velocity vr,d can be approximated as vr,d=vr,gTs1ηvKTs+Ts1,${v_{{\rm{r,d}}}} = {{{v_{{\rm{r,g}}}}T_{\rm{s}}^{ - 1} - \eta {v_{\rm{K}}}} \over {{T_{\rm{s}}} + T_{\rm{s}}^{ - 1}}},$(13)

where vr,g is the radial velocity of the gas, Ts is the dimensionless stopping time of a grain, Ts=ts(vkr),${T_{\rm{s}}} = {t_{\rm{s}}}\left( {{{{v_{\rm{k}}}} \over r}} \right),$(14)

where ts is the stopping time, vk is the Keplerian orbital velocity, and r is the radius in the CPD. The stopping time is ts=(ρgrainρgas)(avth),${t_{\rm{s}}} = \left( {{{{\rho _{{\rm{grain}}}}} \over {{\rho _{{\rm{gas}}}}}}} \right)\left( {{a \over {{v_{{\rm{th}}}}}}} \right),$(15)

where a is the grain size, ρgas is the gas density, ρgrain is the material density of the dust grains, and vth is the thermal velocity of the gas: cs (8/π)0.5 where cs is the speed of sound. The parameter η is η=n(cs2vk2),$\eta = n\left( {{{c_{\rm{s}}^2} \over {v_{\rm{k}}^2}}} \right),$(16)

and n is the power law exponent of the radial pressure gradient (Armitage 2010). We extract the gas density ρgas, sound speed cs, and pressure gradient from our disk models to consistently determine the grain drift velocities for a grain material density ρgrain = 3 g cm−3. The Epstein regime is valid where the dust particle size a ≤ 9λ/4 where λ is the mean free path of the gas particles. In the inner CPD the gas density is sufficiently high that this condition can be violated. For the high-mass CPD this occurs inside of the orbit of Callisto and for the low-mass CPD inside of the orbit of Europa for grains less than 1 cm in size, in which case we transition to the Stokes regime and calculate the drift velocities accordingly (Weidenschilling 1977).

We adopted several simplifying assumptions regarding the radial velocity structure of the gas. The centrifugal radius rc of the CPD is the orbital radius at which the specific angular momentum is equal to the average of the infalling material and where solid material is accreted onto the CPD (Canup & Ward 2002; Sasaki et al. 2010). In the case of a Jupiter-mass planet this lies near the orbit of Callisto. Rather than accreting at precisely this radius infalling matter have some intrinsic spread in angular momentum (Machida et al. 2008). We adopted the prescription of Mitchell & Stewart (2011) that material will accrete between 16 and 28 RJ. The gas falling onto the CPD at rc will viscously spread radially away from where it is accreted (Pringle 1981). Hence interior to rc gas flows towards the planet and exterior to rc it will flow outwards. Recently it has been proposed that to effectively trap dust and allow for planetesimal growth within the CPD, gas may indeed need to flow predominantly away from the planet and thus form a decretion disk (Batygin & Morbidelli 2020).

In our high viscosity models, a parcel of gas that accretes onto the CPD near rc and flows outwards must travel ~0.3 au within 103 yr to be consistent with tvisc and thus has a mean outwards radial velocity on the order of 1 ms−1. For our low viscosity case vr,gas ~ 0.1 m s−1.

3 Results

In our reference model grid (Table 3), we considered the case of a higher mass, optically thick CPD (10−7 M), and lower mass, optically thin CPD (10−8 M) both with d/g = 10−3.3. For each unique mass we considered viscous timescales of 103 and 104 yr. For each of the four resulting CPDs we tested two initial compositions: of chemical inheritance and reset, for a total of eight models.

thumbnail Fig. 5

Rate of water ice formation as a function of time at different radii beyond the snowline in the high-mass, high-viscosity CPD (7–10), to illustrate the two distinct stages of water ice formation. All values are normalized to the maximum formation rate at 10−5 yr, r = 0.08 RH. The times where the decline in the abundance of free O and free H end are indicated by the gray vertical lines. The starting time t0 is defined as when the infalling gas is shock-heated to an atomic and ionized state.

3.1 Timescales of ice formation

In the following section we discuss by which pathways and at what rate water (ice) formation occurs in a chemically reset CPD. Thereafter we contrast these results with that of the CPDs which inherit an initial chemical composition from the circumstellar disk.

3.1.1 Reset

In our “reset” scenario the species initially present are exclusively atomic (with the exception of PAHs) and singly or doubly ionized. All hydrogen is initially present in the form of H+. The reset case is characterized by the formation of ice where it is stable against desorption. Ice formation is suppressed in the innermost region of the CPD due to dust temperatures in excess of 160–180 K. In the outer region of the CPD ice formation is suppressed by the low optical depths and correspondingly high intensity of background radiation which causes ices to photodesorb. Between these two boundaries ices begin to form.

The freeze-out occurs step-wise, with two characteristic timescales on which ice formation proceeds. The rate of water ice formation at different disk radii as a function of time is shown in Fig. 5. Within 0.01–1 yr the first step is completed. This initial, rapid water formation occurs primarily via a path that begins with a three-body recombination reaction important at temperatures below 2800 K (Hidaka et al. 1982; Tsang & Hampson 1986), H+O+MOH+M,${\rm{H}} + {\rm{O}} + {\rm{M}} \to {\rm{OH}} + {\rm{M,}}$(17)

where the third body M = H, H2, or He, and to a lesser extent by H + O → OH + photon. The water then forms by radiative association of the OH with free hydrogen; OH+HH2O+photon,${\rm{OH}} + {\rm{H}} \to {{\rm{H}}_2}{\rm{O}} + {\rm{photon,}}$(18)

after which it adsorbs to a grain. Water formation via the reactions (17) and (18) remains proportional to the declining abundance of free H and ends when it is depleted around 10−1 yr. Typically half of the maximum possible amount of water ice that could form is produced during this stage.

The first stage of water ice formation proceeds inside-out as a result of the strong density-dependence of the collider reaction. The second stage proceeds outside-in due to the pathway’s dependence on intermediate species produced by photo-reactions. This can be seen in Fig. 5 where the inner disk formation rate (red lines) is initially higher while later the outer disk rate (blue lines) is relatively higher. At this lower rate of formation the total mass of water ice doubles at the midplane within ~ 1 Myr, approximately half of which forms by NH2+NON2+H2O,${\rm{N}}{{\rm{H}}_2} + {\rm{NO}} \to {{\rm{N}}_2}{\rm{ + }}{{\rm{H}}_2}{\rm{O,}}$(19)

near the snowline. The second stage of ice formation also involves the freezing-out of NH3 and other more volatiles species. We find that the relatively unstable NH2 exists in abundance at such high densities (nH ~ 1015 cm−3) due to the three-body reaction N + H2 + M → NH2 + M. Although the reaction rate has been experimentally determined k0 = 10−26 cm6 s−1 (Avramenko & Krasnen’kov 1966), it has been noted as being in need of revisiting due to its importance in water formation via NH2 + O → NH + OH (Kamp et al. 2017). We have adopted a rate more typical of collider reactions (10−30 cm6 s−1), which still produces enough NH2 for this path to play an important role. The other half of water ice is formed via the more “classical” route H3O++eH2O+H,${{\rm{H}}_3}{{\rm{O}}^ + } + {{\rm{e}}^ - } \to {{\rm{H}}_2}{\rm{O}} + {\rm{H,}}$(20)

and in the outermost part of the disk where water ice is stable via H3O++H2COH3CO++H2O.${{\rm{H}}_3}{{\rm{O}}^ + } + {{\rm{H}}_2}{\rm{CO}} \to {{\rm{H}}_3}{\rm{C}}{{\rm{O}}^ + } + {{\rm{H}}_2}{\rm{O}}{\rm{.}}$(21)

In the high-mass CPDs water ice formation is typically still ongoing by the end of the viscous timescale, and so the mid-plane ice abundance is not able to converge to the level seen in steady-state within the viscous timescale.

3.1.2 Inheritance

The evolution of the inheritance case is characterized by the desorption of ices in regions where they are not stable against thermal or photo-desorption. Where water ice is stable very little additional water formation occurs within the viscous timescale.

Icy grains interior to the snowline sublimate typically within 10−5 yr, and a “snowline” (water iceline) is clearly established. In some cases the snowline can take longer to stabilize, shifting outwards over time, and sometimes continues to evolve radially for up to 105 yr. This is notable in the (8–11) CPD where the snowline moves from 0.01 RH at 10−5 yr to 0.03 RH after 104 yr. In practice, this implies that there exists a radial span in which the composition of radially drifting icy grains does not immediately begin to reflect the ambient conditions. This is discussed in Sect. 4.2.

In the outer optically thin region of the CPDs, ices are also lost to photodesorption although the process is slower than the thermal desorption occurring in the inner disk. Desorption in this area is typically complete within 103–105 yr and has not always reached steady-state by the end of the viscous timescale.

3.2 The midplane ice mass fraction

While water ice can form efficiently in chemically reset CPDs, the final fice of the solids depends strongly on the total dust mass in the midplane. We explore the role of the global d/g ratio in the reset and inheritance cases in Fig. H.1. A canonical dust-to-gas ratio of 10−2 produces at most grains with an ice mass fraction of <0.1 and is nowhere consistent with the composition of the icy Galilean satellites. In contrast a dust-to-gas ratio of 10−3.3 results in solids with fice both above and below the maximum Galilean fice for the high and low-mass CPDs, respectively. Hereafter the global dust-to-gas ratio of 10−3.3 is considered in discussions of the four reference CPDs.

We show the state of the radial ice mass fraction for the CPDs with d/g = 10−3.3 on their respective viscous timescales to allow for direct comparison between the inheritance and reset cases in Fig. 6. For completeness we include the results where grain-surface chemical reactions are utilized.

The midplane fice is also dependent on the degree of dust sedimentation (settling). A general feature of the fice profiles is an inner maximal peak at the snowline followed by a decline towards the outer edge of the disk. As dust settling is more efficient at larger radii, fice reduces accordingly with increasing radius. Settling of dust to the midplane is counteracted by stochastic advection by turbulent eddies in the gas. We assume that the value of turbulent-α used to determine the degree of settling is equal to the global viscous α used to determine the heating rate by viscous dissipation. In the low viscosity CPDs (7–11) and (8–12) dust settling towards the midplane is thus proportionally enhanced. In the low-viscosity cases the dust density is enhanced in the midplane minimally by a factor ~3 over the global d/g, increasing to a factor ~20 at RH/3. As the degree of settling is also dependent on the adopted surface density power law exponent e we explore the impact of deviation from the assumed ϵ = 1 in Appendix C. Given that we have no reason to believe this value will depart significantly from the range 1.0–1.3, the resulting fice in the inner disk should differ from our reference result by no more than 25–30% interior to the ammonia iceline at ~70 K.

In general the high-mass chemically reset CPDs (7–10) and (7–11) are not able to converge entirely towards a steady-state ice abundance in either the “fast” (103 yr) or “slow” (104 yr) viscous timescales as gas-phase CO is more stable and contains a relatively larger fraction of the total oxygen budget for longer. As a result chemically reset CPDs contain on average less water ice than those which inherit their ices from the circumstellar disk. In contrast, the low-mass chemically reset CPDs (8–11) and (8–12) are able to converge towards the ice abundances seen in the inheritance cases within 100 yr.

The role of surface chemistry. The duration of the initial stage in which water formation is rapid is dependent on the availability of atomic H. When H2 formation is complete this stage ends. The formation of H2 is treated differently with the inclusion of grain surface chemistry. When the chemistry is limited to gas-phase reactions and direct adsorption/desorption only, H2 formation proceeds via the analytic rate determined by Cazaux & Tielens (2002). When surface chemistry is included H2 formation is instead modeled explicitly via reactions involving hydrogenated PAHs (H-PAH; Boschman et al. 2012; Thi et al. 2020a).

A chemical reset poses a scenario in which H2 and H2O formation occur simultaneously. The analytic rate of Cazaux & Tielens (2002) presupposes that chemisorbed H plays a role in H2 formation on silicate or carbonaceous surfaces, in which H goes through an intermediate stage of being chemically, rather than physically, bound to the grain surface. We find that prior to atomic H depletion several (>3) monolayers of H2O have formed on the average sized grain. The formation of H2 via chemisorbed H should thus be suppressed in these regions as the water layers prevent H atoms from chemisorbing to the grain surface (Wakelam et al. 2017). In the absence of chemisorbed H, H2 formation on dust is strongly reduced and H2 formation proceeds primarily via H-PAH. The H2 formation rate under these conditions is less efficient than the analytic rate from Cazaux & Tielens (2002) near the inner edge of the snowline at 150–170 K. The relatively slower formation of H2 and the resulting prolonged availability of atomic H results in a ~30–100% increase in water ice abundance interior to the NH3 iceline prior to tvisc and hence narrows the gap between the inheritance and reset cases in this region. Water ice formation in the inner disk is further enhanced by the inclusion of O2H in the surface chemistry network. Gas-phase three-body reactions with O2 produce O2H which in turn lead to early OH formation. The gas-phase O2 reservoir is thus depleted and efficiently converted via OH into H2O via the three-body reaction O2+H+MO2H+M,${{\rm{O}}_2} + {\rm{H}} + {\rm{M}} \to {{\rm{O}}_2}{\rm{H}} + {\rm{M,}}$(22)

with rates adopted from UMIST2006 (Atkinson et al. 2004; Woodall et al. 2007), which is highly efficient at the densities in the CPD, and thereafter O2H+HOH+OH,${{\rm{O}}_2}{\rm{H}} + {\rm{H}} \to {\rm{OH}} + {\rm{OH,}}$(23) OH+HH2O+photon.${\rm{OH}} + {\rm{H}} \to {{\rm{H}}_2}{\rm{O}} + {\rm{photon}}{\rm{.}}$(24)

In the outer disk, the ice mass fraction can be enhanced relative to the gas-phase chemical network as the freezing-out of more volatile species is facilitated by grain surface reactions. CO2 ice is readily formed on the surface via OH + CO → CO2 + H (Oba et al. 2010; Liu & Sander 2015) for which we adopt an effective barrier of 150 K (Fulle et al. 1996; Ruaud et al. 2016). The formation of carbon bearing ices begins to significantly influence the fice of the chemically reset CPDs only after 103 yr and hence the effect on the high-viscosity CPDs with tvisc = 103 yr is less pronounced. The formation and composition of these ice impurities will be discussed in detail in an accompanying work (Oberg et al., in prep.).

thumbnail Fig. 6

Midplane ice fraction at t = tvisc for the CPDs with d/g = 10−33. The four circles indicate the present day radial location of the Galilean satellites and the uncertainty bars represent their possible range of ice fractions. The thin dotted line indicates the initial radial ice fraction in the inheritance cases. The filled gray region on the right indicates the gravitationally unstable zone outside of RH/3.

thumbnail Fig. 7

Radial drift velocity of particles (upper panel) and radial velocity of the gas (lower panel) in the high-mass, high-viscosity CPD (7–10). The position of the Galilean satellites on the ordinate is arbitrary. The gray region on the right indicates the tidally unstable region beyond RH/3. The striped region on the left indicates the inner cavity. Solid lines indicate the case of a decretion disk where gas falls onto the CPD at the centrifugal radius (0.03 RH). Dashed lines indicate the grain drift and gas radial velocity in the case of a pure accretion disk.

3.3 Grain drift vs. adsorption and desorption

We calculated the velocity of radially drifting grains in the four reference CPDs and showcase the results for the high-mass high-viscosity CPD (7–10) in Fig. 7. We solved for the total time it takes grains of several sizes to reach the inner disk edge. The resulting timescale of grain drift can be seen in Fig. 8 which shows the time for a grain deposited at radius r to reach the CPD inner edge. Grains which become trapped are not included in Fig. 8, as they do not reach the disk inner edge. The regions where thermal desorption, ice formation, and photodesorption predominantly shape grain mantle composition, as well as their respective timescales, are indicated in the figure. These are the timescales with which grain drift competes.

The regions have been defined as follows: the “thermal desorption region” is found interior to the snowline where inherited, initially icy grains lose their icy mantles within the viscous timescale. The “snow border” is the region where the reset and inheritance cases are not able to converge towards a common snowline within 106 yr. Grains here are able to retain their icy mantles but no significant additional ice adsorption occurs. The “ice formation region” is where water begins to adsorb to grains after the CPD has been chemically reset. The “photodesorption region” is the optically thin region exterior to the snowline where inherited grains eventually lose their icy mantles within ~106 yr at most.

We focus on grains of size <10 mm as the adsorption and desorption timescales derived in Sect. 3.1 have only been derived with the thermochemical disk model for grains up to this size. In most cases the grain drift timescale tdrift is longer than the timescales of thermal desorption and the timescale of rapid ice formation. The composition of grains will thus correspond to the fice profiles derived in Sect. 3 except in the case of the low-mass, high-viscosity CPD (8–11). In this CPD icy grains can cross into the “thermal desorption region” but only begin desorbing once they approach the position of Europa.

Dust traps. It is clear from Fig. 7 that in a decreting CPD dust traps are present. Grains deposited near the centrifugal radius drift outwards with the gas until the force of radial advection is balanced by the loss of orbital energy from drag against gas orbiting at sub-Keplerian velocity. Trapped grains thus become radially segregated by size, with smaller grains drifting to the outer edge of the trap and larger grains remaining near the inner edge. In the high-mass high-viscosity CPD (7–10) grains 0.1–1 mm in size become trapped near 0.1–0.2 RH. Grains smaller than 0.05 mm advect with the gas and are able to reach the outermost stable orbit at RH/3 where they would be eventually lost to e.g. tidal stripping.

In general, the dust traps are spatially coincident with the ice formation region and in the lower-mass CPDs also partially with the photodesorption region. Hence continued ice deposition on trapped grains could facilitate grain growth. This phenomena is discussed further in Appendix F.

4 Discussion

We set out to explore the process of ice formation in CPDs with relatively short viscous timescales to constrain their physical, chemical, and dynamical properties. We find that even if infalling gas and ice is fully atomized, re-freezing proceeds quickly in the CPD and solids reach an fice ~ 0.5 by tvisc for an appropriate midplane dust-to-gas ratio. The midplane ice abundance at tvisc is generally insensitive to the initial chemical conditions. Only in the inner disk (r < 0.05–0.1 RH) of the high-mass CPDs is ice formation too slow for the reset and inheritance cases to converge. With our standard chemical network the efficiency of water production in this region is closely tied to the availability of atomic H and thus to the H2 formation rate, which is not well constrained in these conditions. The three-body reaction H + O + M → OH + M is also critical to the process. For this reaction we have adopted the rate coefficients listed in UMIST2006 (Woodall et al. 2007), the value of which originates from a recommendation of Tsang & Hampson (1986) who noted that literature values represented only rough estimates and suggested a factor 10 uncertainty (Baulch 1972; Day et al. 1973). More recent estimates of this rate at high temperatures (>3000 K) suggest this value is accurate to within ~40% (Javoy et al. 2003), but modern low temperature measurements remain desirable. In our expanded grain-surface chemical network, gas-phase O2H plays an important role in accelerating OH formation in the inner disk, diverting atomic oxygen into water rather than O2.

4.1 Constraints on CPD properties

For reasonable assumptions regarding the properties of the CPD the snowline can fall very near the present-day position of Europa. The CPD with mass 10−7 M (10−4 MJ), α = 10−3.6 and d/g = 10−3.3 matches best the compositional gradient of the Galilean moons at their present-day orbits. While this seems like a promising outcome, we emphasize that both inwards gas-driven migration (Canup & Ward 2002; Madeira et al. 2021) and long-term outwards tidally-driven migration (Yoder 1979; Tittemore 1990; Showman & Malhotra 1997; Lari et al. 2020) have potentially played a role in repositioning the satellites both during and post-formation. Other regions of the CPD parameter space are equally capable of producing solids with 0.4 < fice < 0.55, but vary in the position of their snowline. In any case, we believe that it is more meaningful to determine whether a particular CPD can form enough ice on the relevant viscous timescale, rather than at precisely which radii a particular abundance of ice can be found.

To produce solids with a minimum fice ~ 0.5 the global dust-to-gas ratio of the CPD must be ≤10−3 and does not need to be <10−3.7. We suggest that this does not represent a minimum d/g limit as imperfect accretion (Dwyer et al. 2013), (hydrodynamic) proto-atmospheric escape (Kuramoto & Matsui 1994; Bierson & Nimmo 2020), impact-vaporization (Nimmo & Korycansky 2012), or tidal heating (Canup & Ward 2009) may all have played a role in reducing the volatile mass of the satellites either during or post-accretion. A minimum d/g is instead implied by the minimum time required to accrete the solid total mass of the Galilean satellites (~10−7 M) into the CPD. We consider the lifetime of the gaseous circumstellar disk (~10 Myr) to be the maximum time over which the CPD can continue to accrete gas. Assuming Jupiter’s runaway accretion ended <3.94 Myr after CAI formation (Weiss & Bottke 2021) and that moon formation only began at this stage, this leaves ~6 Myr to accrete the refractory material for the moons. We emphasize that this limit is very approximate, as it is possible that circumstellar gas disk lifetimes regularly exceed 10 Myr (Pfalzner et al. 2014). Conversely an even shorter lifetime (<4.89 Myr after CAI formation) has been proposed for the solar nebula on the basis of paleomagnetic measurements (Weiss & Bottke 2021).

The midplane dust-to-gas ratio and thus the fice in the CPD will differ from what has been derived from our models if the grain size distribution is significantly altered in some way. The circumstellar disk gap edge may filter out larger grains from the accretion flow (Rice et al. 2006; Zhu et al. 2012). Grains larger than 100 µm, which settle efficiently, are those primarily responsible for enhancement of the dust mass at the CPD midplane. Massive planets may however vertically stir circumstellar disk midplane dust outside the gap (Binkert et al. 2021). Szulágyi et al. (2022) found that significant vertical stirring of large grains occurs at the gap outer wall in the presence of planets with mass above that of Saturn, resulting in a substantial delivery of mm-sized grains to the CPD. The high-mass tail of the dust distribution could alternatively be depleted by the more rapid inwards drift of these grains in the CPD. We have shown that grains larger than ~1 mm no longer advect with the gas in the CPD. The steady-state dust grain size distribution will thus likely be truncated. In the absence of these grains the limits we have derived on the CPD mass are revised upwards by a factor ~5.

thumbnail Fig. 8

Inwards drift timescale for dust grains (colored lines) relative to the timescales of desorption and ice formation (black dashed lines) in the reference CPDs with d/g = 10−33. From left to right the colored regions correspond to areas where all ices are eventually thermally desorbed (pink), the iceline (light yellow), where grains become icier (light blue) and where all ices are eventually photodesorbed (dark gray). The red line corresponds to the minimum grain size which is not trapped in the CPD (if there is a sign change in the gas radial velocity) and the blue line to the maximum for which the timescales of desorption and ice formation have been verified. The vertical dotted line indicates the approximate outer region of gravitational stability in the CPD.

4.2 Does grain drift erase the radial distribution of ices?

We have tested only the “gas-starved” disk paradigm in which a CPD must over time accrete the solids to form large moons. The sequential formation, episodic growth, and potential loss of migrating moons is a characteristic of this theory (Canup & Ward 2002). In such a dynamical environment the relevancy of the instantaneous radial distribution of icy grains remains to be established. The simplest way in which the chemical properties of dust in the CPD could be imprinted on the final satellite system would be through in-situ formation: the satellites accrete the bulk of their material at fixed radial positions in the CPD (Ronnet et al. 2017). This might occur if the innermost proto-moon were prevented from migrating into Jupiter by the presence of a gas-free magnetically-truncated inner cavity (Takata & Stevenson 1996; Batygin 2018; Shibaike et al. 2019). Additional proto-satellites could then pile-up in a resonant chain and be stabilized against further migration by the proto-moon anchored at the disk inner edge (Ogihara & Ida 2009; Sasaki et al. 2010; Izidoro et al. 2017; Madeira et al. 2021). Drifting grains would still be free to cross the orbit of proto-moons as accretion efficiency remains low when the proto-moons are only a fraction of their final mass (Shibaike et al. 2019). In this paradigm proto-moons at relatively fixed positions, continually accreting grains that drift into their feeding-zone (Ronnet & Johansen 2020). We find that the ice fraction of small (<1 cm) drifting grains in the inner disk will almost always reflect ambient conditions (Sect. 3.3) independently of whether the gas in the CPD flows radially inwards or outwards (the fate of trapped grains is discussed in Appendix F). If the proto-moons (resonantly anchored at fixed radii within the CPD) accrete the majority of their total mass from these drifting grains, their bulk ice fraction would reflect the temperature gradient of the CPD.

5 Conclusions

Circumplanetary disks represent a unique chemical environment characterized by high-densities and a relatively short timescale on which gas and dust are viscously or aerodynamically lost. We aimed to explore the process of ice formation in this environment from sharply contrasting initial chemical conditions, knowing that solids with ice/rock ~1 must be able to form within the viscous timescale of a Jovian CPD. We tested the paradigm in which solids are delivered directly from the circumstellar disk in the form of small grains (<1 cm) to a “gas-starved”, relatively low-mass CPD. We highlight our key conclusions as follows: If infalling material is chemically reset:

  1. High densities in the CPD facilitate three-body “collider” reactions that lead to rapid water ice formation. Roughly half of the water ice is formed within a single year by the hydrogenation of OH.

  2. Solids with the ice fraction of Ganymede or Callisto are produced within tvisc for α ≈ 10−3–10−4 if the midplane is depleted in dust by a factor 20–50 relative to the canonical d/g = 10−2.

If chemical inheritance occurs:

  1. Ices near the planet efficiently sublimate and establish a snowline at a similar location to that of the reset case within tvisc. Additional ice formation is minimal.

In either case:

  1. Icy circumstellar dust grains preserve the majority of their volatile content during gap-crossing if accreted onto the CPD within 100 yr unless the stellar luminosity is >10 L.

  2. The compositional imprint of the CPD temperature profile is not erased by radial dust drift for grains of size a < 1 cm.

  3. Only the “high-mass” CPDs (MCPD = 10−4 MJ) are sensitive to the initial chemical conditions: water ice formation in the inner disk is less efficient if a chemical reset occurs as oxygen tends to remain locked in the gas-phase CO.

In our solar system icy moons are common. No matter whether or not ices sublimate upon incorporation into the CPD, we have demonstrated that ices can be efficiently re-deposited onto dust grains and enable the general ubiquity of icy moons.

Acknowledgements

The research of N.O. and I.K. is supported by grants from the Netherlands Organization for Scientific Research (NWO, grant number 614.001.552) and the Netherlands Research School for Astronomy (NOVA). “This research has made use of NASA’s Astrophysics Data System Bibliographic Services. This research has made extensive use of Numpy (Harris et al. 2020), Matplotlib (Hunter 2007), scipy (Virtanen et al. 2020), and prodimopy https://gitlab.astro.rug.nl/prodimo/prodimopy. The authors would like to thank the anonymous referee for comments that contributed to the accuracy, clarity, and focus of this work.


Appendix A The RT diffusion solver

ProDiMo solves the radiative transfer equation (see Eq. (12) in Woitke et al. 2009) together with the condition of radiative energy conservation, which in general can be written as div F=Γ,${\rm{div}}\,{\bf{F}} = {\rm{\Gamma ,}}$(A.1)

where F = ∫ Fv dv [erg/cm2/s] is the bolometric radiation flux vector and Γ [erg/cm3/s] is the non-radiative heating rate per volume. In the viscous case, we use Γ = Γvis, see Eq. (10) with stellar mass M* and stellar radius R* instead of Mp and Rp for circumstellar discs. The additional non-radiative heating leads to a surplus emission of photon energy according to 4πκvabs(Bv(T)Jv)dv=Γ,$4\pi \int {\kappa _v^{{\rm{abs}}}\left( {{B_v}\left( T \right) - {J_v}} \right)dv = {\rm{\Gamma ,}}} $(A.2)

where is the dust absorption coefficient at frequency v, Bv(T) the Planck function, and Jv the mean intensity. For passive discs, we have Γ = 0, in which case Eq. (A.2) simplifies to the ordinary condition of radiative equilibrium.

The numerical solution method for the radiative transfer (RT) in ProDiMo involves iterations where formal solutions with isotropic scattering are computed along multiple rays to cover the full 4π solid angle as seen from every point in the disc, see Sect. 4 in Woitke et al. (2009). A formal solution results in new Jv(r, z) which are used to update the dust temperatures Tdust(r, z) and source functions. The convergence of this Ʌ-iteration is accelerated by the Ng-method (see Auer 1984). However, despite this acceleration, the convergence is still slow in the midplane, which is a serious problem for all radiative transfer codes for discs, including the Monte-Carlo codes, see Pinte et al. (2009).

Here we describe a method how to avoid this problem. In the diffusion approximation, the bolometric radiation flux F(r,z)=4π3κR(r,z)gradJ(r,z)${\bf{F}}\left( {r,z} \right) = - {{4\pi } \over {3{\kappa _{\rm{R}}}\left( {r,z} \right)}}{\bf{grad}}\,J\left( {r,z} \right)$(A.3)

is given by the gradient of the bolometric mean intensity J(r, z). The Rosseland-mean and Planck-mean opacities are defined as 1κR(r,z)=1κvext(r,z)dBv(T)dTdv/dBv(T)dTdv${1 \over {{\kappa _{\rm{R}}}\left( {r,z} \right)}} = \int {{1 \over {\kappa _v^{{\rm{ext}}}\left( {r,z} \right)}}{{d{B_v}\left( T \right)} \over {dT}}dv/\int {{{d{B_v}\left( T \right)} \over {dT}}dv} } $(A.4) κP(r,z)=κvabs(r,z)Bv(T)dv/Bv(T)dv,${\kappa _{\rm{P}}}\left( {r,z} \right) = \int {\kappa _v^{{\rm{abs}}}\left( {r,z} \right){B_v}\left( T \right)dv/\int {{B_v}\left( T \right)dv,} } $(A.5)

where κvext(r,z)$\kappa _v^{{\rm{ext}}}\left( {r,z} \right)$ is the extinction coefficient and T = Tdust(r, z) the dust temperature at position (r, z) in the disk.

At the beginning of a new RT iteration, the Rosseland and Planck opacities are calculated based on the frequency-dependent disk opacity structure and the current Tdust(r, z). Next, we compute radial and vertical Rosseland optical depths τRoss = ∫ κR ds from every point. When the radial inward, radial outward and vertical upward Rosseland optical depths from that point are all larger than a threshold value (we use a value of 10 here), the point is flagged as being optically thick, and added to the subset of optically thick points M={ (i,j)|(ri,zi,j)  is  optically thick }.$M = \left\{ {\left( {i,j} \right)|\left( {{r_i},{z_{i,j}}} \right)\,\,{\rm{is}}\,\,{\rm{optically}}\,{\rm{thick}}} \right\}.$(A.6)

where i and j are the 2D-indices of a grid point at radius ri and height zi,j. The following method only updates the mean intensities Ji,j and dust temperatures Ti,j on the optically thick grid points (i, j) ∈ M, whereas all other points are regarded as fixed boundary conditions for this problem. To pick up the bolometric mean intensities on the boundary points, we integrate Eq. (A.2) assuming that Jv is close to a Planckian, hence J(r,z)=B(T)Γ(r,z)4πκP(r,z).$J\left( {r,z} \right) = B\left( T \right) - {{{\rm{\Gamma }}\left( {r,z} \right)} \over {4\pi {\kappa _{\rm{P}}}\left( {r,z} \right)}}.$(A.7)

thumbnail Fig. A.1

Volume and areas for the spatial cell around point (i, j). The mean intensities J and Rosseland opacities κR are available on the grid points (i, j) marked with red dots.

where B(T) = σ Tdust(r, z)4/π is the frequency-integrated Planck function. Integration of Eq. (A.1) over the volume associated with grid point (i, j) as sketched in Fig. A.1 results in the following numerical equation Ai½,jverDi½,jJi,jJi1,jΔri½+Ai+½,jverDi+½,jJi,jJi+1,jΔri+½+AihorDi,j½Ji,jJi,j1Δzi,j½+AihorDi,j+½Ji,jJi,j+1Δzi,j+½=Vi,jΓi,j,$\matrix{{A_{i - \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/\kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} ,j}^{{\rm{ver}}}{D_{i - \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/\kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} ,j}}{{{J_{i,j}} - {J_{i - 1,j}}} \over {{\rm{\Delta }}{r_{i - \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/\kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} }}}} + A_{i + \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/\kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} ,j}^{{\rm{ver}}}{D_{i + \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/\kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} ,j}}{{{J_{i,j}} - {J_{i + 1,j}}} \over {{\rm{\Delta }}{r_{i + \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/\kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} }}}}} \hfill \cr { + A_i^{{\rm{hor}}}{D_{i,j - \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/\kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} }}{{{J_{i,j}} - {J_{i,j - 1}}} \over {{\rm{\Delta }}{z_{i,j - \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/\kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} }}}} + A_i^{{\rm{hor}}}{D_{i,j + \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/\kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} }}{{{J_{i,j}} - {J_{i,j + 1}}} \over {{\rm{\Delta }}{z_{i,j + \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/\kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} }}}} = {V_{i,j}}{{\rm{\Gamma }}_{i,j}},} \hfill \cr} $(A.8)

where we note that the vertical fluxes through the cell boundaries involve a scalar product with the slanted normal vector of the surface area, hence Aihor$A_i^{{\rm{hor}}}$ is the cell’s horizontal area after projection onto the vertical direction. The following abbreviations are used for the distances, vertical and horizontal areas, and the volume. They are given by the geometry of the ProDiMo grid points, which are aligned on radial rays on which z/r is constant ri½=riri1${r_{i - \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/\kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} }} = \sqrt {{r_i}{r_{i - 1}}} $(A.9) ri+½=riri+1${r_{i + \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/\kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} }} = \sqrt {{r_i}{r_{i + 1}}} $(A.10) Δri12=riri1${\rm{\Delta }}{r_{i - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}} = {r_i} - {r_{i - 1}}$(A.11) Δri+12=ri+1ri${\rm{\Delta }}{r_{i + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}} = {r_{i + 1}} - {r_i}$(A.12) zi,j12=12(zi,j+zi,j1)${z_{i,j - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}} = {1 \over 2}\left( {{z_{i,j}} + {z_{i,j - 1}}} \right)$(A.13) zi,j+12=12(zi,j+zi,j+1)${z_{i,j + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}} = {1 \over 2}\left( {{z_{i,j}} + {z_{i,j + 1}}} \right)$(A.14) zi12,j12=zi,j12ri12ri${z_{i - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}},j - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}} = {z_{i,j - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}}{{{r_{i - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}}} \over {{r_i}}}$(A.15) zi12,j+12=zi,j+12ri12ri${z_{i - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}},j + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}} = {z_{i,j + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}}{{{r_{i - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}}} \over {{r_i}}}$(A.16) zi+12,j12=zi,j12ri+12ri${z_{i + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}},j - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}} = {z_{i,j - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}}{{{r_{i + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}}} \over {{r_i}}}$(A.17) zi+12,j+12=zi,j+12ri+12ri${z_{i + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}},j + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}} = {z_{i,j + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}}{{{r_{i + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}}} \over {{r_i}}}$(A.18) Ai12,jver=2πri12(zi12,j+12zi12,j12)$A_{i - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}},j}^{{\rm{ver}}} = 2\pi {r_{i - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}}\left( {{z_{i - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}},j + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}} - {z_{i - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}},j - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}}} \right)$(A.19) Ai+12,jver=2πri+12(zi+12,j+12zi12,j12)$A_{i + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}},j}^{{\rm{ver}}} = 2\pi {r_{i + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}}\left( {{z_{i + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}},j + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}} - {z_{i - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}},j - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}}} \right)$(A.20) Aihor=π(ri+122ri122)$A_i^{{\rm{hor}}} = \pi \left( {r_{i + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}^2 - r_{i - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}^2} \right)$(A.21) Vi,j=Aihor(zi,j+12zi,j12)${V_{i,j}} = A_i^{{\rm{hor}}}\left( {{z_{i,j + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}} - {z_{i,j - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}}} \right)$(A.22)

The radiative diffusion coefficients are defined as Di,j=4π3κR(ri,zi,j)${D_{i,j}} = {{4\pi } \over {3{\kappa _{\rm{R}}}\left( {{r_i},{z_{i,j}}} \right)}}$(A.23) Di12,j=Di,jDi1,j${D_{i - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}},j}} = \sqrt {{D_{i,j}}{D_{i - 1,j}}} $(A.24) Di+12,j=Di,jDi,j1${D_{i + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}},j}} = \sqrt {{D_{i,j}}{D_{i,j - 1}}} $(A.25) Di,j12=Di,jDi,j1${D_{i,j - {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}} = \sqrt {{D_{i,j}}{D_{i,j - 1}}} $(A.26) Di,j+12=Di,jDi,j+1${D_{i,j + {\raise0.7ex\hbox{$1$} \!\mathord{\left/{\vphantom {1 2}}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$2$}}}} = \sqrt {{D_{i,j}}{D_{i,j + 1}}} $(A.27)

Equation (A.8) states a system of linear equations for the unknown bolometric mean intensities Ji,j on the optically thick points (i, j) ∈ M of the form A·X=B,$A\,\cdot\,{\bf{X}} = {\bf{B}},$(A.28)

where all quantities in Eqs. (A.9) to (A.27) are constants forming the matrix A, and the volumes Vi,j and heating rates Ti,j are constants forming the rest vector B. The unknowns {Ji,j} at the optically thick points (i, j) ∈ M constitute the solution vector X. All terms in Eq. (A.8) that involve the other Ji,j on the adjacent points are also included into B. The matrix equation to solve (Eq. A.28) has a typical dimension of a few hundred to a few thousand, depending on disk mass, geometry, and dust parameters.

This way we can solve the 2D radiative diffusion problem for the unknown mean intensities in the optically thick region as a linear boundary value problem in one go, where there is one layer of points surrounding the optically thick regions which sets the boundary values. Our method calculates how the disk transports the photon energy through the optically thick core inside of the boundary layer. It is applicable to both cases, passive discs without viscous heating and active discs with Γ > 0.

Once the {Ji,j} on (i, j) ∈ M have been determined, we revert the process described by Eq (A.7) B(T)=J(r,z)+Γ(r,z)4πκp(r,z)$B\left( T \right) = J\left( {r,z} \right) + {{{\rm{\Gamma }}\left( {r,z} \right)} \over {4\pi {\kappa _{\rm{p}}}\left( {r,z} \right)}}$(A.29) Tdust(r,z)=(πσB(T))1/4${T_{{\rm{dust}}}}\left( {r,z} \right) = {\left( {{\pi \over \sigma }B\left( T \right)} \right)^{1/4}}$(A.30)

and multiply the frequency-dependent mean intensities Jv(r, z), as they were determined prior to the application of the diffusion solver, by a constant factor to make Eq. (A.2) valid again, thereby keeping the previously calculated frequency distribution of Jv(r,z).

After having modified Tdust(r,z) and Jv(r,z) this way on all grid points (i, j) ∈ M, the normal RT solution method resumes, which begins by calculating the source functions on all grid points and continues by performing a formal solution.

Figure A.2 shows a benchmark test against the Monte Carlo radiative transfer program MCMax (Min et al. 2011). We consider the disk model that is described in detail by Woitke et al. (2016), see their table 3. The central star is a 2 Myr old T Tauri star with a mass of 0.7 M and a luminosity of 1 L, the disk has a mass of 0.01 M, with a gas/dust ratio of 100. The dust is composed of 60% silicate, 15% amorphous carbon and 25% porosity. The dust grains have sizes between 0.05 µm and 3 mm, with an unsettled powerlaw size-distribution of index −3.5. The dust is settled according to the prescription of Dubrulle et al. (1995) with αsettle = 0.01. In contrast to this standard passive T Tauri model, we use here a mass accretion rate of Macc = 10−8 M/yr to set the viscous heating of the disk according to Eq. (10). We use 140 radial × 100 vertical grid points, and 40 wavelength bins. Figure A.2 shows good agreement.

thumbnail Fig. A.2

Benchmark test for a viscous circumstellar disk with a disk mass of 0.01 M and a mass accertion rate of Macc = 10−8M/yr, see text. The green lines are temperature cuts at selected radii calculated by MCMax, the small wiggles are due to the Monte-Carlo noise. The black dots show the temperatures calculated by ProDiMo.

Figure A.2 reveals a number of interesting features in the disk temperature structure:

thumbnail Fig. A.3

Midplane temperature Tdust(r, 0) as function of the log-distance from the inner rim. Again, the green line is the temperature profile calculated by MCMax, and the small black diamonds show the temperature values calculated by ProDiMo.

  • The dust temperature at the inner rim is not much affected by viscous heating.

  • From top to midplane, the temperature first decreases in the disk shadow, but then the trend is reversed and the temperature re-increases towards the midplane as the viscous heat pumped into the disk needs to flow outward, that is mostly upward, which according to the diffusion approximation requires a negative temperature gradient.

  • There is little effect of viscous heating on rdust outside of the optically thick region which extends outward to about 10 au and upward to about z/r ≈ 0.1 – 0.15 in this model.

  • The temperature profile across the midplane beyond the tapering radius (Rtap = 100 au, see Woitke et al. 2016) shows a deep minimum around the midplane z = 0. This is because of the extreme dust settling that occurs in these diluted outer disk regions, creating more optical thickness along the midplane, which brings down the midplane temperature to only 6 K in this model.

Figure A.3 compares the calculated midplane temperatures between ProDiMo and MCMax, which reveals dust temperatures as high as 2800 K, which is of course questionable because at such temperatures, the dust grains are expected to sublimate.

thumbnail Fig. A.4

Maximum relative temperature change between iterations as function of RT iteration number.

Figure A.4 shows the convergence of the RT method, achieving residual relative temperature changes smaller than 10−4 after about 150 RT iterations. The maximums occuring each 5th iteration are due to the Ng-acceleration algorithm.

Appendix B Adsorption energies

Table B.1

Adsorption energies of the most prevalent molecular ices found in our model CPDs.

The adsorption energies of our most common ices and their respective references are listed in Table B.1.

Appendix C Surface density slope

Our reference CPDs have a surface density powerlaw exponent ϵ = 1. The steady-state solution for a constant-M decretion α-disk is ϵ ≈ 1.25 (Batygin & Morbidelli 2020). The midplane ice mass fraction for a variety of possible values of ϵ is shown in Fig. C.1. For a higher ϵ the NH3 iceline responsible for the “bump” in the fice profile at 0.07 RH moves outwards only negligibly. In the inner disk however the ice mass fraction increases due to a combination of the lower midplane dust-to-gas and more efficient H2O formation.

thumbnail Fig. C.1

Midplane ice mass fraction for a variety of surface density powerlaw exponents ϵ for the (7–11) reset CPD with global d/g = 10−33.

Appendix D Background temperature

Throughout this work we have assumed that the CPD is embedded in a radiation field in which the equilibrium dust temperature is 50 K. The midplane dust temperature at the gap center within the circumstellar disk is 50 ± 2 K, for a solar luminosity 0.83 L, gap AV = 0.008, and heliocentric distance 5.2 au. For earlier formation times with correspondingly higher solar luminosities (2.34–13.6 L), we find gap midplane dust temperatures ranging from 70–120 K at 5.2 au.

We have assumed that the final stage of Jupiter’s accretion and moon formation occured at a radial distance from the sun of 5.2 au. Volatile enrichment in Jupiter’s atmosphere indicates it may have formed further out at circumstellar disk temperatures < 25 K or at radii > 30 au (Öberg & Wordsworth 2019). The Nitrogen abundance of Jupiter, approximately 4× solar, may suggest additional N2 was accreted from the solar nebula near the N2 ice-line (Bosman et al. 2019). In light of this possibility we consider also lower background temperatures down to 20 K. The midplane fice for the reference (7–11) CPD can be seen in Fig.D.1. Inside the optically thick region of the CPD the influence of the background temperature Tback is marginal for temperatures ≤ 70 K. Above 70 K the more volatile NH3 and CO2 are unstable as ices and only water ice remains. Below 40 K the outer disk is able to retain ices at radii where AV < 1 as the photodesorption timescale is in excess of the viscous timescale.

thumbnail Fig. D.1

Midplane ice mass fraction for a variety of background temperatures Tback for the (7–11) reset CPD and global d/g = 10−3.3. The four empty circles with error bars represent the Galilean satellites at their present day locations and composition. The gray striped region on the left represents the inner cavity and the light gray shaded region on the right represents the gravitationally unstable zone.

Appendix E Vertical mixing

We have made the simplifying assumption that material which accretes onto the CPD is instantaneously distributed vertically throughout the disk. The shock front may be found at a few (~ 5) scale heights above the CPD midplane (Tanigawa et al. 2012). At 5 scale heights above the centrifugal radius the dust temperature Tdust = 123 K (relative to 89.5 K at the midplane), and optical extinction AV=0.004 (16.4 at the midplane). The velocity of vertical mixing by turbulent diffusion can be estimated as vz = αcs where cs is the local speed of sound (Heinzeller et al. 2011). We find vz ~ 0.5 – 1 m s−1 in this region for the high-viscosity CPDs, assuming that the magnitude of the turbulence is constant from the midplane up to z = 5H. The resulting vertical diffusion mixing timescale is 0.01 tvisc (10-100 yr). We perform a test in the (7–11) CPD in which a parcel of gas is accreted at z = 5H and iteratively evolve its chemistry in steps as it diffuses towards the midplane over 10 yr to understand the impact of more tenuous and high temperature conditions in the initial stages of ice formation (see Fig. E.1). By the end of the stage of rapid formation of water (1–2 yr), the ice abundance has equalized to the fiducial case at the midplane.

thumbnail Fig. E.1

Evolution of the water ice abundance relative to the total number of hydrogen nuclei at the midplane (blue) and at 5 scale heights above the CPD midplane (red), and in the case of the gas parcel which drifts from z = 5H downwards to the midplane over a period of 100 yr (grey). The period in which the altitude of the gas parcel was iteratively lowered towards the midplane is highlighted in light green.

Appendix F Continued ice deposition on trapped grains

In each of the reference CPDs grains of a certain size range remain trapped within the CPD if we assume that gas actively decretes via the outer edge of the disk. A trapped grain will increase physically in size as ice adsorbs onto its surface and thus alter its aerodynamic properties. Batygin & Morbidelli (2020) proposes that grain growth and sublimation could play a role in trapping and radial cycling of grains with size 0.1–10 mm. The size range of trapped grains is ~ 0.01 – 1 mm in our high-mass CPDs and ~ 10−3 – 0.01 mm in the low-mass CPDs, representing 2.5% and 0.1% of the infalling dust mass that reaches the mid-plane, respectively. A modal icy grain is typically coated in no more than 4000 monolayers of water ice. Assuming a monolayer thickness ~ 0.5 nm (Zangi & Mark 2003) and compact morphology, an icy mantle no more than 1 m thick will form. A grain of size 0.05 m can thus increase in size by a factor 20 and have a density corresponding more closely to water ice rather than silicate. For a 1 m thick mantle the aerodynamic effect of increased cross-section is negated by the corresponding reduction in grain density. If the trapped grain icy mantles continue to grow mantles beyond several 1000 monolayers, the new equilibrium trap radius drifts inwards. Realistically only a fraction (0.01-0.1%) of the total CPD gas mass is accreted per year. Mantle growth for a trapped grain will thus not exceed ~ 2 nm yr−1 on average. The time for a grain to grow an icy mantle that allows it to drift to the trap inner edge is then ~ 106 yr (high-mass CPD) or ~ 105 yr (low-mass CPD) assuming a compact grain structure. Ice deposition is thus unlikely to allow grains to escape traps. This estimate does not take into consideration grain growth by coagulation or fragmentation by collision.

Appendix G Chemical abundances of the 0D “molecular cloud” model

The input parameters of the 0D molecular cloud model can be found in Table G.1. A comparison between the model column densities of several common species with observations of TMC-1 can be found in Fig. G.1. While most of the common species column density agrees relatively well with observations, the abundance of S-bearing species hinge on the uncertainties regarding the S elemental abundance (Ruffle et al. 1999).

Table G.1

Molecular cloud parameters

thumbnail Fig. G.1

Ratio of the column density of several common species relative to H2 in the 0D molecular cloud model (circles with error bars) to observational values for TMC-1 (squares with black border). Errorbars on model values represent a factor 10x uncertainty for illustrative purposes. TMC-1 abundances are taken from (Suutarinen etal. 2011) (OH,CH), (Sakai et al. 2010) (C2H), (Smith et al. 2004; Walsh et al. 2009) otherwise.

Appendix H CPD dust-to-gas ratio

In Fig. H.1 we explore the resulting midplane ice mass fraction fice for possible values of the global dust-to-gas ratio from the canonical 10−2 down to 10−4. The maximum grain size and dust population size distribution is kept constant. A global dust-to-gas ratio of 10−33 results in maximum midplane fice values most consistent with the Galilean satellite bulk compositions and hence is adopted as the reference value throughout this work.

thumbnail Fig. H.1

Midplane radial ice mass fraction fice of the dust-to-gas ratio parameter exploration for the chemically reset (left column) and chemically inherited (right column) reference CPDs. The black dotted, dashed, and solid lines indicate the radial position of Europa, Ganymede, and Callisto. The white dotted, dashed, and solid lines indicate where fice is equivalent to the estimated ice mass fraction of Europa (0.05), Ganymede (0.48), and Callisto (0.52).

References

  1. Adams, F. C. 2010, ARA&A, 48, 47 [NASA ADS] [CrossRef] [Google Scholar]
  2. Alibert, Y., Mousis, O., & Benz, W. 2005, A&A, 439, 1205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21 [NASA ADS] [CrossRef] [Google Scholar]
  4. Aota, T., Inoue, T., & Aikawa, Y. 2015, ApJ, 799, 141 [NASA ADS] [CrossRef] [Google Scholar]
  5. Aoyama, Y., Ikoma, M., & Tanigawa, T. 2018, ApJ, 866, 84 [NASA ADS] [CrossRef] [Google Scholar]
  6. Armitage, P. J. 2007, ArXiv e-prints [arXiv:astro-ph/0701485] [Google Scholar]
  7. Armitage, P. J. 2010, Astrophysics of Planet Formation (Cambridge University Press) [Google Scholar]
  8. Atkinson, R., Baulch, D. L., Cox, R. A., et al. 2004, Atmos. Chem. Phys., 4, 1461 [NASA ADS] [CrossRef] [Google Scholar]
  9. Auer, L. 1984, in Numerical Radiative Transfer, ed. W. Kalkhofen (Cambridge: Cambridge University Press), 101 [Google Scholar]
  10. Avramenko, L. I., & Krasnen’kov, V. M. 1966, Bull. Acad. Sci. USSR, Div. Chem. Sci., 15, 394 [CrossRef] [Google Scholar]
  11. Ayliffe, B. A., & Bate, M. R. 2009, MNRAS, 397, 657 [Google Scholar]
  12. Balbus, S. A. 2003, ARA&A, 41, 555 [Google Scholar]
  13. Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [Google Scholar]
  14. Batygin, K. 2018, AJ, 155, 178 [NASA ADS] [CrossRef] [Google Scholar]
  15. Batygin, K., & Morbidelli, A. 2020, ApJ, 894, 143 [NASA ADS] [CrossRef] [Google Scholar]
  16. Baulch, D. L. 1972, Evaluated kinetic data for high temperature reactions (Butterworths) [Google Scholar]
  17. Bierson, C. J., & Nimmo, F. 2020, ApJ, 897, L43 [NASA ADS] [CrossRef] [Google Scholar]
  18. Bierson, C. J., & Steinbrügge, G. 2021, Planet. Sci. J., 2, 89 [NASA ADS] [CrossRef] [Google Scholar]
  19. Binkert, F., Szulágyi, J., & Birnstiel, T. 2021, MNRAS, 506, 5969 [NASA ADS] [CrossRef] [Google Scholar]
  20. Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Boschman, L., Reitsma, G., Cazaux, S., et al. 2012, ApJ, 761, L33 [NASA ADS] [CrossRef] [Google Scholar]
  22. Bosman, A. D., Cridland, A. J., & Miguel, Y. 2019, A&A, 632, A11 [Google Scholar]
  23. Brown, W. A., & Bolina, A. S. 2006, MNRAS, 374, 1006 [Google Scholar]
  24. Canup, R. M., & Ward, W. R. 2002, AJ, 124, 3404 [Google Scholar]
  25. Canup, R. M., & Ward, W. R. 2009, Origin of Europa and the Galilean Satellites, eds. R. T. Pappalardo, W. B. McKinnon, & K. K. Khurana, 59 [Google Scholar]
  26. Cazaux, S., & Tielens, A. G. G. M. 2002, ApJ, 575, L29 [NASA ADS] [CrossRef] [Google Scholar]
  27. Christiaens, V., Cantalloube, F., Casassus, S., et al. 2019, ApJ, 877, L33 [Google Scholar]
  28. Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133 [NASA ADS] [CrossRef] [Google Scholar]
  29. D’Alessio, P., Cantö, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411 [Google Scholar]
  30. D’Angelo, G., Henning, T., & Kley, W. 2002, A&A, 385, 647 [CrossRef] [EDP Sciences] [Google Scholar]
  31. Day, M., Thompson, K., & Dixon-Lewis, G. 1973, Symp. (Int.) Combust., 14, 47 [CrossRef] [Google Scholar]
  32. de Pater, I., & Lissauer, J. J. 2010, Planet. Sci. (Cambridge, UK: Cambridge University Press) [CrossRef] [Google Scholar]
  33. Draine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269 [Google Scholar]
  34. Drążkowska, J., & Szulágyi, J. 2018, ApJ, 866, 142 [CrossRef] [Google Scholar]
  35. Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237 [Google Scholar]
  36. Dwyer, C. A., Nimmo, F., Ogihara, M., & Ida, S. 2013, Icarus, 225, 390 [NASA ADS] [CrossRef] [Google Scholar]
  37. Fujii, Y. I., Okuzumi, S., Tanigawa, T., & Inutsuka, S.-I. 2014, ApJ, 785, 101 [NASA ADS] [CrossRef] [Google Scholar]
  38. Fulle, D., Hamann, H. F., Hippler, H., & Troe, J. 1996, J. Chem. Phys., 105, 983 [NASA ADS] [CrossRef] [Google Scholar]
  39. Garrod, R. T., & Herbst, E. 2006, A&A, 457, 927 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Gressel, O., Nelson, R. P., Turner, N. J., & Ziegler, U. 2013, ApJ, 779, 59 [NASA ADS] [CrossRef] [Google Scholar]
  41. Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749 [Google Scholar]
  42. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
  43. Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 261, 83 [NASA ADS] [CrossRef] [Google Scholar]
  44. Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167 [Google Scholar]
  45. Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 [Google Scholar]
  46. Hayashi, C. 1981, Progr. Theor. Phys. Suppl., 70, 35 [Google Scholar]
  47. Heinzeller, D., Nomura, H., Walsh, C., & Millar, T. J. 2011, ApJ, 731, 115 [NASA ADS] [CrossRef] [Google Scholar]
  48. Helling, C., Woitke, P., Rimmer, P. B., et al. 2014, Life, 4, 142 [Google Scholar]
  49. Hidaka, Y., Takahashi, S., Kawano, H., Suga, M., & Gardiner, W. C. 1982, J. Phys. Chem., 86, 1429 [CrossRef] [Google Scholar]
  50. Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2005, Icarus, 179, 415 [NASA ADS] [CrossRef] [Google Scholar]
  51. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  52. Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, MNRAS, 470, 1750 [Google Scholar]
  53. Javoy, S., Naudet, V., Abid, S., & Paillard, C. 2003, Exp. Thermal Fluid Sci., 27, 371 [CrossRef] [Google Scholar]
  54. Kamp, I., Tilling, I., Woitke, P., Thi, W. F., & Hogerheijde, M. 2010, A&A, 510, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Kamp, I., Thi, W. F., Woitke, P., et al. 2017, A&A, 607, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2016, PASJ, 68, 43 [NASA ADS] [Google Scholar]
  57. Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Kley, W. 1999, MNRAS, 303, 696 [NASA ADS] [CrossRef] [Google Scholar]
  59. Kruijer, T. S., Burkhardt, C., Budde, G., & Kleine, T. 2017, Proc. Natl. Acad. Sci. U.S.A., 114, 6712 [Google Scholar]
  60. Kuramoto, K., & Matsui, T. 1994, J. Geophys. Res.: Planets, 99, 21183 [NASA ADS] [CrossRef] [Google Scholar]
  61. Kuskov, O. L., & Kronrod, V. A. 2005, Icarus, 177, 550 [NASA ADS] [CrossRef] [Google Scholar]
  62. Lari, G., Saillenfest, M., & Fenucci, M. 2020, A&A, 639, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Liu, Y., & Sander, S. 2015, J. Phys. Chem. A, 119, 1 [NASA ADS] [Google Scholar]
  64. Lodders, K. 2003, ApJ, 591, 1220 [Google Scholar]
  65. Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001 [Google Scholar]
  66. Lunine, J. I., & Stevenson, D. J. 1982, Icarus, 52, 14 [NASA ADS] [CrossRef] [Google Scholar]
  67. Machida, M. N., Kokubo, E., Inutsuka, S.-I., & Matsumoto, T. 2008, ApJ, 685, 1220 [NASA ADS] [CrossRef] [Google Scholar]
  68. Madeira, G., Izidoro, A., & Giuliatti Winter, S. M. 2021, MNRAS, 504, 1854 [NASA ADS] [CrossRef] [Google Scholar]
  69. Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007, ApJ, 655, 541 [Google Scholar]
  70. Martin, R. G., & Lubow, S. H. 2011, MNRAS, 413, 1447 [Google Scholar]
  71. McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. McKinnon, W. B. 1997, Icarus, 130, 540 [NASA ADS] [CrossRef] [Google Scholar]
  73. Min, M., Dullemond, C. P., Kama, M., & Dominik, C. 2011, Icarus, 212, 416 [Google Scholar]
  74. Mitchell, T. R., & Stewart, G. R. 2011, AJ, 142, 168 [Google Scholar]
  75. Morbidelli, A., Szulágyi, J., Crida, A., et al. 2014, Icarus, 232, 266 [Google Scholar]
  76. Mousis, O., & Alibert, Y. 2006, A&A, 448, 771 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Nimmo, F., & Korycansky, D. 2012, Icarus, 219, 508 [NASA ADS] [CrossRef] [Google Scholar]
  78. Oba, Y., Watanabe, N., Kouchi, A., Hama, T., & Pirronello, V. 2010, ApJ, 712, L174 [NASA ADS] [CrossRef] [Google Scholar]
  79. Öberg, K. I., & Wordsworth, R. 2019, AJ, 158, 194 [Google Scholar]
  80. Öberg, K. I., Garrod, R. T., van Dishoeck, E. F., & Linnartz, H. 2009, A&A, 504, 891 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Oberg, N., Kamp, I., Cazaux, S., & Rab, C. 2020, A&A, 638, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Ogihara, M., & Ida, S. 2009, ApJ, 699, 824 [Google Scholar]
  83. Pfalzner, S. 2013, A&A, 549, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Pfalzner, S., Steinhausen, M., & Menten, K. 2014, ApJ, 793, L34 [NASA ADS] [CrossRef] [Google Scholar]
  85. Pinilla, P., Birnstiel, T., Benisty, M., et al. 2013, A&A, 554, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Portegies Zwart, S. 2019, A&A, 622, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  88. Portilla-Revelo, B., Kamp, I., Rab, C., et al. 2022, A&A, 658, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Pringle, J. E. 1981, ARA&A, 19, 137 [NASA ADS] [CrossRef] [Google Scholar]
  90. Quillen, A. C., & Trilling, D. E. 1998, ApJ, 508, 707 [Google Scholar]
  91. Rab, C., Kamp, I., Ginski, C., et al. 2019, A&A, 624, A16 [EDP Sciences] [Google Scholar]
  92. Rafikov, R. R. 2017, ApJ, 837, 163 [NASA ADS] [CrossRef] [Google Scholar]
  93. Rice, W. K. M., Armitage, P. J., Wood, K., & Lodato, G. 2006, MNRAS, 373, 1619 [Google Scholar]
  94. Rivier, G., Crida, A., Morbidelli, A., & Brouet, Y. 2012, A&A, 548, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Röllig, M., Abel, N. P., Bell, T., et al. 2007, A&A, 467, 187 [Google Scholar]
  96. Ronnet, T., & Johansen, A. 2020, A&A, 633, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Ronnet, T., Mousis, O., & Vernazza, P. 2017, ApJ, 845, 92 [NASA ADS] [CrossRef] [Google Scholar]
  98. Ruaud, M., Wakelam, V., & Hersant, F. 2016, MNRAS, 459, 3756 [Google Scholar]
  99. Ruffle, D. P., Hartquist, T. W., Caselli, P., & Williams, D. A. 1999, MNRAS, 306, 691 [Google Scholar]
  100. Sakai, N., Saruwatari, O., Sakai, T., Takano, S., & Yamamoto, S. 2010, A&A, 512, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Sasaki, T., Stewart, G. R., & Ida, S. 2010, ApJ, 714, 1052 [NASA ADS] [CrossRef] [Google Scholar]
  102. Schubert, G., Anderson, J. D., Spohn, T., & McKinnon, W. B. 2004, Interior composition, structure and dynamics of the Galilean satellites, eds. F. Bagenal, T. E. Dowling, & W. B. McKinnon, 281 [Google Scholar]
  103. Schulik, M., Johansen, A., Bitsch, B., Lega, E., & Lambrechts, M. 2020, A&A, 642, A187 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  104. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  105. Shibaike, Y., Okuzumi, S., Sasaki, T., & Ida, S. 2017, ApJ, 846, 81 [NASA ADS] [CrossRef] [Google Scholar]
  106. Shibaike, Y., Ormel, C. W., Ida, S., Okuzumi, S., & Sasaki, T. 2019, ApJ, 885, 79 [NASA ADS] [CrossRef] [Google Scholar]
  107. Showman, A. P., & Malhotra, R. 1997, Icarus, 127, 93 [NASA ADS] [CrossRef] [Google Scholar]
  108. Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593 [Google Scholar]
  109. Smith, I. W. M., Herbst, E., & Chang, Q. 2004, MNRAS, 350, 323 [Google Scholar]
  110. Sohl, F., Spohn, T., Breuer, D., & Nagel, K. 2002, Icarus, 157, 104 [Google Scholar]
  111. Spiegel, D. S., & Burrows, A. 2012, ApJ, 745, 174 [Google Scholar]
  112. Suutarinen, A., Geppert, W. D., Harju, J., et al. 2011, A&A, 531, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  113. Szulágyi, J. 2017, ApJ, 842, 103 [Google Scholar]
  114. Szulágyi, J., & Mordasini, C. 2017, MNRAS, 465, L64 [CrossRef] [Google Scholar]
  115. Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65 [CrossRef] [Google Scholar]
  116. Szulágyi, J., Masset, F., Lega, E., et al. 2016, MNRAS, 460, 2853 [Google Scholar]
  117. Szulágyi, J., Binkert, F., & Surville, C. 2022, ApJ, 924, 1 [CrossRef] [Google Scholar]
  118. Takata, T., & Stevenson, D. J. 1996, Icarus, 123, 404 [NASA ADS] [CrossRef] [Google Scholar]
  119. Tanigawa, T., Ohtsuki, K., & Machida, M. N. 2012, ApJ, 747, 47 [NASA ADS] [CrossRef] [Google Scholar]
  120. Teague, R., Bae, J., & Bergin, E. A. 2019, Nature, 574, 378 [NASA ADS] [CrossRef] [Google Scholar]
  121. Thanathibodee, T., Calvet, N., Bae, J., Muzerolle, J., & Hernández, R. F. 2019, ApJ, 885, 94 [NASA ADS] [CrossRef] [Google Scholar]
  122. Thi, W. F., Woitke, P., & Kamp, I. 2011, MNRAS, 412, 711 [Google Scholar]
  123. Thi, W. F., Hocuk, S., Kamp, I., et al. 2020a, A&A, 634, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  124. Thi, W. F., Hocuk, S., Kamp, I., et al. 2020b, A&A, 635, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  125. Tielens, A. 2021, Molecular Astrophysics (Cambridge University Press) [CrossRef] [Google Scholar]
  126. Tittemore, W. C. 1990, Science, 250, 263 [NASA ADS] [CrossRef] [Google Scholar]
  127. Tsang, W., & Hampson, R. F. 1986, J. Phys. Chem. Ref. Data, 15, 1087 [NASA ADS] [CrossRef] [Google Scholar]
  128. Turner, N. J., Choukroun, M., Castillo-Rogez, J., & Bryden, G. 2012, ApJ, 748, 92 [NASA ADS] [CrossRef] [Google Scholar]
  129. Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2022, ApJ, 930, 11 [NASA ADS] [CrossRef] [Google Scholar]
  130. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  131. Wagner, K., Follete, K. B., Close, L. M., et al. 2018, ApJ, 863, L8 [Google Scholar]
  132. Wakelam, V., Bron, E., Cazaux, S., et al. 2017, Mol. Astrophys., 9, 1 [Google Scholar]
  133. Walsh, C., Harada, N., Herbst, E., & Millar, T. J. 2009, ApJ, 700, 752 [NASA ADS] [CrossRef] [Google Scholar]
  134. Weidenschilling, S. J. 1977, MNRAS, 180, 57 [Google Scholar]
  135. Weiss, B. P., & Bottke, W. F. 2021, AGU Adv., 2, e00376 [NASA ADS] [CrossRef] [Google Scholar]
  136. Woitke, P., Dominik, C., & Sedlmayr, E. 1993, A&A, 274, 451 [NASA ADS] [Google Scholar]
  137. Woitke, P., Kamp, I., & Thi, W.-F. 2009, A&A, 501, 383 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  138. Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  139. Woodall, J., Agúndez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197 [CrossRef] [EDP Sciences] [Google Scholar]
  140. Yoder, C. F. 1979, Nature, 279, 767 [NASA ADS] [CrossRef] [Google Scholar]
  141. Zangi, R., & Mark, A. E. 2003, Phys. Rev. Lett., 91, 025502 [NASA ADS] [CrossRef] [Google Scholar]
  142. Zhou, Y., Bowler, B. P., Wagner, K. R., et al. 2021, AJ, 161, 244 [NASA ADS] [CrossRef] [Google Scholar]
  143. Zhu, Z. 2015, ApJ, 799, 16 [Google Scholar]
  144. Zhu, Z., Nelson, R. P., Dong, R., Espaillat, C., & Hartmann, L. 2012, ApJ, 755, 6 [Google Scholar]
  145. Zhu, Z., Andrews, S. M., & Isella, A. 2018, MNRAS, 479, 1850 [Google Scholar]

All Tables

Table 1

Parameters for the solar circumstellar disk.

Table 2

Conditions at the gap edge “inheritance” point of the circumstellar disk at heliocentric distance 8.2 au, altitude 0.5 au (1 pressure scale height) above the midplane.

Table 3

Variation of parameters for the circumplanetary disk model grid.

Table 4

Parameters common to the CPD models.

Table B.1

Adsorption energies of the most prevalent molecular ices found in our model CPDs.

Table G.1

Molecular cloud parameters

All Figures

thumbnail Fig. 1

2D dust temperature distribution around the circumstellar disk gap in the range 40–120 K (top). Solid, dashed, and dotted black contours indicate the relative UV field strength χ (see Eq. (12)). The white circle with black border represents the position of Jupiter and the horizontal bar indicates the physical extent of the CPD. The white cross with black border represents the location from which the inherited chemistry is extracted at one scale height z = H. Hydrogen nuclei column density (blue line) with unperturbed circumstellar disk density profile for comparison (light blue dashed line; bottom).

In the text
thumbnail Fig. 2

Ice abundance in the circumstellar disk gap normalized to the initial abundance of ice for various stellar luminosities as function of time after the onset of UV exposure. The width of each track represents the range of ice sublimation rates corresponding to variable conditions across the gap region (z = 0–0.5 au, r = 5.2–8.2 au).

In the text
thumbnail Fig. 3

Surface density (red) and midplane FUV field strength in units of the Draine field χ (blue) of the high (10−7 M) and low (10−8 M) mass CPDs. The four circles indicate the present-day radial position of the Galilean satellites and their position on the ordinate is arbitrary. The striped gray region on the left indicates an empty inner magnetic cavity and the light gray region on the right indicates the gravitationally unstable zone outside of RH/3 (where RH = 0.35 au).

In the text
thumbnail Fig. 4

Radial midplane dust temperature profiles of the four reference CPDs (d/g = 10−33), labelled by the model id’s in Table 3. The four circles indicate only the present-day radial position of the Galilean satellites. The striped gray area on the left indicates the inner cavity. The light gray area on the right indicates the gravitationally unstable zone outside of RH/3 (where RH = 0.35 au).

In the text
thumbnail Fig. 5

Rate of water ice formation as a function of time at different radii beyond the snowline in the high-mass, high-viscosity CPD (7–10), to illustrate the two distinct stages of water ice formation. All values are normalized to the maximum formation rate at 10−5 yr, r = 0.08 RH. The times where the decline in the abundance of free O and free H end are indicated by the gray vertical lines. The starting time t0 is defined as when the infalling gas is shock-heated to an atomic and ionized state.

In the text
thumbnail Fig. 6

Midplane ice fraction at t = tvisc for the CPDs with d/g = 10−33. The four circles indicate the present day radial location of the Galilean satellites and the uncertainty bars represent their possible range of ice fractions. The thin dotted line indicates the initial radial ice fraction in the inheritance cases. The filled gray region on the right indicates the gravitationally unstable zone outside of RH/3.

In the text
thumbnail Fig. 7

Radial drift velocity of particles (upper panel) and radial velocity of the gas (lower panel) in the high-mass, high-viscosity CPD (7–10). The position of the Galilean satellites on the ordinate is arbitrary. The gray region on the right indicates the tidally unstable region beyond RH/3. The striped region on the left indicates the inner cavity. Solid lines indicate the case of a decretion disk where gas falls onto the CPD at the centrifugal radius (0.03 RH). Dashed lines indicate the grain drift and gas radial velocity in the case of a pure accretion disk.

In the text
thumbnail Fig. 8

Inwards drift timescale for dust grains (colored lines) relative to the timescales of desorption and ice formation (black dashed lines) in the reference CPDs with d/g = 10−33. From left to right the colored regions correspond to areas where all ices are eventually thermally desorbed (pink), the iceline (light yellow), where grains become icier (light blue) and where all ices are eventually photodesorbed (dark gray). The red line corresponds to the minimum grain size which is not trapped in the CPD (if there is a sign change in the gas radial velocity) and the blue line to the maximum for which the timescales of desorption and ice formation have been verified. The vertical dotted line indicates the approximate outer region of gravitational stability in the CPD.

In the text
thumbnail Fig. A.1

Volume and areas for the spatial cell around point (i, j). The mean intensities J and Rosseland opacities κR are available on the grid points (i, j) marked with red dots.

In the text
thumbnail Fig. A.2

Benchmark test for a viscous circumstellar disk with a disk mass of 0.01 M and a mass accertion rate of Macc = 10−8M/yr, see text. The green lines are temperature cuts at selected radii calculated by MCMax, the small wiggles are due to the Monte-Carlo noise. The black dots show the temperatures calculated by ProDiMo.

In the text
thumbnail Fig. A.3

Midplane temperature Tdust(r, 0) as function of the log-distance from the inner rim. Again, the green line is the temperature profile calculated by MCMax, and the small black diamonds show the temperature values calculated by ProDiMo.

In the text
thumbnail Fig. A.4

Maximum relative temperature change between iterations as function of RT iteration number.

In the text
thumbnail Fig. C.1

Midplane ice mass fraction for a variety of surface density powerlaw exponents ϵ for the (7–11) reset CPD with global d/g = 10−33.

In the text
thumbnail Fig. D.1

Midplane ice mass fraction for a variety of background temperatures Tback for the (7–11) reset CPD and global d/g = 10−3.3. The four empty circles with error bars represent the Galilean satellites at their present day locations and composition. The gray striped region on the left represents the inner cavity and the light gray shaded region on the right represents the gravitationally unstable zone.

In the text
thumbnail Fig. E.1

Evolution of the water ice abundance relative to the total number of hydrogen nuclei at the midplane (blue) and at 5 scale heights above the CPD midplane (red), and in the case of the gas parcel which drifts from z = 5H downwards to the midplane over a period of 100 yr (grey). The period in which the altitude of the gas parcel was iteratively lowered towards the midplane is highlighted in light green.

In the text
thumbnail Fig. G.1

Ratio of the column density of several common species relative to H2 in the 0D molecular cloud model (circles with error bars) to observational values for TMC-1 (squares with black border). Errorbars on model values represent a factor 10x uncertainty for illustrative purposes. TMC-1 abundances are taken from (Suutarinen etal. 2011) (OH,CH), (Sakai et al. 2010) (C2H), (Smith et al. 2004; Walsh et al. 2009) otherwise.

In the text
thumbnail Fig. H.1

Midplane radial ice mass fraction fice of the dust-to-gas ratio parameter exploration for the chemically reset (left column) and chemically inherited (right column) reference CPDs. The black dotted, dashed, and solid lines indicate the radial position of Europa, Ganymede, and Callisto. The white dotted, dashed, and solid lines indicate where fice is equivalent to the estimated ice mass fraction of Europa (0.05), Ganymede (0.48), and Callisto (0.52).

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.