Circumplanetary disk ices

A&A 667, A95 (2022) https://doi.org/10.1051/0004-6361/202244092 © N. Oberg et al. 2022 Astronomy & Astrophysics Circumplanetary disk ices I. Ice formation vs. viscous evolution and grain drift N. Oberg1,2 , I. Kamp1 , S. Cazaux2,3 , P. Woitke4,5 , and W. F. Thi6 1 2 3 4 5 6 Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands e-mail: oberg@astro.rug.nl Faculty of Aerospace Engineering, Delft University of Technology, Delft, The Netherlands University of Leiden, PO Box 9513, 2300 RA, Leiden, The Netherlands Space Research Institute, Austrian Academy of Sciences, Schmiedlstrasse 6, 8042 Graz, Austria Centre for Exoplanet Science, University of St. Andrews, North Haugh, St. Andrews, KY16 9SS, UK Max-Planck-Institut für extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching, Germany Received 23 May 2022 / Accepted 13 August 2022 ABSTRACT Context. The large icy moons of Jupiter formed in a circumplanetary disk (CPD). CPDs are fed by vertically infalling circumstellar gas and dust which may be shock-heated upon accretion. Accreted material is then either incorporated into moons, falls into the planet, or is lost beyond the disk edge on relatively short timescales. If ices are sublimated during accretion onto the CPD we know there must be sufficient time for them to recondense or moons such as Ganymede or Callisto could not form. The chemical timescale to form sufficiently icy solids places a novel constraint on the dynamical behaviour and properties of CPDs. Aims. We aim to explore the process of ice formation in CPDs to constrain which disk properties (such as the mass, viscosity, and dust-to-gas ratio) are consistent with the formation of an icy moon system. Methods. We use the radiation thermochemical code P RO D I M O (Protoplanetary Disk Model) to analyze how the radial ice abundance evolves in CPDs. We consider different initial chemical conditions of the disk to explore the consequences of infalling material being inherited from the circumstellar disk or being reset to atomic conditions by shock-heating. We contrast the timescales of ice formation with disk viscous timescales and radial dust drift. Results. We have derived the radial ice abundance and rate of ice formation in a small grid of model CPDs. Water ice can form very efficiently in the CPD from initially atomic conditions, as a significant fraction is efficiently re-deposited on dust grains within <1 yr. Radial grain drift timescales are in general longer than those of ice formation on grains. Icy grains of size a < 3 mm retain their icy mantles while crossing an optically thin circumstellar disk gap at 5 au for L∗ < 10 L⊙ . Conclusions. Three-body reactions play an important role in water formation in the dense midplane condition of CPDs. The CPD midplane must be depleted in dust relative to the circumstellar disk by a factor 10–50 to produce solids with the ice to rock ratio of the icy Galilean satellites. The CPD snowline is not erased by radial grain drift, which is consistent with the compositional gradient of the Galilean satellites being primordial. Key words. planets and satellites: formation – planets and satellites: composition – accretion, accretion disks – protoplanetary disks – planets and satellites: individual: Jupiter – methods: numerical 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 A95, page 1 of 19 Open 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. A&A 667, A95 (2022) 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 H2 O 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 A95, page 2 of 19 the applicable chemical evolutionary timescales. Hereafter we refer to this as the viscous timescale, defined as MCPD tvisc = , (1) Ṁ where MCPD is the mass of the CPD and Ṁ is the infall rate of circumstellar material onto the CPD. The relatively short tvisc 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 P RO D I M O1 (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. P RO D I M O 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 threebody 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 P RO D I M O 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.). 1 https://www.astro.rug.nl/∼prodimo/ N. Oberg et al.: Circumplanetary disk ices. I. 2.1.1. Ice formation 2.1.2. Properties of the disk models 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 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.46 Myr 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, −1 Rads = 4πa2 vth i i nd S i s , (2) where 4πa2 is a dust grain surface area, a is the grain radius, vth i is the thermal speed (kB T g /2πmi )1/2 , kB is the Boltzmann constant, T g is the gas temperature, mi is the mass of the gas-phase species, nd is the dust number density, and S i is the sticking coefficient. The ice adsorption rate coefficients are then dn#,i −3 −1 = Rads i ni cm s , dt (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 . The Arrhenius formulation for the rate of thermal desorption is then b Rides,th = v0,i e−Ei /kB Td in s −1 , where the pre-exponential (frequency) factor v0,i is s 2Nsurf Eib v0,i = . π2 mi (4) (5) Nsurf is the density of surface binding sites 1.5 × 1015 cm−2 (Hasegawa et al. 1992), and T d 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 i is des,ph Ri = πa2 nd Yi χFDraine s−1 , nact i (6) where Yi is the photodesorption yield, nact i is the concentration of the species in the active layers, nact = ni i if n#, tot 6 4πNsurf a2 nd  = ni Nact /Nlayer  if n#, tot > 4πNsurf a2 nd . (7) The number ofPphysisorbed layers Nlayer = n#, tot /(4πNsurf a2 nd ) and n#, tot = i n#,i is the total number of the physisorbed species, 4πNsurf a2 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 is calculated according to the method of Hasegawa & Herbst (1993). The total desorption rate Rdes i is the sum of the thermal, photo- and cosmic ray induced desdes,ph orption rates Rides,th + Ri + Rides,CR . The desorption rate of the physisorbed species i is then dni act −3 −1 = Rdes i ni cm s , dt (8) fice = mice , mrock + mice (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 αsettle = 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 A95, page 3 of 19 A&A 667, A95 (2022) Table 1. Parameters for the solar circumstellar disk. Parameter Symbol Value Unit Stellar mass Stellar Luminosity Effective temperature UV luminosity X-ray luminosity M∗ L∗ T eff LUV,∗ LX 1.0 0.84 4395 0.01 1030 M⊙ L⊙ K L⊙ erg s−1 Disk mass Disk inner radius Disk outer radius Column density power ind. Reference scale height Flaring index Mdisk Rin Rout ǫ H10au β 0.001 0.1 100 1.5 1 1.15 M⊙ au au – au – Minimum dust size Maximum dust size Dust size power law index Dust-to-gas ratio amin amax apow d/g 0.05 3000 3.5 10−2 ➭m ➭m Dust composition: Mg0.7 Fe0.3 SiO3 Amorphous carbon Vacuum 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. Parameter Hydrogen density Optical extinction Dust-to-gas ratio Dust temperature Gas temperature Ice mass fraction Symbol Value Unit nH,in AV,in d/g in T d,in T g,in fice,in 1010 1.01 10−3.23 57.0 57.3 0.48 cm−3 – – K K – – – 60% 15% 25% Notes. Stellar temperature and luminosity are selected from the premain sequence stellar evolutionary tracks of Siess et al. (2000) for t = 4 Myr. Stellar UV and X-ray luminosities for a representative Class II T Tauri star are adopted from Woitke et al. (2016). 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). 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). of time prior to being accreted onto the CPD. We considered whether ices on grains can survive against thermal- or photodesorption 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 A95, page 4 of 19 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 4 Myr (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 M⊙ 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 N. Oberg et al.: Circumplanetary disk ices. I. Table 3. Variation of parameters for the circumplanetary disk model grid. Model id (7–10) (7–11) (8–11) (8–12) M (M⊙ ) Ṁ (M⊙ yr−1 ) tvisc (yr) 10 10−7 10−8 10−8 10 10−11 10−11 10−12 10 104 103 104 −7 −10 3 α 10−2.7 10−3.6 10−2.7 10−3.6 Notes. Model id’s are of the format (A–B) where the CPD mass is 10−A M⊙ and where the accretion rate onto the CPD is 10−B M⊙ yr−1 . Table 4. Parameters common to the CPD models. Parameter Symbol Value Unit Mp Lp T eff,p LUV,p χ T back 1.0 10−5 1000 0.01 3 × 103 50 MJ L⊙ K Lp – K Disk inner radius Taper radius Disk outer radius Column density power ind. Flaring index Reference scale height Rin,CPD Rtap,CPD Rout,CPD ǫCPD βCPD H0.1 au 0.0015 0.12 0.35 1.0 1.15 0.01 au au au – – au Maximum dust size Dust-to-gas ratio amax,CPD d/g 3000 10−3.3 µm – Planet mass Planetary luminosity Effective temperature UV luminosity Background UV field Background temperature Notes. Parameters which are not common to the CPD models are listed in Table 3. Where not specified the CPD parameters are identical to the circumstellar disk parameters listed in Table 1. Jovian-mass planet these represent planet-disk mass ratios of ∼10−4 and 10−5 , respectively. The CPDs are thus of the “gasstarved” 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 (x–y) where x is related to the CPD mass by MCPD = 10−x M⊙ and y to the mass infall rate (and thus viscosity) by ṀCPD = 10−y M⊙ yr−1 . The list of parameters which are common between the reference CPDs can be found in Table 4. 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 CPD viscosity. weak to explain observations (Shakura & Sunyaev 1973; Pringle 1981). The efficiency of angular momentum transport is parameterized by the dimensionless α-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 P RO D I M O 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 Ṁ = 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 Ṁ = 10−10.8 −10−10.3 M⊙ yr−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 Ṁ = 10−10.9 −10−10.8 M⊙ yr−1 (Zhou et al. 2021). Based on these observational and theoretical constraints we adopt Ṁ = 10−10 M⊙ yr−1 (with a heating rate corresponding to α ≈ 10−2.7 ) and Ṁ = 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 lowmass CPD we adopt Ṁ = 10−11 M⊙ yr−1 and Ṁ = 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) = q 3GMp Ṁ (1 − Rp /r) erg cm−2 s−1 , 8πr3 (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) R ρd (r, z) ρd (r, z′ )dz′ erg cm−3 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 midplane 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 A95, page 5 of 19 A&A 667, A95 (2022) 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). Fig. 4. Radial midplane dust temperature profiles of the four reference CPDs (d/g = 10−3.3 ), 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). 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, bg Iν = χ · 1.71 · Wdil Bν (20 000 K) + Bν (2.7 K), (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 A95, page 6 of 19 provides a minimum χ ∼ 3000 at 5 au (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 (Draż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. 2.2. Dust grain drift within the CPD In P RO D I M O, the chemistry is solved for each grid cell without accounting for radial dust or gas transport. Instead we used properties of the P RO D I M O 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 subKeplerianity can be <1% while in a CPD it could be as high as 20–80% due to significant gas pressure support (Armitage 2007; Draż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 N. Oberg et al.: Circumplanetary disk ices. I. drift velocity vr,d can be approximated as vr,d = vr,g T s−1 − ηvK T s + T s−1 , (13) where vr,g is the radial velocity of the gas, T s is the dimensionless stopping time of a grain, v  k T s = ts , (14) r where ts is the stopping time, vk is the Keplerian orbital velocity, and r is the radius in the CPD. The stopping time is ! ! ρgrain a ts = , (15) ρgas vth 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  2  c  (16) η = n  2s  , vk 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 m s−1 . For our low viscosity case vr,gas ∼ 0.1 m s−1 . 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. 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. 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), 3. Results H + O + M → OH + M, (17) In our reference model grid (Table 3), we considered the case of a higher mass, optically thick CPD (10−7 M⊙ ), and lower mass, 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 A95, page 7 of 19 A&A 667, A95 (2022) association of the OH with free hydrogen; OH + H → H2 O + photon, 3.2. The midplane ice mass fraction (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 photoreactions. 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 + NO → N2 + H2 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 threebody 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 H3 O+ + e− → H2 O + H, (20) and in the outermost part of the disk where water ice is stable via H3 O+ + H2 CO → H3 CO+ + H2 O. (21) In the high-mass CPDs water ice formation is typically still ongoing by the end of the viscous timescale, and so the midplane 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. A95, page 8 of 19 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 grainsurface 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 ǫ 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 H2 O formation occur simultaneously. The analytic rate of Cazaux & Tielens (2002) presupposes that chemisorbed H plays a role in N. Oberg et al.: Circumplanetary disk ices. I. (a) High-Mass High-Viscosity (7-10) (b) High-Mass Low-Viscosity (7-11) (c) Low-Mass High-Viscosity (8-11) (d) Low-Mass Low-Viscosity (8-12) Fig. 6. Midplane ice fraction at t = tvisc for the CPDs with d/g = 10−3.3 . 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. 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 H2 O 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 O2 H in the surface chemistry network. Gas-phase three-body reactions with O2 produce O2 H which in turn lead to early OH formation. The gas-phase O2 reservoir is thus depleted and efficiently converted via OH into H2 O via the three-body reaction O2 + H + M → O2 H + 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 O2 H + H → OH + OH, (23) OH + H → H2 O + photon. (24) In the outer disk, the ice mass fraction can be enhanced relative to the gas-phase chemical network as the freezingout 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 A95, page 9 of 19 A&A 667, A95 (2022) 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. 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. 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.). 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 highviscosity 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 A95, page 10 of 19 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, gasphase O2 H 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 gasdriven migration (Canup & Ward 2002; Madeira et al. 2021) and long-term outwards tidally-driven migration (Yoder 1979; N. Oberg et al.: Circumplanetary disk ices. I. (a) High-Mass High-Viscosity (7-10) (b) High-Mass Low-Viscosity (7-11) (c) Low-Mass High-Viscosity (8-11) (d) Low-Mass Low-Viscosity (8-12) 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−3.3 . 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. 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 dustto-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 A95, page 11 of 19 A&A 667, A95 (2022) 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. 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 gasfree magnetically-truncated inner cavity (Takata & Stevenson 1996; Batygin 2018; Shibaike et al. 2019). Additional protosatellites 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. A95, page 12 of 19 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. References Adams, F. C. 2010, ARA&A, 48, 47 Alibert, Y., Mousis, O., & Benz, W. 2005, A&A, 439, 1205 Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21 Aota, T., Inoue, T., & Aikawa, Y. 2015, ApJ, 799, 141 Aoyama, Y., Ikoma, M., & Tanigawa, T. 2018, ApJ, 866, 84 Armitage, P. J. 2007, ArXiv e-prints [arXiv:astro-ph/0701485] Armitage, P. J. 2010, Astrophysics of Planet Formation (Cambridge University Press) Atkinson, R., Baulch, D. L., Cox, R. A., et al. 2004, Atmos. Chem. Phys., 4, 1461 Auer, L. 1984, in Numerical Radiative Transfer, ed. W. Kalkhofen (Cambridge: Cambridge University Press), 101 Avramenko, L. I., & Krasnen’kov, V. M. 1966, Bull. Acad. Sci. USSR, Div. Chem. Sci., 15, 394 Ayliffe, B. A., & Bate, M. R. 2009, MNRAS, 397, 657 Balbus, S. A. 2003, ARA&A, 41, 555 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 Batygin, K. 2018, AJ, 155, 178 Batygin, K., & Morbidelli, A. 2020, ApJ, 894, 143 Baulch, D. L. 1972, Evaluated kinetic data for high temperature reactions (Butterworths) Bierson, C. J., & Nimmo, F. 2020, ApJ, 897, L43 Bierson, C. J., & Steinbrügge, G. 2021, Planet. Sci. J., 2, 89 Binkert, F., Szulágyi, J., & Birnstiel, T. 2021, MNRAS, 506, 5969 Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 Boschman, L., Reitsma, G., Cazaux, S., et al. 2012, ApJ, 761, L33 Bosman, A. D., Cridland, A. J., & Miguel, Y. 2019, A&A, 632, A11 Brown, W. A., & Bolina, A. S. 2006, MNRAS, 374, 1006 Canup, R. M., & Ward, W. R. 2002, AJ, 124, 3404 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 Cazaux, S., & Tielens, A. G. G. M. 2002, ApJ, 575, L29 Christiaens, V., Cantalloube, F., Casassus, S., et al. 2019, ApJ, 877, L33 Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133 D’Alessio, P., Cantö, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411 D’Angelo, G., Henning, T., & Kley, W. 2002, A&A, 385, 647 Day, M., Thompson, K., & Dixon-Lewis, G. 1973, Symp. (Int.) Combust., 14, 47 de Pater, I., & Lissauer, J. J. 2010, Planet. Sci. (Cambridge, UK: Cambridge University Press) N. Oberg et al.: Circumplanetary disk ices. I. Draine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269 Drażkowska, ˛ J., & Szulágyi, J. 2018, ApJ, 866, 142 Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237 Dwyer, C. A., Nimmo, F., Ogihara, M., & Ida, S. 2013, Icarus, 225, 390 Fujii, Y. I., Okuzumi, S., Tanigawa, T., & Inutsuka, S.-i. 2014, ApJ, 785, 101 Fulle, D., Hamann, H. F., Hippler, H., & Troe, J. 1996, J. Chem. Phys., 105, 983 Garrod, R. T., & Herbst, E. 2006, A&A, 457, 927 Gressel, O., Nelson, R. P., Turner, N. J., & Ziegler, U. 2013, ApJ, 779, 59 Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 261, 83 Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167 Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 Hayashi, C. 1981, Progr. Theor. Phys. Suppl., 70, 35 Heinzeller, D., Nomura, H., Walsh, C., & Millar, T. J. 2011, ApJ, 731, 115 Helling, C., Woitke, P., Rimmer, P. B., et al. 2014, Life, 4, 142 Hidaka, Y., Takahashi, S., Kawano, H., Suga, M., & Gardiner, W. C. 1982, J. Phys. Chem., 86, 1429 Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2005, Icarus, 179, 415 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, MNRAS, 470, 1750 Javoy, S., Naudet, V., Abid, S., & Paillard, C. 2003, Exp. Thermal Fluid Sci., 27, 371 Kamp, I., Tilling, I., Woitke, P., Thi, W. F., & Hogerheijde, M. 2010, A&A, 510, A18 Kamp, I., Thi, W. F., Woitke, P., et al. 2017, A&A, 607, A41 Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2016, PASJ, 68, 43 Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44 Kley, W. 1999, MNRAS, 303, 696 Kruijer, T. S., Burkhardt, C., Budde, G., & Kleine, T. 2017, Proc. Natl. Acad. Sci. U.S.A., 114, 6712 Kuramoto, K., & Matsui, T. 1994, J. Geophys. Res.: Planets, 99, 21183 Kuskov, O. L., & Kronrod, V. A. 2005, Icarus, 177, 550 Lari, G., Saillenfest, M., & Fenucci, M. 2020, A&A, 639, A40 Liu, Y., & Sander, S. 2015, J. Phys. Chem. A, 119, 1 Lodders, K. 2003, ApJ, 591, 1220 Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001 Lunine, J. I., & Stevenson, D. J. 1982, Icarus, 52, 14 Machida, M. N., Kokubo, E., Inutsuka, S.-i., & Matsumoto, T. 2008, ApJ, 685, 1220 Madeira, G., Izidoro, A., & Giuliatti Winter, S. M. 2021, MNRAS, 504, 1854 Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007, ApJ, 655, 541 Martin, R. G., & Lubow, S. H. 2011, MNRAS, 413, 1447 McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36 McKinnon, W. B. 1997, Icarus, 130, 540 Min, M., Dullemond, C. P., Kama, M., & Dominik, C. 2011, Icarus, 212, 416 Mitchell, T. R., & Stewart, G. R. 2011, AJ, 142, 168 Morbidelli, A., Szulágyi, J., Crida, A., et al. 2014, Icarus, 232, 266 Mousis, O., & Alibert, Y. 2006, A&A, 448, 771 Nimmo, F., & Korycansky, D. 2012, Icarus, 219, 508 Oba, Y., Watanabe, N., Kouchi, A., Hama, T., & Pirronello, V. 2010, ApJ, 712, L174 Öberg, K. I., & Wordsworth, R. 2019, AJ, 158, 194 Öberg, K. I., Garrod, R. T., van Dishoeck, E. F., & Linnartz, H. 2009, A&A, 504, 891 Oberg, N., Kamp, I., Cazaux, S., & Rab, C. 2020, A&A, 638, A135 Ogihara, M., & Ida, S. 2009, ApJ, 699, 824 Pfalzner, S. 2013, A&A, 549, A82 Pfalzner, S., Steinhausen, M., & Menten, K. 2014, ApJ, 793, L34 Pinilla, P., Birnstiel, T., Benisty, M., et al. 2013, A&A, 554, A95 Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967 Portegies Zwart, S. 2019, A&A, 622, A69 Portilla-Revelo, B., Kamp, I., Rab, C., et al. 2022, A&A, 658, A89 Pringle, J. E. 1981, ARA&A, 19, 137 Quillen, A. C., & Trilling, D. E. 1998, ApJ, 508, 707 Rab, C., Kamp, I., Ginski, C., et al. 2019, A&A, 624, A16 Rafikov, R. R. 2017, ApJ, 837, 163 Rice, W. K. M., Armitage, P. J., Wood, K., & Lodato, G. 2006, MNRAS, 373, 1619 Rivier, G., Crida, A., Morbidelli, A., & Brouet, Y. 2012, A&A, 548, A116 Röllig, M., Abel, N. P., Bell, T., et al. 2007, A&A, 467, 187 Ronnet, T., & Johansen, A. 2020, A&A, 633, A93 Ronnet, T., Mousis, O., & Vernazza, P. 2017, ApJ, 845, 92 Ruaud, M., Wakelam, V., & Hersant, F. 2016, MNRAS, 459, 3756 Ruffle, D. P., Hartquist, T. W., Caselli, P., & Williams, D. A. 1999, MNRAS, 306, 691 Sakai, N., Saruwatari, O., Sakai, T., Takano, S., & Yamamoto, S. 2010, A&A, 512, A31 Sasaki, T., Stewart, G. R., & Ida, S. 2010, ApJ, 714, 1052 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 Schulik, M., Johansen, A., Bitsch, B., Lega, E., & Lambrechts, M. 2020, A&A, 642, A187 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 Shibaike, Y., Okuzumi, S., Sasaki, T., & Ida, S. 2017, ApJ, 846, 81 Shibaike, Y., Ormel, C. W., Ida, S., Okuzumi, S., & Sasaki, T. 2019, ApJ, 885, 79 Showman, A. P., & Malhotra, R. 1997, Icarus, 127, 93 Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593 Smith, I. W. M., Herbst, E., & Chang, Q. 2004, MNRAS, 350, 323 Sohl, F., Spohn, T., Breuer, D., & Nagel, K. 2002, Icarus, 157, 104 Spiegel, D. S., & Burrows, A. 2012, ApJ, 745, 174 Suutarinen, A., Geppert, W. D., Harju, J., et al. 2011, A&A, 531, A121 Szulágyi, J. 2017, ApJ, 842, 103 Szulágyi, J., & Mordasini, C. 2017, MNRAS, 465, L64 Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65 Szulágyi, J., Masset, F., Lega, E., et al. 2016, MNRAS, 460, 2853 Szulágyi, J., Binkert, F., & Surville, C. 2022, ApJ, 924, 1 Takata, T., & Stevenson, D. J. 1996, Icarus, 123, 404 Tanigawa, T., Ohtsuki, K., & Machida, M. N. 2012, ApJ, 747, 47 Teague, R., Bae, J., & Bergin, E. A. 2019, Nature, 574, 378 Thanathibodee, T., Calvet, N., Bae, J., Muzerolle, J., & Hernández, R. F. 2019, ApJ, 885, 94 Thi, W. F., Woitke, P., & Kamp, I. 2011, MNRAS, 412, 711 Thi, W. F., Hocuk, S., Kamp, I., et al. 2020a, A&A, 634, A42 Thi, W. F., Hocuk, S., Kamp, I., et al. 2020b, A&A, 635, A16 Tielens, A. 2021, Molecular Astrophysics (Cambridge University Press) Tittemore, W. C. 1990, Science, 250, 263 Tsang, W., & Hampson, R. F. 1986, J. Phys. Chem. Ref. Data, 15, 1087 Turner, N. J., Choukroun, M., Castillo-Rogez, J., & Bryden, G. 2012, ApJ, 748, 92 Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2022, ApJ, 930, 11 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 Wagner, K., Follete, K. B., Close, L. M., et al. 2018, ApJ, 863, L8 Wakelam, V., Bron, E., Cazaux, S., et al. 2017, Mol. Astrophys., 9, 1 Walsh, C., Harada, N., Herbst, E., & Millar, T. J. 2009, ApJ, 700, 752 Weidenschilling, S. J. 1977, MNRAS, 180, 57 Weiss, B. P., & Bottke, W. F. 2021, AGU Adv., 2, e00376 Woitke, P., Dominik, C., & Sedlmayr, E. 1993, A&A, 274, 451 Woitke, P., Kamp, I., & Thi, W.-F. 2009, A&A, 501, 383 Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103 Woodall, J., Agúndez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197 Yoder, C. F. 1979, Nature, 279, 767 Zangi, R., & Mark, A. E. 2003, Phys. Rev. Lett., 91, 025502 Zhou, Y., Bowler, B. P., Wagner, K. R., et al. 2021, AJ, 161, 244 Zhu, Z. 2015, ApJ, 799, 16 Zhu, Z., Nelson, R. P., Dong, R., Espaillat, C., & Hartmann, L. 2012, ApJ, 755, 6 Zhu, Z., Andrews, S. M., & Isella, A. 2018, MNRAS, 479, 1850 A95, page 13 of 19 A&A 667, A95 (2022) Appendix A: The RT diffusion solver P RO D I M O 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 divF = Γ , (A.1) R where F = Fν dν [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 Z   4π κνabs Bν (T ) − Jν dν = Γ , (A.2) where κνabs is the dust absorption coefficient at frequency ν, Bν (T ) the Planck function, and Jν 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 P RO D I M O 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 Jν (r, z) which are used to update the dust temperatures T dust (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π grad J(r, z) 3κR (r, z) (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 Z . Z dBν (T ) 1 1 dBν (T ) = dν dν (A.4) κR (r, z) dT κext (r, z) dT Z ν Z . κP (r, z) = κνabs (r, z)Bν (T ) dν Bν (T ) dν , (A.5) where κνext (r, z) is the extinction coefficient and T = T dust (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 frequencydependent disk opacity structure and the current T dust (r, z). Next, we R 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} . (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 T i, 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 A95, page 14 of 19 Zi+1/2,j+1/2 Zi,j+1/2 Zi-1/2,j+1/2 Zi+1/2,j-1/2 Zi,j Zi,j-1/2 Zi-1/2,j-1/2 ri-1 ri-1/2 ri ri+1/2 ri+1 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. mean intensities on the boundary points, we integrate Eq. (A.2) assuming that Jν is close to a Planckian, hence J(r, z) = B(T ) − Γ(r, z) . 4π κP (r, z) (A.7) where B(T ) = σ T dust (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 Ji, j − Ji−1, j Ji, j − Ji+1, j + Aver (A.8) 1/2, j Di+ 1/2, j i+ ∆ri−1/2 ∆ri+1/2 Ji, j − Ji, j−1 Ji, j − Ji, j+1 + Ahor + Ahor = Vi, j Γi, j , i Di, j−1/2 i Di, j+1/2 ∆zi, j−1/2 ∆zi, j+1/2 Aver i−1/2, j Di−1/2, j 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 Ahor i 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 P RO D I M O grid points, which are aligned on radial rays on which z/r is constant √ (A.9) ri−1/2 = ri ri−1 √ 1 ri+ /2 = ri ri+1 (A.10) ∆ri−1/2 = ri − ri−1 (A.11) ∆ri+1/2 = ri+1 − ri (A.12) 1 zi, j−1/2 = (zi, j + zi, j−1 ) (A.13) 2 1 (A.14) zi, j+1/2 = (zi, j + zi, j+1 ) 2 ri−1/2 zi−1/2, j−1/2 = zi, j−1/2 (A.15) ri ri−1/2 zi−1/2, j+1/2 = zi, j+1/2 (A.16) ri ri+1/2 zi+1/2, j−1/2 = zi, j−1/2 (A.17) ri ri+1/2 zi+1/2, j+1/2 = zi, j+1/2 (A.18) ri Aver (A.19) i−1/2, j = 2π ri−1/2 (zi−1/2, j+1/2 − zi−1/2, j−1/2 ) Aver i+1/2, j = 2π ri+1/2 (zi+1/2, j+1/2 − zi+1/2, j−1/2 ) Ahor i = Vi, j = 2 2 π(ri+ 1/2 − ri− 1/2 ) Ahor i (zi, j+1/2 − zi, j−1/2 ) (A.20) (A.21) (A.22) N. Oberg et al.: Circumplanetary disk ices. I. The radiative diffusion coefficients are defined as 4π Di, j = 3κR (ri , zi, j ) p Di−1/2, j = Di, j Di−1, j p Di+1/2, j = Di, j Di+1, j p Di, j−1/2 = Di, j Di, j−1 p Di, j+1/2 = Di, j Di, j+1 (A.23) (A.24) (A.25) (A.26) (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.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 Γi, 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) Γ(r, z) 4π κP (r, z) π 1/4 T dust (r, z) = B(T ) σ B(T ) = J(r, z) + Fig. A.2: Benchmark test for a viscous circumstellar disk with a disk mass of 0.01 M⊙ and a mass accertion rate of Ṁacc = 10−8 M⊙ /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 P RO D I M O. Figure A.2 reveals a number of interesting features in the disk temperature structure: – The dust temperature at the inner rim is not much affected by viscous heating. (A.29) (A.30) and multiply the frequency-dependent mean intensities Jν (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 Jν (r, z). After having modified T dust (r, z) and Jν (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 Ṁacc = 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. Fig. A.3: Midplane temperature T dust (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 P RO D I M O. A95, page 15 of 19 A&A 667, A95 (2022) – 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 T dust 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 P RO D I M O 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. 100 Appendix B: Adsorption energies Table B.1: Adsorption energies of the most prevalent molecular ices found in our model CPDs. Ice H2 O NH3 NH2 C2 H6 C2 H4 C2 H2 CO2 CH3 OH OH Eads [K] 4800 5534 3956 2300 3487 2587 2990 5534 2850 ref. a b b c b d e a b (a) (Brown & Bolina 2006) (b) (Garrod & Herbst 2006) (c) (Öberg et al. 2009) (d) (Collings et al. 2004) (e) (McElroy et al. 2013) The adsorption energies of our most common ices and their respective references are listed in Table B.1. T-change 10-1 Appendix C: Surface density slope 10-2 10-3 10-4 10-5 0 50 100 RT iterations 150 Our reference CPDs have a surface density powerlaw exponent ǫ = 1. The steady-state solution for a constant- Ṁ 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 H2 O formation. 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. 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−3.3 . A95, page 16 of 19 N. Oberg et al.: Circumplanetary disk ices. I. 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 iceline (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 T back 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. Fig. D.1: Midplane ice mass fraction for a variety of background temperatures T back 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 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. 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 T dust = 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. A95, page 17 of 19 A&A 667, A95 (2022) 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 midplane, 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 Parameter Value Unit Hydrogen density 104 cm−3 Gas temperature 10.0 K Dust temperature 10.0 K Optical extinction 10.0 Mean grain radius 0.1 ➭m Cloud Lifetime 1.7 × 105 yr Note: Parameter values of the molecular cloud model and integration time are chosen according to the method of Helling et al. (2014) recommended for TMC-1 (Taurus Molecular Cloud) by McElroy et al. (2013). Initial atomic abundances intended to represent typical diffuse interstellar medium conditions are also adopted from (McElroy et al. 2013). A95, page 18 of 19 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). Error bars on model values represent a factor 10x uncertainty for illustrative purposes. TMC-1 abundances are taken from (Suutarinen et al. 2011) (OH,CH), (Sakai et al. 2010) (C2 H), (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−3.3 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. N. Oberg et al.: Circumplanetary disk ices. I. (a) High-Mass High-Viscosity (7-10) reset (b) High-Mass High-Viscosity (7-10) inherit (c) High-Mass Low-Viscosity (7-11) reset (d) High-Mass Low-Viscosity (7-11) inherit (e) Low-Mass High-Viscosity (8-11) reset (f) Low-Mass High-Viscosity (8-11) inherit (g) Low-Mass Low-Viscosity (8-12) reset (h) Low-Mass Low-Viscosity (8-12) inherit 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). A95, page 19 of 19