Issue 
A&A
Volume 686, June 2024



Article Number  A103  
Number of page(s)  10  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202449237  
Published online  03 June 2024 
Thermal structure of circumbinary discs: Circumbinary planets should be icy, not rocky
^{1}
Laboratoire d’astrophysique de Bordeaux, Univ. Bordeaux, CNRS,
B18N, allée Geoffroy SaintHilaire,
33615
Pessac,
France
email: arnaud.pierens@ubordeaux.fr
^{2}
Astronomy Unit, Queen Mary University of London,
Mile End Road,
London
E1 4NS,
UK
Received:
15
January
2024
Accepted:
7
March
2024
The process of forming a circumbinary planet is thought to be intimately related to the structure of the nascent circumbinary disc. It has been shown that the structure of a circumbinary disc depends strongly on threedimensional effects and on detailed modelling of the thermodynamics. Here, we employ threedimensional hydrodynamical simulations, combined with a proper treatment of the thermal physics using the RADMC3D radiation transport code, to examine the location of the snow line in circumbinary discs. The models have application to the circumbinary planets that have been discovered in recent years by the Kepler and TESS transit surveys. We find that the snow line is located in a narrow region of the circumbinary disc, close to the inner cavity that is carved out by the central binary, at typical orbital distances of ~ 1.5–2 au for the system parameters considered. In this region, previous work has shown that both grain growth and pebble accretion are likely to be inefficient because of the presence of hydrodynamical turbulence. Hence, in situ planet formation interior to the snow line is unlikely to occur and circumbinary planets should preferentially be icy, not rocky.
Key words: accretion, accretion disks / hydrodynamics / methods: numerical / planets and satellites: formation / protoplanetary disks / planetdisk interactions
© The Authors 2024
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.
1 Introduction
The discovery of circumbinary planets (CBPs) such as Kepler16b (e.g. Doyle et al. 2011) has sparked new interest in how planets can form in highly perturbed environments such as circumbinary discs. Due to the constraints resulting from gravitational perturbations induced by the central binary, circumbinary discs are ideal laboratories for testing planet formation theories. The different scenarios for planet formation in circumbinary discs are generally based on the assumption that CBPs were either formed in situ, at the location where they are currently observed, or at large distances from the binary. In the latter case, an episode of inward migration must be invoked to explain the present orbits of CBPs, which are generally found to orbit close to the zone of dynamical stability defined by Holman & Wiegert (1999).
Forming CBPs in situ has been proven to be difficult in the context of standard planetesimal accretion (e.g. Paardekooper et al. 2012). It has been shown that large planetesimal eccentricities can result from the tidal interaction with the binary or with the eccentric circumbinary disc (Marzari et al. 2008; Kley & Nelson 2010), resulting in collision velocities that are too high to enable planetesimal growth. Growing a massive planetesimal seed in situ through pebble accretion also appears to be very challenging. The inner disc is prone to a parametric instability that generates hydrodynamical turbulence (Papaloizou 2005; Barker & Ogilvie 2014), the main impact of which is to stir up the pebble layer, resulting in inefficient pebble accretion (Pierens et al. 2020, 2021).
Taken together, the above constraints seem to favour the migration scenario, although this scenario has its own drawbacks. In particular, the stopping locations of CBPs migrating in 2D disc models can be a poor match to the locations of some of the known CBPs (Pierens & Nelson 2013; Penzlin et al. 2021), especially for higheccentricity binaries like Kepler34 (Welsh et al. 2012). This arises because a wide and eccentric tidally induced inner cavity is formed in the circumbinary disc, and the planet itself tends to acquire significant eccentricity, causing its inward migration to stall too far from the central binary. In contrast, recent studies have shown that when the disc is circularized by a planet opening a partial gap (Mutter et al. 2017; Penzlin et al. 2021), or because of dust feedback onto the gas when the inner regions of the disc become enriched by inwarddrifting dust (Coleman et al. 2022), a better fit to the observed orbital parameters can be obtained. Compared to 2D discs, Pierens & Nelson (2023) also find a systemic trend for the size and eccentricity of the inner cavity to be smaller in 3D models that include the effects of stellar heating and radiative cooling on the thermal disc structure. These latter results underline the importance of considering 3D effects combined with a proper treatment of the disc thermodynamics to obtain realistic structures for circumbinary discs.
In this context, we present the results of 3D hydrodynamical simulations of circumbinary discs that include irradiation from the central stars and cooling to estimate the location of the snow line in circumbinary discs, with application to several of the systems known to host CBPs. There have been several previous studies that have focused on the position of the snow line. Using a static model of a circumbinary disc that is heated by stellar irradiation and viscous dissipation, Clanton (2013) found that ice lines tend to lie within the zone of instability, and suggested that rocky planets should not form in these systems. The effect of tidal dissipation was examined in globally viscous circumbinary discs by Shadmehri & Khajenabi (2015) and Vartanyan et al. (2016), who found that the combination of viscous dissipation and the damping of density waves driven by the binary tends to push the ice line out beyond 3–5 au. With respect to these previous studies, our work selfconsistently computes the density and temperature structure of the disc using 3D hydrodynamical simulations coupled with Monte Carlo radiative transfer calculations. We also examine the influence of localised viscous heating that arises because of the dissipation of the binaryinduced hydrodynamical turbulence described above. We work under the assumption that the low ionisation fraction of the circumbinary discs prevents the development of magnetised turbulence within them (e.g. Gammie 1996), such that the only source of heating far from the tidally induced inner cavity is stellar irradiation. We investigate the dependence of the results on the parameters of the central binary, and apply the results to the putative circumbinary discs of several systems known to host CBPs.
This paper is structured as follows. We describe the numerical model in Sect. 2 and we present the setup employed for the simulations in Sect. 3. In Sect. 4, we present the results of our 3D radiationhydrodynamical simulations. In Sect. 5, we discuss our results and draw conclusions regarding the location of the snow line in circumbinary discs and the implications that this has for the likely composition of CBPs.
2 Numerical model
2.1 Threedimensional hydrodynamical simulations
2.1.1 Equations of motion
We solved the hydrodynamical equations for the conservation of mass, momentum, and internal energy in spherical coordinates (r, θ, ϕ) (radial, polar, azimuthal), with the origin of the frame located at the centre of mass of the binary. In particular, the energy equation that is incorporated in the code is given by $$\frac{\partial e}{\partial t}+\nabla \cdot (ev)=(\gamma 1)e\nabla \cdot v+{Q}_{\text{cool}}+{Q}_{\text{bulk}}^{+},$$(1)
where v is the velocity, and e = ρc_{v}T is the internal energy per unit volume, with ρ the density, c_{v} the specific heat capacity at constant volume, and T the temperature. In the previous equation, γ is the adiabatic index, which was set to γ = 1.4.
We considered a cooling scheme where the temperature is relaxed towards a reference temperature, T_{irr}, on a cooling timescale, t_{cool}. In Eq. (1), the cooling rate, Q_{cool}, is therefore given by $${Q}_{\text{cool}}=\rho {c}_{\upsilon}\frac{T{T}_{\text{irr}}}{{t}_{\text{cool}}}.$$(2)
Here, we assumed that T_{irr} corresponds to the temperature set by stellar irradiation, augmented by viscous heating if this is included in the model. During the simulations, this temperature was regularly updated using the RADMC3D Monte Carlo radiative transfer code (Dullemond et al. 2012), applying the procedure described below in Sect. 2.2. The cooling timescale, t_{cool}, was calculated at every time step by considering the timescale for radiative loss of energy from a Gaussian sphere with a characteristic length scale corresponding to the gas scale height, H (Lyra et al. 2016). It is given by $${t}_{\text{cool}}=\frac{\rho {c}_{\upsilon}H{\tau}_{\text{eff}}}{3{\sigma}_{\text{SB}}{T}^{3}},$$(3)
where σ_{SB} is the StefanBoltzmann constant and τ_{eff} is the effective optical depth given by $${\tau}_{\text{eff}}=\frac{3}{8}\tau +\frac{\sqrt{3}}{4}+\frac{1}{4\tau}.$$(4)
Following Lyra et al. (2016), the optical depth, τ, was calculated as 1/τ = 1/τ_{upper} + 1/τ_{lower}, with $${\tau}_{\text{upper}}={{{\displaystyle \int}}^{\text{}}}_{z}^{{z}_{\mathrm{max}}}\rho \left({z}^{\prime}\right)\kappa \left({z}^{\prime}\right)\text{d}{z}^{\prime}$$(5)
and $${\tau}_{\text{lower}}={{{\displaystyle \int}}^{\text{}}}_{{z}_{\text{min}}}^{z}\rho \left({z}^{\prime}\right)\kappa \left({z}^{\prime}\right)\text{d}{z}^{\prime},$$(6)
where the opacity, κ, was computed using the Rosseland mean opacity of Zhu et al. (2009).
Finally, the heating source term, ${Q}_{\text{bulk}}^{+}$, is provided by artificial viscous heating resulting from shock damping of the density waves excited by the binary.
2.1.2 Numerical methods
The simulations presented in this paper were performed using the multifluid version of FARGO3D (BenítezLlambay & Masset 2016). Computational units were chosen such that the total mass of the binary is M_{★} = 1, the gravitational constant G = 1, and the radius, R = 1 in the computational domain, corresponds to the binary semimajor axis, a_{bin}. When presenting the simulation results, unless otherwise stated we use the binary orbital period, ${T}_{\text{bin}}=2\pi \sqrt{{a}_{\text{bin}}^{3}/G{M}_{\star}}$, as the unit of time.
The computational domain in the radial direction extends from R_{in} = 1.13 a_{bin} to R_{out} = 50 a_{bin} and we employed 684 logarithmically spaced grid cells. In the azimuthal direction the simulation domain extends from 0 to 2π, with 384 uniformly spaced grid cells. In the meridional direction, the simulation domain covers 3.5 disc pressure scale heights above and below the disc midplane, and we adopted 48 uniformly spaced grid cells. We note that this resolution allows us to compute the quasisteady structure of 3D circumbinary disc models for a range of binary parameters, given the computational resources at our disposal. It is not, however, sufficient to directly capture the hydrodynamical turbulence that is excited through the parametric instability in the inner regions of inviscid circumbinar discs (see Pierens et al. 2021). It is for this reason that we included a viscous heating term in the radiative transfer calculations described below, since we want to estimate how the dissipation of the turbulence affects the position of the snow line.
2.1.3 Boundary conditions
We adopted an outflow boundary condition at the inner edge to allow mass to accrete onto the binary. Here, the velocity was set to 0 if directed into the computational domain, such that gas can leave the domain but not flow into it. The values for the gas density and internal energy in the ghost zones are the same as in the first active zone. In the context of 2D models, it has been shown that employing an open boundary together with a location of the inner disc edge, R_{in} ≈ a_{bin}, leads to a quasistationary disc structure after O(10^{4}) binary orbits. Compared to other boundary conditions, adopting an open boundary also appears to be less sensitive to numerical issues (Thun et al. 2017).
At the outer radial boundary, we employed a wavekilling zone for R > 0.88 R_{out} to avoid wave reflection at the disc’s outer edge (de ValBorro et al. 2006). The impact of this wavekilling zone on the disc shape and eccentricity is expected to be small since the disc structure near the wavekilling zone remains close to the initial one.
At the meridional boundaries of the 3D domain, an outflow boundary condition was also employed, with the exception that for the gas density we followed Bae et al. (2016) and maintained vertical stratification by solving the following condition for hydrostatic equilibrium: $$\frac{1}{\rho}\frac{\partial}{\partial \theta}\left({c}_{\text{s}}^{2}\rho \right)=\frac{{\upsilon}_{\varphi}^{2}}{\mathrm{tan}\theta},$$(7)
where c_{s} is the sound speed and v_{φ} the azimuthal velocity.
2.2 Radiative transfer calculations
2.2.1 Radiation sources
The radiation transfer calculations were performed using RADMC3D (Dullemond et al. 2012), in which radiation from multiple stars can be incorporated, as well as an extra internal source of energy such as viscous heating. We assumed that both stars are located at the centre of the domain and did not take into account the timedependent effects arising from their orbital motion. Regarding the emission from the stars, we first considered their blackbody radiation fields using N = 10^{8} photons. In some models, we also took into account localised viscous heating driven by turbulent viscosity, using the modified random walk method (Min et al. 2009; Robitaille 2010). In this case, we employed a reduced number of photons N = 10^{7} due to the higher computational cost required to run these models. The corresponding heat source that was employed in RADMC3D is given by $${\text{\Gamma}}_{v}=\frac{9}{4}\frac{G{M}_{\star}}{{R}^{3}}\rho v,$$(8)
where the kinematic viscosity is given by v = αH^{2}Ω. in the context of the α model (Shakura & Sunyaev 1973). Previous work has shown that the most plausible source of turbulence in the inner regions of circumbinary discs is a parametric instability driven by the resonant interaction between inertialgravity waves and an eccentric mode in the disc (Papaloizou 2005; Barker & Ogilvie 2014). Because it is directly related to the disc eccentricity, turbulence generated in this way is more vigorous close to the inner cavity where the disc eccentricity is the highest, and operates much less efficiently when moving away from the inner cavity. A direct consequence is that the α parameter is not constant, with a peak value of α ≈ 8 × 10^{−3} at the cavity edge and α ≈ 2 × 10^{−3} at a distance of R ~ 6 a_{bin}.·In this work, we adopted the following functional dependence of α on radius: $$\alpha (R)=0.01\frac{\mathrm{exp}\left({R}_{\text{gap}}R\right)}{1+\mathrm{exp}\left(\frac{R{R}_{\text{gap}}}{0.05{R}_{\text{gap}}}\right)},$$(9)
where R_{gap} is the estimated gap size (see below). This form was chosen to reproduce the α profile found in Pierens et al. (2021) (see their Fig. 4), and is plotted in Fig. 1.
Fig. 1 α viscous parameter as a function of distance from the binary. The dependence of α on the radius was chosen so as to reproduce the α profile found in Pierens et al. (2021; see their Fig. 4), where hydrodynamical turbulence operates. 
2.2.2 Grain properties and dust density profiles
The chosen grain size distribution is composed of 10 size bins, logarithmically distributed between a minimum grain size of α_{min} = 0.01 µm and a maximum grain size of α_{max} = 1 mm. The fraction of dust mass in each size bin was calculated assuming a dusttogas ratio є = Σ_{d}/Σ_{g} = 0.01, where Σ_{g} and Σ_{d} are the gas and dust surface densities, respectively. The grain sizes are described by a MRN distribution, for which the number density of grains in the range [α, α + dα] is given by n(a)da∝ a^{−3.5} (Mathis et al. 1977). We note that the RADMC3D calculations were performed using these ten grain sizes as ten different dust species, each having its own opacity. The opacity of each grain size was calculated using optool (Dominik et al. 2021) and by considering the DSHARP opacities (Birnstiel et al. 2018). These assume that the dust is composed of water ice, astronomical silicates, troilite, and refractory organic material, with volume fractions being 36%, 17%, 3%, and 44%, respectively. This results in a bulk density of the dust mixture of ρ_{p} = 1.675 g cm^{−3}. The dust volume density was defined by imposing a Gaussian vertical profile: $${\rho}_{\text{d}}(R,z,a)=\frac{{\text{\Sigma d}}_{}(R,a)}{\sqrt{2\pi}{H}_{\text{d}}(R,a)}\mathrm{exp}\left(\frac{{z}^{2}}{2{H}_{\text{d}}{(R,a)}^{2}}\right).$$(10)
To account for dust settling, the particle scale height, H_{d}(R, a), varies with grain size, a, and is given by $${H}_{\text{d}}(R,a)={H}_{\text{g}}(R)\sqrt{\frac{\alpha}{\text{St}+\alpha}},$$(11)
where St is the Stokes number calculated in the Epstein regime: $$\text{St}=\frac{{\rho}_{\text{p}}a}{\rho {c}_{\text{s}}}\text{\Omega}.$$(12)
For low levels of turbulence, particles tend to be strongly confined in the disc midplane, resulting in a very optically thick region. To avoid this, we employed a residual viscosity with α_{res} = 3 × 10^{−5} that limits the scale height of mmsized grains to values of H_{d} = 0.1−0.2H_{g}, consistent with the values reported in Pierens et al. (2021; see their Fig. 6).
3 Model setup
3.1 Initial conditions
The initial gas surface density profile is given by $${\text{\Sigma g}}_{}={f}_{\text{gap}}{\text{\Sigma}}_{\text{g},0}{\left(\frac{R}{{R}_{0}}\right)}^{1/2},$$(13)
where Σ_{0} is the surface density at R_{0} = a_{bin}, which was defined such that the mass of the gas disc is 5% of the total binary mass within 40 au. ƒ_{gap} is a gap function used to initiate the disc with an inner cavity (assumed to be created by the binary), and is given by $${f}_{\text{gap}}={\left(1+\text{exp\hspace{0.17em}}\left[\frac{R{R}_{\text{gap}}}{0.1{R}_{\text{gap}}}\right]\right)}^{1},$$(14)
where R_{gap} is the estimated gap size, which was set to R_{gap} = 2.5 a_{bin}. The initial azimuthal velocity was set to the Keplerian velocity, whereas both the initial radial and meridional velocities were set to zero.
The initial temperature in the disc corresponds to the stellar irradiationdominated temperature, $${T}_{\text{irr}}(R)={\left(\frac{{L}_{\star}}{4\pi {R}^{2}{\sigma}_{\text{SB}}}\right)}^{1/4},$$(15)
where L_{⋆} is the total stellar luminosity. Using Eqs. (13) and (15), one can deduce the 3D gas density, which, combined with the dust volume density given by Eq. (10), can be subsequently employed in RADMC3D to compute the dust temperature. The new gas temperature was then simply obtained by considering the gas and dust temperatures to be equal, which is a reasonable assumption for µmsized grains that tend to be strongly coupled to the gas. After this radiative transfer calculation, the new gas density was then calculated by FARGO3D and we iterated this procedure every binary orbit until we converged on a quasiequilibrium solution, which is typically obtained after ~5 iterations. Over longer evolution timescales, we increase the interval between two radiative transfer calculations to ten and then 100 binary orbits until a quasistationary state is reached.
3.2 Binary parameters
It can be reasonably expected that over the ~10^{7} yr of circumbinary disc evolution, the characteristics of the central binary (temperature, luminosity, etc.) change significantly as well. In the context of CBP formation, this is an important issue, as it may influence the amount of stellar heating received by the disc, resulting in a change in the structure of the inner cavity. This is illustrated in Fig. 2, where we show evolutionary tracks on a HR diagram for the circumbinary systems that we will focus on in this paper. The stellar luminosity continuously decreases as it follows the Hayashi track, until it reaches a value close to the luminosity at the zeroage main sequence (ZAMS). The overplotted stellar isochrone at 1 Myr confirms that at the formation epoch of a CBP, the binary star’s combined luminosity was probably much higher than it is at the present day. For the 11 known CBP systems, the values of the luminosities after 1 Myr of evolution are listed in Table 1 and the orbital parameters are listed in Table 2. Representing the binary as a point source located at the centreofmass, the last column gives the corresponding irradiation flux, F_{in}, at a distance of a_{bin}. In terms of the values of F_{in}, the circumbinary systems that are presented in Fig. 2 span the range of values that are obtained across the whole set of known circumbinary systems. Hence, we focus on these six particular systems in the rest of this paper.
We note that over an evolution timescale of 1 Myr, the binary separation is also expected to be impacted by its interaction with the circumbinary disc (e.g. Turpin & Nelson 2024). For simplicity, we here assumed a fixed value for a_{bin} but it should be kept in mind that compared to its present value, the binary semimajor axis was also probably slightly larger when CBPs formed.
Fig. 2 Evolutionary tracks on the HR diagram for different circumbinary systems. The upper dashed line shows the stellar isochrone at 1 Myr, whereas the lower dashed line corresponds to the ZAMS. Isochrones were computed using MIST evolutionary models (Choi et al. 2016). 
4 Results
For each system considered, Fig. 3 shows the azimuthally averaged surface density, midplane temperature, and disc aspect ratio profiles at time=10^{3} T_{bin}. This corresponds to the evolution timescale over which the circumbinary disc structure reached a quasisteady state in 3D models (Pierens & Nelson 2023). We see that inviscid discs and those with viscous heating have similar profiles in all quantities (especially the surface density), except for Kepler34 where the size of the cavity in the model with viscous heating is slightly larger. As is discussed below, this largely arises because the viscous dissipation is localised to the regions just outside of the inner cavity, but also because the localised viscous heating rate is similar to the shock dissipation and irradiation heating rates, so that viscosity does not play a dominant role in determining the density and temperature structure of the inner regions of circumbinary discs.
The filled circles and triangles in Fig. 3, representing the locations of the snow lines at the disc midplanes, show that these are located close to the cavity edge at distances of between ~0.7– 1.8 au from the central binary, depending on which system is being considered and whether or not viscous heating has been included. As was expected, the effect of including viscous heating is to move the snow line outwards, but in general this is a small effect. In the models with viscous heating, the snow line sits out beyond the location where the viscosity is at its maximum value. Hence, it is not just the local viscous heating that plays a role but also the radiative diffusion radially outwards of the viscously generated heat that is important^{1}.
Looking at the H/R profiles in Fig. 3, we observe slightly larger disc aspect ratios in the innermost regions of the discs with viscous heating included. In all models, the disc aspect ratio exhibits significant bumps near the inner cavity where all the heating processes operate most effectively. This leads to the formation of a selfshadowed region in the outer disc. Compared to an inviscid disc, a viscous disc tends to have a higher aspect ratio in the vicinity of the density peak because viscous heating scales with gas density and the viscous dissipation is concentrated in this region. The resulting enhanced selfshadowing effect leads to a slightly smaller temperature and aspect ratio in the outer disc in comparison to an inviscid disc.
In Fig. 4, we show the azimuthally averaged temperature distributions in the (R–Z) plane for each system considered. Again, only small differences between models with or without viscous heating are observed. Compared to a viscous disc, the regions beyond the inner cavity and density maximum of an inviscid disc tend to have a slightly higher temperature at the surface, presumably because more stellar radiation is intercepted in the inner disc regions in the viscous models because of the higher aspect ratios there.
We also observe that in the vicinity of the density peak, and out just beyond it, the temperature of an inviscid disc can reach T ~ 500–600 K in the midplane. This is a consequence of the shock dissipation of density waves excited by the binary, which provides an important source of additional heating, as can be seen in the contours in Fig. 4, where the wavelike structures illustrate the heating contribution from shock heating. The effect is particularly pronounced for the Kepler34 run, and Fig. 5, which shows the different contributions to the change in specific internal energy, shows that the heating efficiency due to shocks can be of the same order as the viscous and stellar heating combined near the tidally truncated cavity.
The importance of the dissipation of spiral shocks in setting the disc thermodynamical structure of circumbinary discs has also been emphasised by Vartanyan et al. (2016). In their models, which involve globally viscous discs in which the torque from the central binary prevents inward mass accretion at the inner edge of the cavity, and in which heating is provided by viscosity, shock dissipation, and stellar irradiation, the snow line is typically found to lie beyond 36 au, depending on the model and the evolutionary stage under consideration. These values are clearly much larger than those obtained by the models presented here, and arise because of the fundamentally different underlying disc models that they considered.
The locations of the snow lines for the different circumbinary systems considered here are shown in Fig. 6 by the green circles overplotted on the surface density images. As was described above, the snow line is located slightly further away from the binary in viscous models compared to inviscid ones, typically at radial distances ≲ 1.5 au. More importantly, we find that in both inviscid and viscous discs, the snow line tends to be located close to the density maximum and inner cavity. Hence, there is only a narrow range of radial locations in the discs where the midplane temperature is ≳ 170 K, such that it is difficult to envisage rocky planets that contain no icy material being able to form. This statement is reinforced by the fact that high levels of hydrodynamical turbulence are expected to be present near the inner cavity, such that pebble accretion onto seed protoplanets would proceed in the inefficient, 3D regime in which the pebble layer has a scale height that exceeds the Hill sphere radius of the accreting body. In situ grain growth is also expected to be difficult in this region due to the high collision velocities that are also expected to arise because of the turbulence. Indeed, Pierens et al. (2021) found that at distances ≲8–9 a_{bin}, which, for example, would correspond to ~2 au for Kepler16 and ~0.7 au for Kepler47, collision velocities are higher than the fragmentation velocities of silicate aggregates. Taken together, these factors suggest that planets forming in circumbinary discs are likely to form exterior to the snow line, and hence will be composed of a significant water ice fraction.
Binary stellar parameters for known CBP hosts taken from Georgakarakos et al. (2021).
Binary orbital parameters for all known CBP hosts.
Fig. 3 Surface density, midplane temperature, and aspect ratio profiles of the different circumbinary systems that we considered. Solid lines correspond to models that include viscous heating, whereas dashed lines are for models without viscous heating. All models include irradiation and shock heating. In the surface density profiles, the markers indicate the location of the ice line, which is also represented as a horizontal dashed line in the temperature plots. The dashed black lines in the surface density plots show the values of α in the viscously heated models – see the righthand vertical axis. 
Fig. 4 Temperature distributions of each of the irradiated and viscously heated models that we considered. 
Fig. 5 Contribution of the different heating terms to the rate of specific internal change for Kepler34. 
5 Discussions and conclusions
In this paper, we have presented the results of 3D hydrodynamical simulations that are tightly coupled to Monte Carlo radiative transfer calculations to determine the location of the snow line in circumbinary discs. We chose binary parameters that match those of binary systems that are known to host CBPs (e.g. Kepler16, 34).
The simulations include a cooling scheme where the disc temperature is relaxed towards an irradiation temperature on a finite cooling timescale that was calculated using realistic opacities. The irradiation temperature is determined by performing Monte Carlo radiative transfer calculations using RADMC3D, in which localised viscous heating due to hydrodynamical turbulence was incorporated as an additional heating source for a subset of the models. Our method involves evolving the hydrodynamical and thermodynamical states of the disc in lockstep by performing radiative transfer calculations at different stages during the hydrodynamical simulations. Hence, the evolution of the density and temperature structures of the disc feedback on each other as a quasisteady state is established, and the effects of stellar irradiation, shock heating due to the dissipation of binaryinduced spiral density waves, and viscous dissipation (when switched on) are included in the models.
We find that the snow line lies close to the density maximum that sits just outside of the tidally truncated cavity in all the circumbinary discs. In discs that include viscous heating, the ice line lies further out, but only by a small amount. The similar disc structures obtained in viscous and inviscid models are a consequence of shock heating, stellar irradiation, and viscous dissipation all providing a similar amount of heating in the disc, and because these heat sources are quite localised in the inner regions of the disc close to the cavity.
In the outer disc, heating at shocks and through viscous dissipation are obviously much smaller. The latter depends on the level of turbulence operating in the disc and how this is spatially distributed. The origin of the turbulence is a parametric instability that only operates close to the inner cavity where the disc becomes visibly eccentric, and hence viscous heating is only important in a narrow region close to the surface density maximum. This results in typical snow line locations of ≲1.5 au from the binary. At these distances, pebble accretion is believed to be inefficient due to the vertical stirring of dust and pebbles induced by the turbulence. Similarly, the high relative collision velocities there reduce the likelihood of in situ grain growth to sizes sufficient to trigger the streaming instability.
These considerations indicate that it will be difficult to form rocky planets that are devoid of significant fractions of water ice in circumbinary discs, and we suggest that CBPs should preferentially be icy and not rocky. A similar conclusion was reached by Clanton (2013) using a simpler disc model than the one presented in this paper. This statement should be particularly true for planets evolving on orbits that are coplanar to the binary orbit. This is the case for most of the CBPs detected so far but we note that this may be partly due to an observational bias, as misaligned CBPs are more difficult to detect due to their nonperiodic transits (Martin & Triaud 2014). Nevertheless, the observation of misaligned circumbinary discs with ALMA (e.g. Czekala et al. 2019), and the detection of Kepler413b (Kostov et al. 2014) and Kepler453b (Welsh et al. 2015) whose orbital planes are slightly misaligned with respect to the binary orbital plane, indicate that misaligned CBPs do exist. It has even been proposed that CBPs around polar orbits may form with comparable efficiency to coplanar planets (Childs & Martin 2021). Such systems may form in circumbinary discs that are polar aligned to the binary orbit. For an eccentric binary, a circumbinary disc with a sufficient initial inclination with respect to the binary orbit may achieve a polar configuration due to viscous dissipation (Martin & Lubow 2017). Compared to the coplanar case, however, the tidal torque exerted on a disc in a polar state is much weaker, resulting in an inner disc edge that lies closer to the binary (Martin & Lubow 2018; Franchini et al. 2019), plausibly inside the ice line. This suggests that if rocky planets can form in circumbinary discs, these planets would preferentially be found on orbits that are highly inclined with respect to the binary orbit.
The prediction that CBPs should be icy, not rocky receives tentative support from current observations of CBPs. All CBPs discovered so far have radii ≥3 R_{⊕} and appear to have significant hydrogenhelium envelopes, so in this regard they are more similar to the ice and gas giant planets of the Solar System than the rocky terrestrial planets. Furthermore, the known population of CBPs differs significantly from the exoplanet population that orbits around single stars, which is dominated by smaller superEarths and miniNeptunes.
Attempts have been made to search for additional transiting CBPs orbiting eclipsing binaries in the Kepler data, using customised algorithms that account for the large transit timing variations inherent in such systems (Martin & Fabrycky 2021, 2022), but these searches have failed to find any new planets, even though they should have been sensitive to finding planets smaller than 3 R_{⊕} if they were present in the data (Martin & Fabrycky 2022).
While it is currently not possible to determine the compositions of the cores of the known CBPs, it is tempting to infer that the lack of smaller CBPs arises because their formation involves the incorporation of a substantial fraction of icy material that augments the masses of their cores, leading to the accretion of significant hydrogen envelopes. Future observations by TESS, and discoveries by the PLATO mission, will enhance the population of known CBPs and will provide direct measurements of their radii. Radial velocity surveys such as the Binaries Escorted By Orbital Planets (BEBOP; Martin et al. 2019) survey are also detecting and discovering CBPs using the radial velocity method (Triaud et al. 2022; Standing et al. 2023). Hence, it is hoped that highprecision radial velocity measurements of future planets discovered by transit surveys will provide measurements of their density, and by extension constraints on their bulk compositions. Eventually, this will allow a CBP massradius diagram to be constructed, providing insight into the compositional diversity of the CBP population.
Fig. 6 Surface density distribution for each of the simulations. The location of the ice line is indicated by the green circle, while the locations of the stars are represented as white circles. 
Acknowledgments
Computer time for this study was provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Université de Bordeaux and by HPC resources of Cines under the allocation A0110406957 made by GENCI (Grand Equipement National de Calcul Intensif). RPN acknowledges support from the Leverhulme Trust through grant number RPG2018418 and from STFC through grants ST/X000931/1 and ST/T000341/1.
References
 Bae, J., Nelson, R. P., & Hartmann, L. 2016, ApJ, 833, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Barker, A. J., & Ogilvie, G. I. 2014, MNRAS, 445, 2637 [NASA ADS] [CrossRef] [Google Scholar]
 BenítezLlambay, P., & Masset, F. S. 2016, ApJS, 223, 11 [Google Scholar]
 Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45 [CrossRef] [Google Scholar]
 Childs, A. C., & Martin, R. G. 2021, ApJ, 920, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102 [Google Scholar]
 Clanton, C. 2013, ApJ, 768, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Czekala I., Chiang E., Andrews S. M., et al. 2019, ApJ, 883, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Coleman, G. A. L., Nelson, R. P., & Triaud A. H. M. J. 2022, MNRAS, 513, 2563 [NASA ADS] [CrossRef] [Google Scholar]
 de ValBorro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 [Google Scholar]
 Dominik, C., Min, M., & Tazaki, R. 2021, Astrophysics Source Code Library [record ascl:2104.010] [Google Scholar]
 Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Science, 333, 1602 [NASA ADS] [CrossRef] [Google Scholar]
 Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, Astrophysics Source Code Library [record ascl:1202.015] [Google Scholar]
 Franchini, A., Lubow, S. H., & Martin, R. G. 2019, ApJ, 880, L18 [Google Scholar]
 Gammie, C. F. 1996, ApJ, 457, 355 [Google Scholar]
 Georgakarakos, N., Eggl, S., & DobbsDixon, I. 2021, BAAS [Google Scholar]
 Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621 [Google Scholar]
 Kley, W., & Nelson, R. P. 2010, Planets Binary Star Syst., 366, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Kostov, V. B., McCullough, P. R., Hinse, T. C., etal. 2013, ApJ, 770, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Kostov, V. B., McCullough, P. R., Carter, J. A., et al. 2014, ApJ, 784, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Kostov, V. B., Orosz, J. A., Welsh, W. F., et al. 2016, ApJ, 827, 86 [NASA ADS] [CrossRef] [Google Scholar]
 Kostov, V. B., Orosz, J. A., Feinstein, A. D., et al. 2020, AJ, 159, 253 [Google Scholar]
 Lyra, W., Richert, A. J. W., Boley, A., et al. 2016, ApJ, 817, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, D. V., & Fabrycky, D. C. 2021, AJ, 162, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, D. V., & Fabrycky, D. 2022, BAAS [Google Scholar]
 Martin, R. G., & Lubow, S. H. 2017, ApJ, 835, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, R. G., & Lubow, S. H. 2018, MNRAS, 479, 1297 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, D. V., & Triaud, A. H. M. J. 2014, A&A, 570, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Martin, D. V., Triaud, A. H. M. J., Udry, S., et al. 2019, A&A, 624, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marzari, F., Thébault, P., & Scholl, H. 2008, ApJ, 681, 1599 [NASA ADS] [CrossRef] [Google Scholar]
 Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [Google Scholar]
 Min, M., Dullemond, C. P., Dominik, C., et al. 2009, A&A, 497, 155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mutter, M. M., Pierens, A., & Nelson, R. P. 2017, MNRAS, 465, 4735 [NASA ADS] [CrossRef] [Google Scholar]
 Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012a, Science, 337, 1511 [Google Scholar]
 Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012b, ApJ, 758, 87 [Google Scholar]
 Papaloizou, J. C. B. 2005, A&A, 432, 743 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paardekooper, S.J., Leinhardt, Z. M., Thébault, P., & Baruteau, C. 2012, ApJ, 754, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Penzlin, A. B. T., Kley, W., & Nelson, R. P. 2021, A&A, 645, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pierens, A., & Nelson, R. P. 2013, A&A, 556, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pierens, A., & Nelson, R. P. 2023, A&A, 670, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pierens, A., McNally, C. P., & Nelson, R. P. 2020, MNRAS, 496, 2849 [CrossRef] [Google Scholar]
 Pierens, A., Nelson, R. P., & McNally, C. P. 2021, MNRAS, 508, 4806 [CrossRef] [Google Scholar]
 Robitaille, T. P. 2010, A&A, 520, A70 [CrossRef] [EDP Sciences] [Google Scholar]
 Socia, Q. J., Welsh, W. F., Orosz, J. A., et al. 2020, AJ, 159, 94 [Google Scholar]
 Schwamb, M. E., Orosz, J. A., Carter, J. A., et al. 2013, ApJ, 768, 127 [Google Scholar]
 Shadmehri, M., & Khajenabi, F. 2015, MNRAS, 447, 1439 [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Socia, Q. J., Welsh, W. F., Orosz, J. A., et al. 2020, AJ, 159, 94 [Google Scholar]
 Standing, M. R., Sairam, L., Martin, D. V., et al. 2023, NatAs, 7, 702 [NASA ADS] [Google Scholar]
 Thun, D., Kley, W., & Picogna, G. 2017, A&A, 604, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Triaud, A. H. M. J., Standing, M. R., Heidari, N., et al., 2022, MNRAS, 511, 3561 [NASA ADS] [CrossRef] [Google Scholar]
 Turpin, G. A., & Nelson, R. P. 2024, MNRAS, 528, 7256 [NASA ADS] [CrossRef] [Google Scholar]
 Vartanyan, D., Garmilla, J. A., & Rafikov, R. R. 2016, ApJ, 816, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Welsh, W. F., Orosz, J. A., Carter, J. A., et al. 2012, Nature, 481, 475 [Google Scholar]
 Welsh, W. F., Orosz, J. A., Short, D. R., et al. 2015, ApJ, 809, 26 [Google Scholar]
 Zhu, Z., Hartmann, L., Gammie, C., et al. 2009, ApJ, 701, 620 [CrossRef] [Google Scholar]
All Tables
Binary stellar parameters for known CBP hosts taken from Georgakarakos et al. (2021).
All Figures
Fig. 1 α viscous parameter as a function of distance from the binary. The dependence of α on the radius was chosen so as to reproduce the α profile found in Pierens et al. (2021; see their Fig. 4), where hydrodynamical turbulence operates. 

In the text 
Fig. 2 Evolutionary tracks on the HR diagram for different circumbinary systems. The upper dashed line shows the stellar isochrone at 1 Myr, whereas the lower dashed line corresponds to the ZAMS. Isochrones were computed using MIST evolutionary models (Choi et al. 2016). 

In the text 
Fig. 3 Surface density, midplane temperature, and aspect ratio profiles of the different circumbinary systems that we considered. Solid lines correspond to models that include viscous heating, whereas dashed lines are for models without viscous heating. All models include irradiation and shock heating. In the surface density profiles, the markers indicate the location of the ice line, which is also represented as a horizontal dashed line in the temperature plots. The dashed black lines in the surface density plots show the values of α in the viscously heated models – see the righthand vertical axis. 

In the text 
Fig. 4 Temperature distributions of each of the irradiated and viscously heated models that we considered. 

In the text 
Fig. 5 Contribution of the different heating terms to the rate of specific internal change for Kepler34. 

In the text 
Fig. 6 Surface density distribution for each of the simulations. The location of the ice line is indicated by the green circle, while the locations of the stars are represented as white circles. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.