Dynamical Effects on the Habitable Zone for Earth-like Exomoons

With the detection of extrasolar moons (exomoons) on the horizon, it is important to consider their potential for habitability. If we consider the circumstellar Habitable Zone (HZ, often described in terms of planet semi-major axis and orbital eccentricity), we can ask,"How does the HZ for an Earth-like exomoon differ from the HZ for an Earth-like exoplanet?"For the first time, we use 1D latitudinal energy balance modelling to address this question. The model places an Earth-like exomoon in orbit around a Jupiter mass planet, which in turn orbits a Sun-like star. The exomoon's surface temperature is evolved under the action of stellar insolation, atmospheric cooling, heat diffusion, eclipses and tidal heating. We use this model to carry out two separate investigations. In the first, four test cases are run to investigate in detail the dependence of the exomoon climate on orbital direction, orbital inclination, and on the frequency of stellar eclipse by the host planet. We find that lunar orbits which are retrograde to the planetary orbit exhibit greater climate variations than prograde orbits, with global mean temperatures around 0.1 K warmer due to the geometry of eclipses. If eclipses become frequent relative to the atmospheric thermal inertia timescale, climate oscillations become extremely small. In the second investigation, we carry out an extensive parameter study, running the model many times to study the habitability of the exomoon in the 4-dimensional space composed of the planet semi-major axis and eccentricity, and the moon semi-major axis and eccentricity. We find that for zero moon eccentricity, frequent eclipses allow the moon to remain habitable in regions of high planet eccentricity, but tidal heating severely constrains habitability in the limit of high moon eccentricity, making the habitable zone a sensitive function of moon semi-major axis.


