Introduction
Almost since the first detection of extrasolar planets (exoplanets, Mayor and Queloz Reference Mayor and Queloz1995), extrasolar moons (exomoons) have been the subject of intense scientific inquiry. In the Solar System, moons act as tracers of planet formation and evolution. In some cases, moons such as Europa, Enceladus and Ganymede may even host subsurface habitats (Melosh et al. Reference Melosh, Ekholm, Showman and Lorenz2004; Iess et al. Reference Iess, Stevenson, Parisi, Hemingway, Jacobson, Lunine, Nimmo, Armstrong, Asmar, Ducci and Tortora2014; Thomas et al. Reference Thomas, Tajeddine, Tiscareno, Burns, Joseph, Loredo, Helfenstein and Porco2015; Saur et al. Reference Saur, Duling, Roth, Jia, Strobel, Feldman, Christensen, Retherford, McGrath, Musacchio, Wennmacher, Neubauer, Simon and Hartkorn2015). It is quite possibly the case that subsurface-habitable moons could greatly outnumber habitable planets in the Milky Way (Scharf Reference Scharf2006; Heller and Pudritz Reference Heller and Pudritz2015). It is also quite possibly the case that exomoons are massive enough to host a substantial atmosphere like the Earth, and possess a similar biosphere.
Earth-like exomoons have been thought to be unlikely based on Solar System evidence. The moons of the giant planets possess a satellite-to-planet mass ratio 
$\epsilon \lessapprox10^{-4}$, which is consistent with models of moon formation in a circumplanetary disc (Mosqueira and Estrada Reference Mosqueira and Estrada2003a,Reference Mosqueira and Estradab; Ward and Canup Reference Ward and Canup2010; Canup and Ward Reference Canup and Ward2006). For a Jupiter-mass planet to host an Earth-like exomoon, this would require ε ≈ 3 × 10−3, i.e. at least an order of magnitude higher.
Recently, Teachey et al. (Reference Teachey, Kipping and Schmitt2017) described evidence pointing to a candidate exomoon in orbit of the Jupiter-sized transiting exoplanet Kepler-1625b. The host star, Kepler-1625, is G type, of approximately one solar mass and has recently evolved off the main sequence (Berger et al. Reference Berger, Huber, Gaidos and van Saders2018).
The exomoon candidate was observed in three transits of Kepler-1625b over the 4-year primary mission of the Kepler Space Telescope, out of a maximum of five transits (the orbital period of Kepler-1625b being 287 days). Transit T2 appears to indicate the satellite performing an early ingress of the transit, with T4 showing the satellite in late egress. Interestingly, T5 shows both early ingress and late egress, where the moon lags some 10 h behind the planet in exiting the stellar disc.
Amongst other measurements, this gave initial constraints on not only the projected separation of the two bodies, but also the orbital period of the moon – around P ps ~ 72 h, and a separation of around 19 planetary radii (R P). Hence by Kepler's third law, the barycentric mass for the planet moon component was derived to be around 
$17.6^{+2.1}_{-1.9}\, M_{\rm Jup}$.
The mass for Kepler-1625b was not well-constrained, and hence the exomoon candidate (Kepler-1625b-i) also has a poorly constrained mass.Footnote 1 Teachey et al. (Reference Teachey, Kipping and Schmitt2017) suggested that the planet is 10 M Jup, with the moon being around 17 
$M_{\rm \oplus}$, giving ε = 5.34 × 10−3. Heller (Reference Heller2017) notes that the barycentric mass could be shared differently than Teachey et al. (Reference Teachey, Kipping and Schmitt2017)'s description. Indeed, Kepler-1625b may be a brown dwarf or a low mass star, and the exomoon candidate mass would be closer to 1
$M_{\rm \oplus }$.
Subsequent observation with the Hubble Space Telescope (HST) found further evidence in favour of the exomoon's presence (Teachey and Kipping Reference Teachey and Kipping2018), favouring the 10M Jup planet, Neptune-mass exomoon interpretation. However, this inference is now dominated by the single HST observation, and the authors advocate further monitoring to confirm the moon-like signal.
In any case, if this candidate is confirmed by subsequent observations, this would indicate that satellite systems with ε ~ 10−3 can indeed exist, and that massive exomoons with Earth-like properties are indeed possible.
Massive exomoons raise the stakes for habitability even further. One-dimensional (1D) climate modelling of Earth-like exomoons in the orbit of giant planets around Sun-like stars show that the habitable zone for such objects occupies a significant volume of orbital parameter space (Forgan and Kipping Reference Forgan and Kipping2013; Heller and Barnes Reference Heller and Barnes2013). This parameter space is generally larger in volume than that of Earth-like exoplanets due to the additional sources and sinks of radiation.
Planetary illumination allows the habitable zone to move further away from the host star, especially if the planet is sufficiently large that it is significantly self-luminous in thermal radiation (Heller and Barnes Reference Heller and Barnes2013). Typically, the thermal flux dominates over reflected starlight, suggesting a slowly varying illumination, but it is worth noting that moons on eccentric orbits will still experience strong variations in illumination (Hinkel and Kane Reference Hinkel and Kane2013).
On the other hand, moons orbiting close to their host planet are more likely to experience relatively long eclipses of the star by the planet, which can result in a net loss of radiation (Heller Reference Heller2012). While the loss from a single eclipse can usually be buffered by the thermal inertia of an Earth-like atmosphere, if the eclipses are frequent and sufficiently long in duration, moons can enter a snowball state which they cannot exit even after leaving the eclipse (Forgan and Yotov Reference Forgan and Yotov2014).
Tidal heating can allow the moon's habitability to be almost independent of the star, but by the same token can be strong enough to render the moon uninhabitable. The complex nature of tidal heating demands careful modelling of the tidal force, its resultant deformation of the moon's interior, and the release of this stress as heat (Dobos and Turner Reference Dobos and Turner2015; Forgan and Dobos Reference Forgan and Dobos2016).
In this letter, we compute the habitable zone for Earth-like exomoons orbiting Kepler-1625b assuming Teachey and Kipping (Reference Teachey and Kipping2018)'s derived parameters, and our exomoon climate models (Forgan and Kipping Reference Forgan and Kipping2013; Forgan and Yotov Reference Forgan and Yotov2014; Forgan and Dobos Reference Forgan and Dobos2016). We briefly describe the three versions of the climate model used to compute exomoon habitability in section ‘Methods’, discuss the resulting habitable zones in section ‘Results and Discussion’, and conclude in section ‘Conclusions’.
Methods
Simulation setup
We fix the stellar and planetary parameters, in accordance with Teachey and Kipping (Reference Teachey and Kipping2018) as follows. The star mass 
$M_\ast= 1.04 M_{\rm \odot }$. The planet's mass M p = 10M Jup, and has radius R p = 1.015 R Jup. The planet's orbital semi-major axis and eccentricity are also fixed at a p = 0.98 and e p = 0, respectively. This fixes the planet's Hill radius
$$R_{\rm H,p} = a_{\rm p} \left({M_{\rm p}\over 3M_{\ast}}\right)^{1/3} =0.1422\, \rm{AU} \approx 293.7 \,R_{\rm p}.$$ We assume that the moon is Earth-like (
$M_{\rm s} = 1M_{\rm \oplus}$, 
$R_{\rm s} = 1R_{\rm \oplus}$). We further assume that the planet resides at the barycentre of the moon-planet system, which is satisfactory given that ε ≪ 1. The inclination of the planet relative to the stellar equator, i p = 0 (i.e. the planet orbits in the x−y plane).
The inclination of the moon relative to the planet's equator, i.e. the inclination of the moon relative to the x−y plane, i m, is also zero (unless otherwise stated). The orbital longitudes of the planet and moon are defined such that ϕ p = ϕ m = 0 corresponds to the x-axis. We also assume that the moon's obliquity has been efficiently damped by tidal evolution, and we therefore set it to zero.
We allow the moon's semi-major axis and eccentricity (a m and e m, respectively) to vary. For each model, we run 500 simulations, where the lunar semi-major axis takes a range of values: a m = (0.05, 0.3) R H,p, and eccentricities range from e m = (0, 0.08). Note that the exomoon candidate Kepler-1625b-i has an estimated orbital separation of a m,i = 0.153 R H,p.
Latitudinal energy balance modelling
We compute the habitability of Earth-like exomoons using the Latitudinal Energy Balance Model (LEBM) approach (North and Coakley Reference North and Coakley1979; North et al. Reference North, Mengel and Short1983).
The core equation of the LEBM is a diffusion equation, solved over latitude λ ∈ ( −90°, 90°). In practice, we solve the equation using the convenience variable x ≡ sin λ:
$$ \hskip-1pc C {\partial T \over \partial t} - {\partial \over \partial x}\left(D(1-x^2) {\partial T \over \partial x}\right) = (S+S_{\rm p})\left[1-A(T)\right] + \zeta - I(T),$$where T = T(x, t) is the temperature at time t, and the boundary condition dT/dx = 0 at the poles.
C is the atmospheric heat capacity, the diffusion coefficient D controls latitudinal heat redistribution, S and S p are the stellar and planetary illumination, respectively, ζ is the surface heating generated by tides in the moon's interior, I is the atmospheric infrared cooling and A is the albedo.
We produce five measures of the habitable zone, using three different versions of the climate model, as described in Table 1. Essentially, these represent increasing realism for the energy balance model, as we add in the effect of planetary illumination, and the carbonate-silicate cycle. They also consider the current luminosity of Kepler-1625 as it evolves off the main sequence, as well as its previous luminosity while on the main sequence.
Table 1. The three-model runs in this paper
The control runs (CN).
In the first set of three-model runs (corresponding to the model used in Forgan and Kipping (Reference Forgan and Kipping2013)), we use the following prescriptions.
The atmospheric heat capacity C depends on what fraction of the moon's surface is ocean, f ocean = 0.7, what fraction is land f land = 1.0 − f ocean, and what fraction of the ocean is frozen f ice:
$$C = f_{\rm land}C_{\rm land} + f_{\rm ocean}\left[(1-f_{\rm ice})C_{\rm ocean} + f_{\rm ice} C_{\rm ice}\right].$$The heat capacities of land, ocean and ice-covered areas are
$$C_{\rm land} = 5.25 \times 10^9 \rm{erg \, cm^{-2} K^{-1}},$$
$$C_{\rm ocean} = 40.0C_{\rm land},$$
$$C_{\rm ice} = \left \{ \matrix{9.2C_{\rm land} \hfill & 263 {\rm K} \lt T \lt 273 {\rm K} \hfill\cr 2C_{\rm land} \hfill & T \lt 263 {\rm K} \hfill\cr 0.0 & T \gt 273 {\rm K}.\hfill}.\right.$$These parameters assume a wind-mixed ocean layer of 50m (Williams and Kasting Reference Williams and Kasting1997). Increasing the assumed depth of this layer would increase C ocean (see e.g. North et al. (Reference North, Mengel and Short1983) for details). The albedo function is
$$A(T) = 0.525 - 0.245 \tanh \left[{T-268\, {\rm K} \over 5\, {\rm K}} \right].$$This produces a rapid shift from low albedo (~0.3) to high albedo (~0.75) as the temperature drops below the freezing point of water, producing highly reflective ice sheets. Figure 1 of Spiegel et al. (Reference Spiegel, Menou and Scharf2008) demonstrates how this shift in albedo affects the potential for global energy balance, and that for planets in circular orbits, two stable climate solutions arise, one ice-free, and one ice-covered. Spiegel et al also show that such a function is sufficient to reproduce the annual mean latitudinal temperature distribution on the Earth.
Fig. 1. The habitable zone for exomoons orbiting Kepler-1625b, for the three control model runs CN. Top: The habitable zone for i m = 0°, assuming stellar luminosity consistent with the main sequence (CN0). Middle: The same, but for i m = 45° (CN45). Bottom: The habitable zone for i m = 0, using estimates of Kepler-1625b's current luminosity, rather than that derived from main sequence relations (CNL). Note that the candidate exomoon Kepler-1625b-i orbits at 0.153 Hill radii. If Kepler-1625b-i exists, then any other exomoon orbiting with a m < 0.17 Hill radii is likely to be dynamically unstable.
The diffusion constant D is calibrated such that a fiducial Earth–Sun climate system reproduces Earth's observed latitudinal temperature (see e.g. North et al. (Reference North, Cahalan and Coakley1981); Spiegel et al. (Reference Spiegel, Menou and Scharf2008)). Planets that rotate rapidly experience inhibited latitudinal heat transport, due to Coriolis forces truncating the effects of Hadley circulation (cf Farrell (Reference Farrell1990); Williams and Kasting (Reference Williams and Kasting1997)). The partial pressure of CO2 also plays a role. We follow Williams and Kasting (Reference Williams and Kasting1997) by scaling D according to:
$$D=5.394 \times 10^2 \left({\omega_{\rm d} \over \omega_{\rm d,\oplus}}\right)^{-2} \left({P_{\rm{CO_2}} \over P_{\rm{CO_2},\oplus}}\right),$$where ω d is the rotational angular velocity of the planet, and 
$\omega _{\rm d,\oplus}$ is the rotational angular velocity of the Earth, and 
$P_{\rm {CO_2},\oplus }=3.3 \times 10^{-4} \,\rm {bar}$. In this run, the partial pressure of CO2 is fixed: 
$P_{\rm {CO_2}} = P_{\rm {CO_2},\oplus}$.
The stellar insolation flux S is a function of both season and latitude. At any instant, the bolometric flux received at a given latitude at an orbital distance r is
$$S = q_0 \left({M \over M_{\rm \odot}}\right)^4\cos Z \left({1 AU \over r}\right)^2,$$where q 0 is the bolometric flux received from a 
$1M_{\rm \odot }$ star at a distance of 1 AU, and we have assumed a standard main sequence luminosity relation. Z is the zenith angle:
$$q_0 = 1.36\times 10^6\left({M_{\ast} \over M_{\rm \odot}}\right)^4 {\rm erg} \,{\rm s}^{-1}\, {\rm cm}^{-2}$$
$$\cos Z = \mu = \sin \lambda \sin \delta + \cos \lambda \cos \delta \cos h.$$δ is the solar declination, and h is the solar hour angle. As stated previously, we set the moon's obliquity δ 0 to zero. The solar declination is calculated as:
$$\sin \delta = - \! \sin \delta_0 \cos(\phi _{\ast{\rm m}}-\phi_{\rm peri,m}-\phi_{\rm a}),$$where ϕ 
$ _{\ast{\rm m}}$ is the current orbital longitude of the moon relative to the star, ϕ peri,m is the longitude of periastron, and ϕ a is the longitude of winter solstice, relative to the longitude of periastron. We set ϕ peri,m = ϕ a = 0 for simplicity.
We must diurnally average the solar flux:
$$S = q_0 \bar{\mu}.$$This means we must first integrate μ over the sunlit part of the day, i.e. h = [ −H, +H], where H is the radian half-day length at a given latitude. Multiplying by the factor H/π (as H = π if a latitude is illuminated for a full rotation) gives the total diurnal insolation as
$$S = q_0 \left({H \over \pi}\right) \bar{\mu} = {q_0 \over \pi} \left(H \sin \lambda \sin \delta + \cos \lambda \cos \delta \sin H\right).$$The radian half day length is calculated as
$$\cos H = - \! \tan \lambda \tan \delta.$$We allow for eclipses of the moon by the planet (where the insolation S = 0). We detect an eclipse by computing the angle α between the vector connecting the moon and planet, s, and the vector connecting the moon and the star 
$ {\bi s}_*$:
$$\cos \alpha = \hat{\bi s}.\hat{\bi s}_{\ast}.$$It is straightforward to show that an eclipse is in progress if
$$\left|s_{\ast}\right| \sin \alpha < R_{\rm p}.$$We do not model the eclipse ingress and egress, and instead simply set S to zero at any point during an eclipse. A typical eclipse duration in these runs is approximately 6h (for an exomoon in a circular orbit around Kepler-1625b at a m = 0.1 RH). Our simulation timestep includes a condition to ensure that any eclipse must be resolved by at least ten simulation timesteps.
We use the following infrared cooling function:
$$ I(T) = {\sigma_{SB}T^4 \over 1 +0.75 \tau_{IR}(T)},$$where the optical depth of the atmosphere
$$\tau_{IR}(T) = 0.79\left({T \over 273\,{\rm K}}\right)^3.$$Tidal heating is calculated by assuming the tidal heating per unit area is (Peale et al. Reference Peale, Cassen and Reynolds1980; Scharf Reference Scharf2006):
$$\left({{\rm d}E \over {\rm d}t}\right)_{\rm tidal} = {21 \over 38}{\rho^2_{\rm m} R^{5}_{\rm m} e^2_{\rm m} \over \Gamma Q}\left({GM_{\rm p} \over a^3_{\rm m}}\right)^{5/2},$$where Γ is the moon's elastic rigidity (which we assume to be uniform throughout the body), R m is the moon's radius, ρ m is the moon's density, M p is the planet mass, a m and e m are the moon's orbital semi-major axis and eccentricity (relative to the planet), and Q is the moon's tidal dissipation parameter. We assume terrestrial values for these parameters, hence Q = 100, Γ = 1011 dyne cm−2; (appropriate for silicate rock) and ρ m = 5 g cm−3;.
We assume that this heating occurs uniformly across the moon's surface, which is a large approximation but a necessity of 1D LEBM models.
The planetary illumination S p is set to zero for all three control runs. Two of our three runs consider the habitable zone for i m = 0° (CN0) and i m = 45° (CN45), assuming the stellar luminosity is set by main sequence relations, giving
$${L \over L_{\rm \odot}}= \left({M_{\ast} \over M_{\rm \odot}}\right)^4 = 1.16.$$These runs consider the main sequence phase of Kepler-1625b. However, as previously stated Kepler-1625b is now evolving off the main sequence. Therefore, in the last of the control runs, we reset i m = 0° and replace the main sequence luminosity with an estimate of Kepler-1625b's current luminosity, which is calculated assuming a blackbody and an effective temperature of T eff = 5, 548 K (Mathur et al. Reference Mathur, Huber, Batalha, Ciardi, Bastien, Bieryla, Buchhave, Cochran, Endl, Esquerdo, Furlan, Howard, Howell, Isaacson, Latham, MacQueen and Silva2017). This results in a much greater luminosity of 
$L = 2.68 L_{\rm \odot }$.
Planetary illumination (IL).
In this run, we use the same inputs as CN (with i m = 0 and main sequence luminosity), but now implement planetary illumination and eclipses of the moon by the planet. We use the same prescriptions as Forgan and Yotov (Reference Forgan and Yotov2014) (see also Heller and Barnes (Reference Heller and Barnes2013)). Illumination adds a second insolation source to the system, dependent on both the starlight reflected by the planet, and on the planet's own thermal radiation.
We assume the planet's orbit is not synchronous, and that the temperature of the planet T p is uniform across the entire surface. We also fix the planetary albedo α p = 0.3. The insolation due to the planet is therefore
$$S_{\rm p}(t)=f_{\rm t}(t)+f_{\rm r}(t).$$The thermal flux f t as a function of latitude λ is:
$$f_{\rm t}(\lambda,t)= {2R_{\rm p}^2 \sigma_{SB} \over a_{\rm m}^2} (\cos \lambda) T_{\rm p}^4,$$And the reflected flux f r is:
$$f_{\rm r}(\lambda,t)={2L_{\ast} \over 4\pi r_{{\rm p}{\ast}}^2} {R_{\rm p}^2 \pi \alpha_{\rm p} \over a_{\rm m}^2} \cos { \lambda} {\islant -20% \Xi} (t).$$Where Ξ is the fraction of dayside visible from the lunar surface (see Forgan and Yotov (Reference Forgan and Yotov2014) for more details). Calculating the ratio f r/f t for Kepler-1625b shows that thermal flux dominates the planetary illumination budget, with around 1% of the total illumination being produced by reflected starlight. We should therefore expect S p to be roughly constant over the moon's orbit (as T p is uniform).
The carbonate-silicate-cycle (CS).
In this final run, we now use a simple piecewise function to determine 
$P_{\rm {CO_2}}$ as a function of local temperature (Spiegel et al. Reference Spiegel, Raymond, Dressing, Scharf and Mitchell2010):
$$P_{\rm{CO_2}} = \left \{ \matrix{ 10^{-2} \, \rm{bar} \hfill& T \leq 250 \, {\rm K} \hfill \cr 10^{-2 - (T-250)/27} \,\rm{bar} \hfill& 250 \,{\rm K} \lt T \lt 290 \,{\rm K} \hfill \cr P_{\rm{CO_2}, \oplus} \hfill& T \geq 290 \, {\rm K}. }\right.$$Our prescription now allows D to vary with latitude, depending on the local temperature. This is not guaranteed to produce Hadley circulation (see e.g. Vladilo et al. (Reference Vladilo, Murante, Silva, Provenzale, Ferri and Ragazzini2013) for details on how D can be modified to achieve this). As we allow partial pressure of CO2 to vary, we now replace equation (18) with Williams and Kasting (Reference Williams and Kasting1997)'s prescription for the cooling function, 
$I(T,P_{\rm {CO_2}})$:
$$\eqalign{I = & \,9.468980 -7.714727 \times 10^{-5} \beta - 2.794778T \cr &- 3.244753 \times 10^{-3} \beta T -3.4547406 \times 10^{-4}\beta^2 \cr &+ 2.212108 \times 10^{-2} T^2 + 2.229142 \times 10^{-3} \beta^2 T \cr &+ 3.088497 \times 10^{-5} \beta T^2 - 2.789815 \times 10^{-5} \beta^2 T^2 \cr &- 3.442973 \times 10^{-3} \beta^3 - 3.361939 \times 10^{-5} T^3 \cr &+ 9.173169 \times 10^{-3} \beta^3 T - 7.775195 \times 10^{-5} \beta^3 T^2 \cr &- 1.679112 \times 10^{-7} \beta T^3 + 6.590999 \times 10^{-8} \beta^2 T^3 \cr & +1.528125 \times 10^{-7} \beta^3 T^3 - 3.367567 \times 10^{-2} \beta^4 \cr & -1.631909 \times 10^{-4} \beta^4 T + 3.663871 \times 10^{-6} \beta^4 T^2 \cr & -9.255646 \times 10^{-9} \beta^4 T^3}$$where we have defined
$$\beta = \log \left({P_{\rm{CO_2}} \over P_{\rm{CO_2},\oplus}}\right).$$Simulation timestep
The diffusion equation is solved using a simple explicit forward time, centre space finite difference algorithm. A global timestep was adopted, with constraint
$$\delta t < {\left(\Delta x\right)^2C \over 2D(1-x^2)}.$$This timestep constraint ensures that the first term on the left-hand side of equation (2) is always larger than the second term, preventing the diffusion term from setting up unphysical temperature gradients. The parameters are diurnally averaged, i.e. a key assumption of the model is that the moons rotate sufficiently quickly relative to their orbital period around the primary insolation source. This is broadly true, as the star is the principal insolation source, and the moon rotates relative to the star on timescales of a few days.
Mapping the exomoon habitable zone.
For each run, we simulate climate models over a range of lunar semi-major axes a m and eccentricities e m, mapping out the habitable zone as a function of (a m, e m) by classifying the resulting climate according to its habitability function ξ:
$$\xi(\lambda,t) = \left\{ \matrix{1 \hfill& 273 \,\hbox{K} \lt T ( \lambda , t) \lt 373 \,\hbox{K} \hfill \cr 0 \hfill & \hbox{otherwise}.\hfill} \right.$$We average this over latitude to calculate the fraction of habitable surface at any timestep:
$$\xi(t) = {1 \over 2} \int_{-\pi/2}^{\pi/2}\xi(\lambda,t)\cos \lambda \, {\rm d}\lambda.$$Each simulation evolved until it reaches a steady or quasi-steady state, and the final 10 years of climate data are used to produce a time-averaged value of ξ(t), 
$\bar {\xi }$. Along with the sample standard deviation, σξ, we can classify each simulation as follows:
1. Habitable moons - these moons possess a time-averaged
$ \bar {\xi } \gt 0.1$, and $\sigma _{\xi } \lt 0.1\bar {\xi }$, i.e. the fluctuation in habitable surface is less than 10% of the mean.2. Hot moons - these moons have average temperatures above 373 K across all seasons, and are therefore conventionally uninhabitable, and
$\bar {\xi } \lt 0.1$.3. Snowball moons - these moons have undergone a snowball transition to a state where the entire moon is frozen, and are therefore conventionally uninhabitableFootnote 2.
4. Transient moons - these moons possess a time-averaged
$\bar {\xi } \gt 0.1$, and $\sigma _{\xi } \gt 0.1\bar {\xi }$, i.e. the fluctuation in habitable surface is greater than 10% of the mean.
Results and discussion
Control runs
In Figure 1, we display results for the CN0, CN45 and CNL runs. For CN0 and CN45, orbits with a semi-major axis above 0.1 R H,p are habitable for a wide range of eccentricities. There is no outer circumplanetary habitable edge, despite all bodies rotating in the same plane, maximising the effect of eclipses. The effect of inclination results in some moons near the habitable zone inner boundary experiencing large temperature variations, but otherwise maintaining a habitable surface.
Notably, there are no habitable moons when Kepler-1625b's luminosity is increased from its main sequence value to its current estimated value. Extra runs considering larger a m also fail to find any habitable solutions.
Planetary illumination and the carbonate-silicate cycle
In Figure 2, we show the results from the IL and CS runs (with i m = 0). Planetary illumination makes little appreciable impact on the habitable zone, which is consistent with previous calculations (Forgan and Yotov Reference Forgan and Yotov2014). Ideally, as the planetary illumination is typically in the infrared, it should be subject to a different albedo than the stellar illumination (cf Heller and Pudritz Reference Heller and Pudritz2015). Experimenting with different albedos for S p make little appreciable difference to the results.
Fig. 2. The habitable zone for exomoons orbiting Kepler-1625b, for model runs IL (top) and CS (bottom). Note that the candidate exomoon Kepler-1625b-i orbits at 0.153 Hill radii. If Kepler-1625b-i exists, then any other exomoon orbiting with a m < 0.17 Hill radii is likely to be dynamically unstable.
When CO2 partial pressure is allowed to vary (run CS), a runaway greenhouse takes effect in all runs, preventing the moon from sustaining a habitable surface. Further runs have demonstrated that this is independent of tidal heating, and that extending the parameter sweep to larger a m does not yield a habitable solution.
Orbital stability of moons due to Kepler-1625b-i
Moons remain orbitally stable for semi-major axes a m < 0.3 − 0.5 R H,p (Domingos et al. Reference Domingos, Winter and Yokoyama2006), and we can therefore expect all moons simulated in this work to remain on stable orbits for long timescales, provided they do not orbit too close to Kepler-1625b-i (if it exists). Taking Teachey and Kipping (Reference Teachey and Kipping2018)'s estimate of Kepler-1625b-i as approximately Neptune-mass 
$\approx 17M_{\rm \oplus}$, then we should expect an inner orbital stability limit defined roughly by the Hill radius of the satellite candidate:
$$R_{\rm H,i} = a_{\rm m,i} \left({M_{\rm s} \over 3 M_{\rm p}}\right)^{1/3} = 0.019 R_{\rm H,p},$$where a m,i = 0.153 R H,p. Any satellite which orbits within a few mutual Hill radii of the candidate will be subject to dynamical instability. We can therefore expect that any Earth-like exomoon orbit will only be stable provided 
$a_{\rm m} \gtrapprox 0.17-0.2 R_{\rm H,p}$. This still leaves a large range of orbits for Earth-like exomoons that are both habitable (depending on the input climate model) and dynamically stable.
We note that while many Earth-like exomoons can orbit stably in the Kepler-1625b system despite the presence of Kepler-1625b-i, we have not computed the cyclic variations in orbital eccentricity we might expect as a Neptune-sized body exerts its gravitational field on a neighbour. The large mass and inclination of Kepler-1625b-i is likely to drive strong gravitational perturbations. If the eccentricity of a body is driven too high, tidal heating could rapidly render an Earth-like exomoon uninhabitable. We are currently running climate calculations using OBERON (Forgan Reference Forgan2016a) that includes the gravitational interaction between bodies to explore this further.
Clouds are not explicitly considered in this model (although they are implicitly accounted for in run CS). Clouds can modify both the albedo and optical depth of the system significantly. Also, we assume that both stellar and planetary flux are governed by the same lunar albedo, which in truth is not likely to be the case (see e.g. Heller and Barnes (Reference Heller and Barnes2015)). Planetary flux at infrared wavelengths is more easily absorbed by ice sheets and produces more efficient melting (see e.g. Shields et al. Reference Shields, Meadows, Bitz, Pierrehumbert, Joshi and Robinson2013). However, our simulations do not produce large quantities of ice in any of the model runs. It is possible that more efficient absorption of IR radiation might move the inner habitable zone further outwards in a m, but this can be equally offset by cloud cover.
In several runs, we find that an Earth-like exomoon at the calculated location of Kepler-1625b-i would possess a habitable climate. Given that Kepler-1625b-i appears to be a relatively massive exomoon, we can consider the possibility that this exomoon candidate could have its own satellite (a ‘moon–moon’). If we assume an Earth-like satellite of the exomoon candidate, we can compute the minimum and maximum permitted orbital semi-major axies. The Hill radius of the satellite candidate is approximately 5 R p, giving an outer stability limit of approximately 2.5R p ≈ 5.7R s (Domingos et al. Reference Domingos, Winter and Yokoyama2006). The inner stability limit is defined by tidal disruption. Setting 
$M_{\rm ss}=1M_{\rm \oplus}$ to be the mass of the moon–moon, and 
$R_{\rm ss}=1R_{\rm \oplus}$, we can compute the Roche radius:
$$R_{\rm Roche, ss} = R_{\rm ss}\left({M_{\rm s} \over M_{\rm ss}}\right)^{1/3} = 0.5 R_{\rm s}$$There therefore exists a reasonable range of distances, between say [2, 5.7]R s, at which an Earth-like moon–moon could orbit the exomoon candidate. If this moon–moon was Earth-like, our models suggest it would have been habitable during the star's main sequence phase, and its climate would be quite analogous to the climates of planets in binary systems (cf Forgan (Reference Forgan2012); Kaltenegger and Haghighipour (Reference Kaltenegger and Haghighipour2013)).
Conclusions
We have applied several different exomoon climate models to the Kepler-1625b system, to investigate the morphology of the exomoon habitable zone around this Jupiter-radius object.
We find that for a range of assumptions, Earth-like exomoons were likely to have been habitable while Kepler-1625 occupied the main sequence, for a wide range of orbital semi-major axes and eccentricities. These exomoons would remain orbitally stable even if the exomoon candidate Kepler-1625b-i is indeed present. However, as Kepler-1625 evolved off the main sequence, the luminosity increased to levels that generally destroy the habitability of any Earth-like exomoons possibly present.
Exomoon detection remains a challenging observational endeavour – as such, the ability to determine the habitability of detected exomoons is an even greater challenge. The four classes of exomoon climate displayed in this work will be extremely difficult to distinguish between. The best approaches will require some form of spectroscopic data to assess atmospheric composition and thermal state. This could be obtained if the moon is sufficiently bright compared with the planet at wavelengths of interest. Proposed techniques for obtaining exomoon spectral data involve spectroastrometry (Agol et al. Reference Agol, Jansen, Lacy, Robinson and Meadows2015), oscillations in the combined exoplanet–exomoon phase curve (Forgan Reference Forgan2017), or radial velocity measurements of the exoplanet using high dispersion spectroscopy (Brogi & Forgan, in prep.).
If Kepler-1625b-i is real, it is likely massive enough to possess its own satellite (a ‘moon–moon'). If the said moon–moon was Earth-like, it could have resided in a moon–moon habitable zone during Kepler-1625's main sequence phase. The morphology of moon–moon habitable zones are not yet explored, but will share similarities with that of S-type binary star systems (Kaltenegger and Haghighipour Reference Kaltenegger and Haghighipour2013; Cuntz Reference Cuntz2014; Forgan Reference Forgan2016b).
As an aside, we should note that 1D calculations are now giving way to full 3D global circulation models of exomoon atmospheres (Haqq-Misra and Heller Reference Haqq-Misra and Heller2018). The 3D aspect of exomoon climates is crucial, as planetary illumination heats the top of the atmosphere, and tidal effects heat the surface, resulting in complex heat redistribution patterns. For example, planetary illumination amplifies warming at the moon's poles, an effect not seen in 1D calculations. We look forward to further 3D studies of exomoon atmospheres that explore the habitable zone's orbital parameters as defined by studies such as this work.
We suggest that future studies of habitability of rocky worlds should continue to explore what have until now been considered rather unusual regimes in parameter space, as it seems likely the Universe will continue to deliver surprising configurations for celestial objects.
Author ORCIDs
Duncan H. Forgan, 0000-0003-1175-4388
Acknowledgments
The author gratefully acknowledges support from the ECOGAL project, grant agreement 291227, funded by the European Research Council under ERC-2011-ADG. The author warmly thanks the reviewer for their comments and suggestions. This research has made use of NASA's Astrophysics Data System Bibliographic Services.