INTRODUCTION
The increasing rate of discovery of extrasolar planets (exoplanets) has given astronomers licence to consider their potential as habitats for extraterrestrial life. To date, more than 860 exoplanets have been confirmed 1 , with over 2,300 exoplanet candidates identified by the Kepler mission (Batalha et al. 2013). Of the confirmed exoplanets, several are in the circumstellar Habitable Zone (HZ), an annular region surrounding the star where planets with orbits inside it may be expected to possess liquid water on their surface.
The HZ is derived for planets of terrestrial mass and composition -the inner edge is determined by water loss (via hydrogen escape and photolysis), and the outer edge is deter-tle viscosity, which suppresses convection and as a result has implications for heat transfer and the planet's ability to degas (Stamenkovic, Spohn & Breuer 2010;Liu & Zhong 2010).
This being the case, the current crop of exoplanets in the HZ should not be ruled out as a source of habitable niches. Extrasolar moons (exomoons) belonging to these planets may be terrestrial in nature, and as a result may be habitable (Williams, Kasting & Wade 1997). This is true even in our own Solar System: outside the Earth, the best chances for life in the Solar System are believed to be on tidally heated icy moons which are thought to possess liquid subsurface oceans, such as Europa (Melosh et al. 2004) and Enceladus (Parkinson et al. 2007). Titan, the largest moon in the Solar System, has a thick atmosphere, with evidence of precipitation and lakes constituting a hydrological cycle (Stofan et al. 2007).
Exomoons are yet to be detected, but Earth-mass exomoons are in principle detectable indirectly with Kepler (Kipping, Fossey & Campanella 2009). Studying an exoplanet's transit timing variations and transit duration variations (TTV and TDV respectively) allows observational constraints on the moon's mass, and its orbital semi-major axis around the host exoplanet. Additionally, auxiliary transits of the moon and planet-moon mutual events lead to unique transit features revealing the moon's radius and other orbital parameters (Kipping 2011). Such efforts to realise a detection are underway as part of the Hunt for Exomoons with Kepler (HEK) project (Kipping et al. 2012).
A large component of habitability studies is the study of the energy balance associated with the celestial body in question. The various sources and sinks of heat are catalogued and compared, and the system is evolved over time to understand their interplay, either by analytical calculation or numerical simulation. In the case of exomoons, the list of sources is longer than that of exoplanets. Alongside direct stellar insolation, exomoons receive starlight reflected from the planet (as well as the planet's own radiation), and they are tidally heated. Also, the moon is more prone to stellar eclipses as a result of the host planet, providing an extra (effective) energy sink.
The study of habitable exomoons has grown apace in recent years: early calculations by Reynolds, McKay & Kasting (1987) demonstrated that Europa is a viable niche for marine life such as that found under Antartic ice, and proposed a tidally heated circumplanetary habitable zone. In the early phase of exoplanet detection, Williams, Kasting & Wade (1997) suggested that two of the then 9 known giant exoplanets could host habitable moons. They raised the three key obstacles to exomoon habitability, namely: (i) Tidal locking of the moon with respect to the planet, resulting in extreme weather conditions (Joshi 1997).
(ii) Increased EUV and X-Ray radiation from the host planet's magnetosphere leading to loss of the moon's atmosphere (Kaltenegger 2000), and (iii) A poorly proportioned distribution of volatile abundances, due to the differing environments of planet and moon formation. Kaltenegger (2000) notes that Solar System moons may be overabundant with water to be habitable, if placed in the inner Solar System, as they would lack the large areas of land required to sustain a carbonate-silicate cycle. Without this, the body would lack a crucial regulatory system to ameliorate the effect of the star's growing luminosity as it evolves along the Main Sequence. Scharf (2006) investigated the then-known exoplanet population for the potential to host habitable exomoons. Estimating that almost 30% of the total population could host small moons with icy mantles, he goes on to estimate that tidal heating could increase the habitable zone's extent in planet semi-major axis by a factor of around 2 (although he notes that the tidal heating required can become several orders of magnitude larger than that experienced by Io, one of the more extreme cases in our Solar System). Heller & Barnes (2013) points out that exomoons may in some cases provide better niches than exoplanets, for several reasons: (i) Exoplanets in the HZ of low mass stars (e.g. Dressing & Charbonneau 2013) are likely to be tidally locked, whereas moons of such exoplanets are likely to be tidally locked to the planet, not the star (Kaltenegger & Selsis 2010;Sasaki, Barnes & O'Brien 2012). This reduces the strong fluctuation of surface insolation (and potential atmospheric collapse) that a moon would suffer if it was tidally locked to the star (Joshi 1997;Kite, Gaidos & Manga 2011). This may also provide energy to drive the moon's internal dynamo, and generate a magnetosphere to guard against atmospheric loss.
(ii) Exoplanet obliquity will erode due to tidal evolution, as demonstrated by the axial orientations of Venus and Mercury. This erosion is resisted more easily by more massive planets, if its orbital semi major axis is not too low and the planet does not possess relatively massive neighbours (Heller, Leconte & Barnes 2011;Sasaki, Barnes & O'Brien 2012). As a result, exomoons tidally locked to their host planet can retain the obliquity of their host planet and not suffer such strong erosion.
In their work, Heller & Barnes (2013) calculate the time dependent flux upon a moon due to the star's insolation, reflected starlight and thermal emission from the planet's surface, and add to this the tidal heating received. They use this total effective flux to calculate a parameter space of habitability in the moon's orbit about the host planet, and the host planet's mass, which they investigate for Kepler-22b and KOI211.01, both of which reside in the circumstellar HZ. By averaging over the moon's orbital period, they produce a corrective term for the circumstellar HZ depending on the planet's albedo and the moon's orbital semi-major axis relative to the planet, as well as surface temperature maps for the moons. However, while these calculations are among the most sophisticated to date, they do not allow for the transfer of heat due to an atmosphere. This is clearly a crucial component of any climate model, particularly in the case of exomoon systems, which present many opportunities to set up strong temperature gradients that will generate heat transport.
In this paper, we implement a 1D latitudinal energy balance model (LEBM) to study the climate of an Earth-like exomoon in orbit around a Jupiter mass planet, which in turn orbits a Sun-like star. Our use of a LEBM allows us to go beyond analytical calculations of detailed energy balance by allowing the moon to transfer heat across its surface (via a simple diffusion approximation) as well as adding some of the sources and sinks of energy described above.
Rather than investigate the circumstellar HZ separately from the circumplanetary habitable region, we choose instead to investigate them together as a combined four dimensional manifold, constituted by the planet's semi-major axis and eccentricity relative to the star (ap, ep) and the moon's semi-major axis and eccentricity relative to the planet (am, em). We wish to investigate how the extra dimensions of this manifold produces features unique to exomoon habitability, and how other factors, such as the moon's orbital direction, affect the resulting climate.
We investigate this problem in two parts: in the first, we simulate four simple test cases to study in detail the effect of orbital dynamics on the resulting climate. In the second, we carry out a large parameter study of simulations to investigate the four dimensional habitable manifold.
In section 2 we describe the LEBM and its implementation; in section 3 we display results from our investigations; in section 4 we discuss the results, and in section 5 we summarise the work.

Simulation Setup
In all simulations, we fix the star mass at M * = 1M⊙, the mass of the host planet at Mp = 1MJup, and the mass of the moon at 1M⊕. This system has been demonstrated to be dynamically stable on timescales comparable to the Solar System lifetime (Barnes & OBrien 2002).
The planet's orbit around the star is specified by its semi-major axis ap and its eccentricity ep, and the moon's orbit around the planet is specified by its semi-major axis am, eccentricity em and inclination relative to the planet's orbital plane im (see Figure 1). The orbital longitudes of the planet and moon are defined such that φp = φm = 0 corresponds to the x-axis. Note that as the moon mass is much less than the planet mass, we assume am is equal to the distance of the moon from the planet, as this is approximately equal to the distance from the moon to the barycentre of the planetmoon system.
Comparing this to the more traditional setup for exoplanet LEBMs, there are two principal differences: (i) The epicyclic motion of the moon relative to the star (as described by the distance between the star and the moon, r * m), and (ii) The potential for stellar eclipses by the planet.
A third difference is the potential for an extra energy source in the form of tidal heating, which we approximate using the expressions in Scharf (2006). While our calculations are of a lower dimension than those of Heller & Barnes (2013)By allowing for diffusion of heat across latitudes, we also differ from Heller & Barnes (2013). We describe the equations that constitute the LEBM in the following section.

A One Dimensional Latitudinal Energy Balance Model with Tidal Heating
In this paper, we adapt the one dimensional LEBM applied to planets to function for moons with an atmosphere of similar composition to the Earth. LEBMs solve the following diffusion equation: where T = T (x, t) is the temperature at time t, x = sin λ, and λ is the latitude (between −90 • and 90 • ). This equation is evolved with the boundary condition dT dx = 0 at the poles. The (1 − x 2 ) term is a geometric factor, arising from solving the diffusion equation in spherical geometry.
C is the effective heat capacity of the atmosphere, D is a diffusion coefficient that determines the efficiency of heat redistribution across latitudes, S is the insolation received from the star, I is the atmospheric infrared cooling and A is the moon's albedo. In the above equation, C, S, I and A are functions of x (either explicitly, as S is, or implicitly through T (x, t)).
D is defined such that an Earth-like planet at 1 AU around a star of 1M⊙, with diurnal period of 1 day will reproduce the average temperature profile measured on Earth, see e.g. Spiegel, Menou & Scharf (2008). Planets that rotate rapidly experience inhibited latitudinal heat transport, due to Coriolis effects (see Farrell 1990). We follow Spiegel, Menou & Scharf (2008) by scaling D according to: where ω d is the rotational angular velocity of the planet, and ω d,⊕ is the rotational angular velocity of the Earth. This is a simplified expression, which does not describe the full effects of rotation, but allows for rapid computation without severely compromising the model's accuracy. The diffusion equation is solved using a simple explicit forward time, centre space finite difference algorithm (as was done in Forgan 2012). A global timestep was adopted, with constraint . ( 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. This is clearly not applicable for certain orbital parameters, as tidal locking will play a significant role at low values of am (Joshi 1997).
The atmospheric heat capacity depends on what fraction of the planet's surface is ocean, focean, what fraction is land f land = 1.0 − focean, and what fraction of the ocean is frozen fice: The heat capacities of land, ocean and ice covered areas are The infrared cooling function is where the optical depth of the atmosphere The albedo function is This produces a rapid shift from low albedo to high albedo as the temperature drops below the freezing point of water, producing highly reflective ice sheets. It is this transition that makes the outer habitable zone extremely sensitive to changes to various orbital and planetary parameters, and makes LEBMs an important tool in studying short-term temporal evolution of planetary climates. The 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 where q0 is the bolometric flux received from the star at a distance of 1 AU, and Z is the zenith angle: We have assumed a simple main sequence scaling for the luminosity. δ is the solar declination, and h is the solar hour angle. The solar declination is calculated from the obliquity δ0 as: where φ * 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: 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 The radian half day length is calculated as In the interest of computational expediency, we make a simple approximation for tidal heating, by firstly assuming the tidal heating per unit area is (Peale, Cassen & Reynolds 1980;Scharf 2006): where Γ is the moon's elastic rigidity (which we assume to be uniform throughout the body), Rm is the moon's radius, ρm is the moon's density, Mp is the planet mass, am and em 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, Γ = 10 11 dyne cm −2 (appropriate for silicate rock) and ρm = 5 g cm −3 .
To calculate how this heating is distributed latitudinally, we assume the heating rate is at its maximum at the subplanetary point.
We reuse equation (16), substituting q0 forėT , δ for the appropriate δ of the moon relative to the planet (which in this case is equal to zero as we assume δ0, i.e the relative obliquity between the moon and the planet, is zero), and fixing H = π (i.e. tidal heating occurs throughout the diurnal period of the moon). This is very much an approximation -indeed, the multi-dimensional nature of tidal heating prohibits latitudinal models from improving much on approximations such as this -habitability studies typically assume tidal heating is uniformly distributed throughout the body (see e.g. Heller & Barnes 2013). We calculate habitability indices in the same manner as Spiegel, Menou & Scharf (2008). The habitability function ξ is 2 : We then average this over latitude to calculate the fraction of habitable surface at any timestep: Each simulation is allowed to evolve until it reaches a steady or quasi-steady state, and the final ten years of climate data are used to produce a time-averaged value of ξ(t),ξ, and the sample standard deviation, σ ξ . We use these two parameters to classify each simulations as follows: (i) Habitable Moons -these moons possess a time-averagedξ > 0.1, and σ ξ < 0.1ξ, i.e. the fluctuation in habitable surface is less than 10% of the mean.
(ii) Hot Moons -these moons have temperatures above 373 K across all seasons, and are therefore conventionally uninhabitable (ξ < 0.1).
(iii) Snowball Moons -these moons have undergone a snowball transition to a state where the entire moon is frozen, and are therefore conventionally uninhabitable (ξ < 0.1).

RESULTS
Exomoon systems present a parameter space of high dimension. We therefore take two approaches to exploring subsets of this space. The first approach fixes the planet's orbit, and considers four different sets of moon orbital parameters as test cases.
The second is a larger survey of parameter space, systematically varying ap, ep, am, em (with the other orbital parameters held fixed). This allows us to map out a four-dimensional manifold to represent the habitable zone, a subset of the true, higher dimensional habitable zone manifold (which will depend on extra parameters such as the obliquity of the moon and its chemical composition, which we do not vary).

Four Test Cases
We fix the orbit of the planet with ap = 0.9 AU with ep = 0, and consider four simple test cases to probe the effects of epicylic motion and eclipses on exomoon habitability.
(i) A prograde circular orbit (i.e. where the orbital angular momenta of the planet and moon are aligned), with am = 0.023AU = 481Rp (approximately 0.3 Hill Radii), with im = 90 • (ii) As 1., but with a retrograde orbit (where the orbital angular momenta are anti-aligned, and im = 270 • ).
These parameter choices are not entirely physically motivated -for example, it is unlikely that the orbit described in case 3 would remain stable on long time scales (Weidner & Horne 2010), not to mention that synchronous rotation would become important to the moon's climate. We choose these cases to illustrate more clearly the dynamical effects of changing the orbital parameters of the moon. We select ap = 0.9 AU to guarantee that all four test cases are habitable (see following section). All four moons are entirely habitable (ξ = 1). This is easily established by studying the mean, minimum and maximum surface temperatures ( Figure 2). The minimum temperature does not decrease beneath 290 K throughout the orbit of the moon and the planet. The differences between prograde and retrograde moons are immediately obvious (top left and top right panel respectively). While the mean temperatures appear to be similar, the maximum and minimum temperatures fluctuate much more strongly in the retrograde case. Case 3 is very similar to case 1, and case 4 also has a similar mean temperature (with a more pronounced oscillation in the maximum and minimum temperatures).
Why should we see such a large change by changing something as simple as the direction of orbit, and yet see little change by reducing the moon's orbital distance by a factor of 5? The answer lies in considering the mean temperature more carefully (Figure 3). The mean temperature fluctuates by < 0.1 K, in all cases, but it is the nature of the fluctuations that are telling. The retrograde orbit is consistently hotter than the prograde orbit, and the frequency of temperature fluctuations is also higher. Case 4, where the exomoon orbits perpendicular to the plane, exhibits a superposition of frequencies, and case 3, with low am, exhibits essentially no variation at all.
We can find an explanation for these phenomena if we consider the orbital dynamics. Let us consider first the difference between the prograde and retrograde cases. The planet's orbit is circular, and as such has no apastron or periastron. The moon's orbit is also circular relative to the planet, and hence has no apojove or perijove. But, relative to the star it undergoes epicyclic motion, and hence the moon has well-defined periastron and apastron. Figure 4 shows the location of the moon's periastron over the course of the planet's orbit.
The longitude of moon periastron, φper, precesses with a period equal to the planet's orbital period, P pl . We define the x-axis as corresponding to φ = 0, and therefore initially, φper = π radians. So, φper(t) = π +φ pl t mod 2π.
If the planet were static (φ pl = 0) then the longitude of periastron would not precess, and the time taken for the moon to move from apastron to periastron would be whereφm = 2π/Pm = nm, and Pm is the moon's orbital period. However,φ pl = 0, and therefore the longitude of periastron moves during this interval. If the moon's orbit is prograde, the angular distance between apastron and periastron increases, and the timescale becomes τP = π +φ pl τṖ φm .
Rearranging the above equations gives and τR = π φm +φ pl respectively. Therefore, the ratio of one to the other is The orbital periods of both bodies are related by (Appendix B of Kipping 2009a): where χ is am normalised by the moon's Hill Radius. This gives the more digestible For the object masses and semi-major axes used here, this corresponds to a ratio of 0.758. Therefore, the forcing timescale of retrograde moons is faster than prograde exomoons, and hence we might expect to see more variations in temperature for these orbits. Comparing the period of oscillations for cases 1 and 2, we find the ratio is indeed approximately 0.758. We should also note the amplitude of the oscillations, with the prograde displaying an amplitude of around 0.14 K, compared to the retrograde amplitude of 0.1K. This is again sensible considering the orbital dynamics: τP > τR, the prograde moon will spend more time at larger and smaller distances from the star, allowing the mean temperature to increase and decrease appropriately. The retrograde moon will encounter apastron and periastron more frequently, but spend less time there. This will not allow the mean temperature to rise as much, but will facilitate some latitudes to increase and decrease their temperature more readily (which gives an explanation for the high oscillations seen in the maximum and minimum temperatures in Figure 2).
Case 3 has a significantly lower orbital period (approximately 3.6 days), compared to the 41 day orbit of the other moons. As the diurnal and orbital periods are extremely close, and the periastron and apastron distances are extremely similar, the effect of epicyclic motion is extremely weak, and the climate becomes difficult to distinguish from that of an Earth-like planet orbiting with parameters (ap, ep).
The counter-intuitive behaviour of Case 4 becomes clear when the alignment of the two orbital angular momentum vectors, Lp and Lm are considered. Lp is aligned with the z-axis, and Lm is aligned with the x-axis in the negative direction (see Figure 5). At t = 0, the vector product Lp × Lm = L ′ is aligned with the negative y direction. As there are no external torques acting on the system, L ′ maintains this alignment throughout the orbit, and as such the moon's periastron and apastron distances become a function of the orbital longitude of the planet.
At t = 0, the periastron and apastron distances are equal, and the distance between the star and the moon r * m will not change during its orbit. After 0.25 planetary orbital periods have elapsed, the periastron and apastron distances are now equal to that in the other 3 cases. Similarly, at t = 0.5 planetary orbital periods the periastron and apastron distances are the same, and at t = 0.75 planetary orbital periods they reach their maximum values once more. This provides the amplitude modulation seen in the bottom right panel of Figure 3. Note the frequency of the oscillation is the same as that of case 1, as they are both prograde orbits.

Parameter study
The previous section has shown us that while the global properties of the exomoon are relatively insensitive to the orbital architecture, the detailed behaviour of the exomoon's climate depends on the orbital direction and period. We shall now look at how the global properties depend on the parameters more easily determinable by exoplanet and exomoon searches: the semi major axes of both objects relative to their host, and their eccentricities. To do this, we carry out over 3500 separate climate simulations, and classify each one as hot, cold, habitable or transient according to the classification system outlined previously.

A Control -The Circumstellar Planetary Habitable Zone
To provide grounds for comparison, we ran a separate series of 440 simulations with an Earth-like planet orbiting a Sun-like star to map out the traditional circumstellar HZ according to our classification system. The results of this are displayed in Figure 6. Green points represent simulations where the planet is classified as habitable; purple points where the planet is classified as a transient; and the red and blue points indicate planets that are too hot and too cold respectively.
The model does not contain a carbonate-silicate cycle, unlike e.g. , which modified the atmospheric CO2 pressure in accordance with temperature dependent weathering rates. Lower temperatures produce lower weathering rates, and as a result the atmospheric CO2 (produced by volcanic outgassing) cannot be sequestered. Therefore, cooler planets can be expected to have higher atmospheric concentrations of CO2, boosting the greenhouse effect and moving the outer edge of the HZ to higher semi-major axes than we see in Figure 6 (cf Kasting, Whitmire & Reynolds 1993).
The extension of the HZ at low ep to semi-major axes as low as 0.7 AU is a reflection of our (fairly lenient) criterion for habitability -namely, that 10% of the planet's surface remains habitable over a ten year period, with a standard deviation less than 10% of the mean habitable area. As ep is increased, σ ξ increases quickly, producing planets which are habitable on a seasonal basis only. If we are to compare to traditional habitability studies, then we should infer that the inner edge of the HZ is at approximately 0.8 AU for ep = 0. Equally, many of the transient classifications in this study would normally have been considered to be uninhabitable (as many of these simulations undergo seasonal periods when the habitable surface fraction is close to zero, but this is not sufficient to label them as hot or cold planets). We should bear this in mind as we consider the habitable zone for exomoons in the following sections.  Figure 6. The habitable zone for an Earth-like planet around a Sun-like star, as calculated from a LEBM using the classification system outlined in this paper, for comparison with the exomoon habitable zones displayed in later figures. Each point represents a simulation run with these parameters, and the colour of the point indicates its outcome. Red points produce hot planets with no habitable surface; blue points produce cold planets with no habitable surface; green points represent warm planets with at least ten percent of the surface habitable and low seasonal fluctuations; purple points represent warm moons with high seasonal fluctuations.  Figure 7. The exomoon habitable zone, as a function of host planet semi-major axis ap and eccentricity ep, for moons with zero orbital eccentricity. Each point represents a simulation run with these parameters, and the colour of the point indicates its outcome. Red points produce hot moons with no habitable surface; blue points produce cold moons with no habitable surface; green points represent warm moons with at least ten percent of the surface habitable, and low seasonal fluctuations; purple points represent warm moons with high seasonal fluctuations. In each plot, the moon's orbital semi-major axis relative to the host planet, am is fixed at a value between 0.1 and 0.3 Hill Radii.

Zero Moon Orbital Eccentricity
We first consider the habitable zone manifold in the case where em = 0, (and consequently the tidal heating rate is zero) and we consider several different values of am. This allows us to construct ap−ep slices of the HZ manifold to compare with the typical planet HZ shown in the previous section. Figure 7 shows four ap − ep slices, for four different values of am. In each simulation, we fix am relative to the Hill Radius at the given value of ap, so the absolute value of am will increase with increasing ap.
Comparing the top left panel of Figure 7 with Figure 6, we can immediately see that the habitable zone exists in general at lower semi-major axis. This is especially true at high ep: the inner edge of the HZ at e = 0.5 is now at 0.84 AU, as compared to 0.89 AU in the planetary case. As am is increased, the high ep component of the inner HZ boundary shifts to higher and higher ap. We therefore surmise that the effect of frequent eclipses allows the moon to remain cooler, and damps the fluctuations in temperature experienced as the planet eclipses the star.
In all cases, the low ap, low ep habitable tail observed in the planetary case disappears in the exomoon case. While frequent eclipses can cool the climate and damp oscillations, the habitable surface fraction of these moons is generally low (typically close to the minimum 0.1). Therefore, a relatively small increase in the value of σ ξ due to eclipses is large enough to ensure σ ξ > 0.1ξ, and as a result be classified as transient moons rather than habitable moons.
Comparing the outer edge of the HZ in the planet and moon cases, we can see it takes on a significantly different shape. While the planet outer HZ moves steadily outward as ep is increased, the moon outer HZ changes direction as ep is increased. Initially habitable at 0.99 AU for ep = 0, the outer boundary moves inward as ep increases, finding its minimum value of ∼ 0.95 AU at ep = 0.2, and then moves outward again. This appears to be due to a combination of planet and moon motion relative to the star. The moon's epicyclic motion at ep = 0 changes its distance from the star by only a small amount, and hence the climate remains stable enough to be classified as habitable. As ep is increased, the epicyclic motion of the moon is modulated by the motion of the planet as it moves between its own periastron and apastron. This extra modulation is sufficient to trigger a positive feedback in albedo, and a subsequent snowball effect. As ep increases to large values, the strong heating the moon experiences as the planet passes through its periastron is sufficient to keep it warm throughout the rest of the orbit, and the outer HZ boundary remains at ap ∼ 0.99 AU. As the eclipse duration increases with increasing am, the outer HZ boundary is pushed inwards for larger ep as am is increased.

High Moon Orbital Eccentricity
We now repeat the above analysis, but increase the moon's eccentricity to ep = 0.1. This is quite a large value to assume for moon eccentricity -without eccentricity pumping, it is unlikely that moons will possess eccentricities as large as this after tidal evolution with their host planet -we present these results as an extreme case. Intermediate eccentricities will produce results somewhere between the results of this section and the results of the previous section.
Indeed, Figure 8 illustrates how extreme this case is by the effect on the climate of the simulations. As tidal heating increases as am is decreased, it is clear that tidal heating dominates the proper-ties of highly eccentric moons for am < 0.2 Hill Radii (top panels), and all simulations are classified as hot.
As am is increased, a habitable zone appears. This zone is much narrower than seen previously, and is restricted to large values of ap, as tidal heating allows for warmer solutions when insolation is weak. In the case where am = 0.21 Hill Radii (bottom left panel of Figure 8), the HZ exists at ap > 0.95 AU, and potentially extends beyond 1.1 AU at high ep. The HZ is typically around 0.1 AU in width, although at low ep it is as narrow as 0.03 AU.
When am is increased to 0.26 Hill Radii, the HZ begins to more closely resemble the HZ seen for em = 0 moons. While still around 0.1 AU in width, it now displays the knee in the outer HZ boundary seen previously, and it is centred on similar values of ap as before. At this value of am, the tidal heating has become small compared to the traditional heating from insolation. This being the case, tidal heating still makes its presence felt at low ap, where simulations that would have been habitable for ep = 0 are now seasonally habitable due to the extra forcing from tidal heating producing stronger climate oscillations.

Limitations of the Model
As with all studies using 1D LEBMs, the nature of the assumptions made limit the applicability of the results obtained in this analysis. The models are insensitive to the longer term climate variations produced by oceanic circulation, as well as the compensatory effects of the carbonate-silicate cycle which can halt the albedodriven snowball effect. These effects, with timescales ranging between a few thousand to a million years  are much longer than the orbital periods of both the moon and the host planet, and as such are of less consequence (in the short term) than the dynamical forcing due to eclipses and tidal heating.
Assuming a value for the diffusion constant, D limits the application of LEBMs to exomoons in a different sense. D is calibrated such that a fiducial Earth-Sun climate model produces the correct latitudinal temperature distribution as obtained from real data (Spiegel, Menou & Scharf 2008). Therefore, D contains a great deal of hidden information about the atmospheric and thermodynamic properties of the Earth. Modifying these properties, e.g. the atmospheric pressure, has significant consequences for the habitability of the object Vladilo et al. 2013) and it is less clear how D ought to be modified for frozen moons such as Europa, or moons with fundamentally different chemistry such as Titan, without more data on their own latitudinal temperature distributions. Even with this data, a full understanding of how tidal heating affects the temperature balance at any latitude is essential for the correct calibration of D. Equally, a more accurate implementation of tidal heating in the LEBM is vital for future studies of exomoons.
Tidal forces are expected to affect land and ocean according to their elastic rigidities. While we have assumed a fixed ocean fraction of 0.7 throughout this analysis, the LEBM does not contain information as to how these oceans are distributed across the surface. Instead, the ocean and land components are assumed to be uniformly mixed over all latitudes. Future calculations of tidal heating will require the LEBM to be more specific about the distribution of landmasses. Ocean fraction may therefore become an important parameter in future studies of the exomoon HZ manifold (which should also investigate the effects of changing other  Figure 8. The exomoon habitable zone, as a function of host planet semi-major axis ap and eccentricity ep, for moons with eccentricity em = 0.1. Each point represents a simulation run with these parameters, and the colour of the point indicates its outcome. Red points produce hot moons with no habitable surface; blue points produce cold moons with no habitable surface; green points represent warm moons with at least ten percent of the surface habitable, and low seasonal fluctuations; purple points represent warm moons with high seasonal fluctuations. In each plot, the moon's orbital semi-major axis relative to the host planet, am is fixed at a value between 0.1 and 0.3 Hill Radii. parameters such as obliquity). Besides its effect on tidal heating, the ocean fraction also affects the thermal relaxation timescale of the moon -reducing the ocean fraction leaves the moon more susceptible to larger climate oscillations (Spiegel, Menou & Scharf 2008;Abe et al. 2011;Forgan 2012).
Tied to this problem of limited surface data is the limited dimension of the model itself -ignoring longitudinal information requires us to assume the moon is rapidly rotating relative to the insolation source. This will remain true for exomoons which are tidally locked to the planet -as they will still rotate with respect to the primary insolation source, the star -but the more complicated nature of the lunar day compared to a planetary day may mean that the rapid rotation assumption is satisfied more weakly at some points in the orbit compared to others.
We have also ignored one important heat source -illumination of the moon by the planet. Heller & Barnes (2013) show calculations for illumination from a tidally locked planet, incorporating the planet's thermal radiation and reflected starlight, with reflected starlight typically dominating unless the planet's albedo is very low. Adding the reflected starlight alone can increase the outer HZ boundary by a few percent -if this was incorporated into the LEBM, it could act to prevent the "knee" in the outer boundary at moderate values of ep.
Finally, we assume orbital parameters for all objects involved, and we do not evolve these parameters under the presence of any gravitational field. Incorporating the LEBM component into a more involved computation of the tidal evolution of the star-planet-moon system (e.g. Laskar & Robutel 1993;Sasaki, Barnes & O'Brien 2012) would allow an investigation of a unique set of Milankovitchesque cycles produced by the moon's orbital evolution (see e.g. Spiegel et al. 2010), with an accurate computation of the moon's obliquity evolution and tidal heating as a bonus.

Implications for Observations
Our simulations suggest that determining the orbital eccentricity of any exomoons detected in or near the conventional HZ will be crucial in assessing whether the moon itself is also habitable, much more so in fact than determining the host planet's eccentricity. Exomoons with even a moderate amount of eccentricity of ∼0.1 lead to a more limited habitable-zone with respect to the other system parameters as a result of tidal heating dominating. One may consider that the HZ extends outwards in aP for such scenarios but longer period planets are considerably rarer in the known transiting exoplanet catalogue, due to both geometric bias and lower composite signal-to-noise.
The orbital eccentricity of a moon will affect the position and duration of moon transits in the light curve (Kipping 2011) and is, in principle, an observable quantity. It is therefore quite reasonable that constraints on exomoon eccentricity could be determined in the case of a confirmed detection. It is also worth noting that exomoons are expected to rapidly tidally circularize even if starting with a large initial eccentricity due to a capture, for example Porter & Grundy (2011);Williams (2013).
To maintain a large exomoon eccentricity over Gyr would likely require forcing due to resonant massive satellites, but the occurence of multiple captures of Earth-mass moons around a single planet is surely less than that of single captures of Earth-mass moons.
We also find that the sense of orbital motion affects the climate of moons and in principle this may also be determined observationally (Kipping 2009b). However, in general such a determination is highly challenging with Kipping et al. (2012) finding even optimstic synthetic data is marginal for uniquely determining this to high-confidence using Kepler. Of course, once a confirmed signal has been found, follow-up with larger instrumentation such as JWST would greatly increase our ability to make this assesment.
Finally, in this work we have established an LEBM capable of modelling hundreds of realizations of exomoon configurations. Due to the likely low signal-to-noise of any future exomoon detection one should expect fuzzy and complex parameter posteriors, as demonstrated in synthetic tests by Kipping et al. (2012). We therefore suggest coupling an LEBM, such as that presented here, directly with representative realizations drawn from the joint posterior distribution of the relevant parameters. This would enable us to compute a marginalised distribution for the habitability of an exomoon (ξ).

CONCLUSIONS
In this work, we have used a one dimensional latitudinal energy balance model (LEBM) to investigate the evolution of climate on an Earth-like exomoon. We have incorporated the effects of stellar insolation, tidal heating, eclipses of the star by the planet, atmospheric cooling and heat diffusion, and investigated how the dynamical circumstances of exomoons affect the resulting climate.
We simulated four test cases to study the exomoon climates in detail, finding that exomoons in orbits retrograde to their host planet's orbit experience stronger climate oscillations than their prograde equivalents, and are around 0.1K warmer. In the case of moons with sufficiently short orbital periods, we find that the finite thermal relaxation time of the atmosphere allows it to act as a buffer against the otherwise strong rapid temperature oscillations produced by frequent eclipses.
We then went on to carry out an extensive parameter study of the four dimensional space constituted by the planet's semi-major axis and eccentrictity (ap, ep) and the moon's semi-major axis and eccentricity relative to the planet (am, em). Comparing to the habitable zone of an Earth-like planet orbiting the same star, we find that if the exomoon's orbit is circular, the exomoon habitable manifold tends to exist at lower ap thanks to the cooling effect of eclipses, and the inner edge of the habitable zone can exist at lower ap for higher ep, again as eclipses keep the moon sufficiently cool. If the exomoon's orbit is highly eccentric, the heating budget is dominated by tidal effects, and the habitable manifold is much narrower in ap, and exists at greater distance from the star.
The simplicity of LEBMs, plus their (albeit limited) ability to be augmented by new heat sources and sinks, makes them a useful tool for studying Earth-like exomoon climates. As they are computationally inexpensive, it is straightforward to run them multiple times, whether it is to map out their behaviour in a set parameter space as done here, or to assess the habitability of detected exomoons, by using realisations of the joint posterior distribution derived from observations. We are keen to hone and utilise this modelling technique for both of these uses.