Habitability of Exoplanetary Systems

TL;DR

This dissertation investigates habitability in exoplanetary systems by examining planets and moons, focusing on habitable zones, tidal heating effects, and methods to detect ice-covered satellites through photometric observations.

Contribution

It provides a comprehensive analysis of habitability beyond planets, including moons, and explores new detection techniques for icy satellites in exoplanetary systems.

Findings

Moons can be habitable due to tidal heating even if their planets are not.

Tidal heating can sustain subsurface oceans on icy moons, supporting potential life.

Photometric methods can identify ice-covered satellites by their high reflectance during occultation.

Abstract

The aim of my dissertation is to investigate habitability in extra-Solar Systems. Most of the time, only planets are considered as possible places where extraterrestrial life can emerge and evolve, however, their moons could be inhabited, too. I present a comprehensive study, which considers habitability not only on planets, but on satellites, as well. My research focuses on three closely related topics. The first one is the circumstellar habitable zone, which is usually used as a first proxy for determining the habitability of a planet around the host star. The word habitability is used in the sense that liquid water, which is essential for life as we know it, may be present on the planetary surface. Whether the planet is habitable or not, its moon might have a suitable surface temperature for holding water reservoirs, providing that tidal heating is in action. Tidal heating is generated inside the satellite and its source is the strong gravitational force of the nearby planet. The second topic of my research explores tidal heating and the habitability of extra-solar moons with and without stellar radiation and other related energy sources. Life is possible to form even on icy planetary bodies, inside tidally heated subsurface oceans. The third topic probes the possibility of identifying an ice-covered satellite from photometric observations. A strong indication of surface ice is the high reflectance of the body, which may be measured when the moon disappears behind the host star, so its reflected light is blocked out by the star.

Tables

Table 1. Table 3.1: Plateau duration for different planetary masses at 3 AU distance from a Sun-like star.

		plateau duration	plateau duration	plateau duration
Name	$a$ [km]	with a 1 $M_J$ mass	with a 5 $M_J$ mass	with a 10 $M_J$ mass
		planet [h]	planet [h]	planet [h]
Enceladus	237,378	3.83	2.00	1.42
Europa	670,900	10.83	9.52	6.74

Table 2. Table 4.1: Constants of the fitted equations for calculating the stellar temperature ( T 𝑇 T ) and the luminosity ( L 𝐿 L ).

	fit1			fit2
	$a$	$b$	$c$	$p$	$q$
$T$	$1379\pm 232$	$6219\pm 494$	$-1915\pm 262$	$5639\pm 19$	$0.438\pm 0.016$
$L$	$0.465\pm 0.102$	$-2.345\pm 0.409$	$2.950\pm 0.330$	$1.050\pm 0.045$	$4.266\pm 0.284$

Table 3. Table 4.2: Reduced chi-square ( χ 2 superscript 𝜒 2 \chi^{2} ) values for the five different models. The first column describes the name of the used method; the second, third, fourth and fifth columns show the χ 2 superscript 𝜒 2 \chi^{2} values for the stellar temperature, luminosity, inner and outer boundaries of the HZ, respectively.

Method	${\chi }_T^2$	${\chi }_L^2$	${\chi }_{\mathrm{inner}}^2$	${\chi }_{\mathrm{outer}}^2$
Zaninetti	$2.11$	$8.13$	$9.33$	$9.03$
Razzaque SPL	$2.26$	$7.58$	$10.94$	$10.40$
Razzaque BPL	$6.48$	$7.35$	$10.53$	$9.93$
fit1	$1.53$	$8.01$	$9.20$	$8.85$
fit2	$1.66$	$7.14$	$9.20$	$8.84$

Table 4. Table 4.3: Number of used equations for calculating the effective temperature ( T 𝑇 T ), the luminosity ( L 𝐿 L ) and the boundaries of the HZ.

Method	$T$	$L$	HZ
Zaninetti	3.3	3.4	3.1a, 3.1b
Razzaque SPL	3.5	3.6	3.1a, 3.1b
Razzaque BPL	3.7	3.8	3.1a, 3.1b
fit1	4.1a	4.1b	3.1a, 3.1b
fit2	4.1b	4.1b	3.1a, 3.1b

Table 5. Table 5.1: Q ​ μ 𝑄 𝜇 Q\mu values for Solar System-like, rocky moons. The radii and densities are those of the corresponding Solar System bodies. The reference for these values is Murray & Dermott ( 1999 , Appendix A) .

label	$R_m\left[\mathrm{km}\right]$	$\rho \left[\mathrm{kg}/m^3\right]$	$e$	$P\left[\mathrm{days}\right]$	$T_{\mathrm{surf}}\left[K\right]$	${\log }_{10}\left(Q\mu \right)$ [Pa]
Earth-like	6378	5515	0.01	1	281	14.0
				2	213	13.0
				3	180	12.4
			0.1	1	378	15.5
				2	291	14.4
				3	249	13.8
Mars-like	3394	3933	0.01	1	256	12.5
				2	194	11.5
				3	163	10.9
			0.1	1	342	14.0
				2	263	13.0
				3	225	12.4
Moon-like	1738	3342	0.01	1	232	11.1
				2	176	10.1
				3	n.a.	n.a.
			0.1	1	313	12.6
				2	240	11.5
				3	205	10.9
Io-like	1821	3532	0.01	1	235	11.2
				2	177	10.2
				3	n.a.	n.a.
			0.1	1	315	12.7
				2	242	11.7
				3	207	11.1

Table 6. Table 6.1: Calculated values of the MO depth for different bodies around a Sun-like star. A g subscript 𝐴 g A_{\mathrm{g}} : geometric albedo, R m subscript 𝑅 m R_{\mathrm{m}} : radius of the exomoon. The distance from the central star is 3 AU in each case. The geometric albedo values are obtained from Pang et al. ( 1981 ); Verbiscer et al. ( 2007 , 2013 ) .

Name	$A_g$	$R_m$ [km]	$y_3-y_4$ [ppm]
Enceladus	1.38	252	$4.3\cdot {10}^{-7}$
Europa	1.02	1561	$1.2\cdot {10}^{-5}$
Io	0.74	1821	$1.2\cdot {10}^{-5}$
Earth	0.37	6371	$7.5\cdot {10}^{-5}$
icy Earth	1.38	6371	$2.8\cdot {10}^{-4}$

Table 7. Table 6.2: Estimated albedo values at different photometric bands, based on the spectrum of Enceladus (water ice surface, Verbiscer et al. ( 2006 , 2007 ) ) and Eris (methane ice surface, Alvarez-Candal et al. ( 2011 ) ).

Surface	V	R	I	J	H	K
Water ice	1.375	n.a.	1.00	0.84	0.52	0.43
Methane ice	0.97	1.02	0.98	0.65	0.45	0.42

Table 8. Table 6.3: Instrument parameters for calculating the photon noise level. A 𝐴 A : aperture, Q 𝑄 Q : quantum efficiency, T campaign subscript 𝑇 campaign T_{\mathrm{campaign}} : observation period per target field for survey instruments. In most cases rough estimations were used. Reference numbers shown in the last column are explained in the footnote. In the cases of TESS and CHEOPS the quantum efficiencies are estimated based on the Q ​ ( λ ) 𝑄 𝜆 Q(\lambda) curves provided by the manufacturers. Asterisk indicates that instead of the quantum efficiency, the system throughput is presented.

Instrument	$A$ [m]	Filter	${\lambda }_f$ [nm]	${\lambda }_{\min }$ [nm]	${\lambda }_{\max }$ [nm]	$Q$ [%]	$T_{\mathrm{campaign}}$	Ref.
TESS	0.1	–	800	600	1000	80	80 days	1
CoRoT	0.27	–	650	400	900	63	150 days	2
CHEOPS	0.3	–	750	400	1100	67	–	3, 4
PLATO 2.0	0.68	–	750	500	1000	50*	2 years	5
Kepler	0.95	–	660	420	900	50*	4 years	6
HST ACS WFC	2.4	F555W	541	458	621	35*	–	7
HST WFC3 IR	2.4	F125W	1250	1100	1400	52*	–	8
JWST NIRCam	6.5	F115W	1200	1000	1400	40*	–	9
JWST NIRCam	6.5	F140M	1400	1300	1500	45*	–	9
E-ELT MICADO	39	J	1258	1181	1335	60*	–	10

Table 9. Table 6.4: Summary of various factors that influence the successful albedo estimation for an exomoon. The most influencing ones are listed in the left column, the less important and indifferent ones are listed in the right column.

Substantial influence	moderate influence	minimal influence
$\bullet$ stellar distance	$\bullet$ eccentricity of the exomoon	$\bullet$ direct/retrograde orbital
$\bullet$ exomoon size	$\bullet$ inclination of the exomoon	direction
$\bullet$ albedo	$\bullet$ hemispherical asymmetry	$\bullet$ reflected light from the
$\bullet$ star–planet size ratio	of the moon	planet
$\bullet$ planetary distance	$\bullet$ mass of the planet

Equations

$${\dot{E}}_{\mathrm{tidal}}=\frac{21}{2}\frac{k_2}{Q}\frac{G{M}_{\mathrm{p}}^2{R}_{\mathrm{m}}^{5}n{e}^2}{a^6},$$

$$k_2=\frac{3/2}{1+\frac{19\mu }{2\rho g{R}_{\mathrm{m}}}},$$

$$r_{\mathrm{inner}}=\sqrt{\frac{L}{S_{\mathrm{inner}}}},$$

$$r_{\mathrm{outer}}=\sqrt{\frac{L}{S_{\mathrm{outer}}}},$$

$$S_{\mathrm{inner}}=1.296-2.139\cdot 10^{-4}T+4.19\cdot 10^{-8}{T}^2,$$

$$S_{\mathrm{outer}}=0.234-1.319\cdot 10^{-5}T+6.19\cdot 10^{-9}{T}^2$$

$$T=5760\cdot M^{0.48}\mathrm{if}0.3M_{\odot }<M<18.5M_{\odot },$$

$$L=1.15\cdot M^{3.43}.$$

$$T=5700\cdot M^{0.5}\mathrm{if}M\le 10M_{\odot },$$

$$L=M^{3.6},$$

$$T=5700\cdot M^{0.8}\mathrm{if}0.1M_{\odot }\le M<2M_{\odot },$$

$$L=M^{4.8}\mathrm{if}M<2M_{\odot }.$$

$${\chi }^2=\frac{1}{N}\sum _{i=1}^N\frac{{\left[y_i-f(x_i)\right]}^2}{{\sigma }_i^2},$$

$${\dot{E}}_{\mathrm{tidal}}=-\frac{21}{2}\mathrm{Im}(k_2)\frac{R_{\mathrm{m}}^5{n}^5{e}^2}{G}.$$

$$-\mathrm{Im}(k_2)=\frac{57\eta \omega }{4\rho g{R}_{\mathrm{m}}\left[1+{\left(1+\frac{19\mu }{2\rho g{R}_{\mathrm{m}}}\right)}^2\frac{{\eta }^2{\omega }^2}{{\mu }^2}\right]},$$

$$\eta ={\eta }_0\exp \left(\frac{E}{\mathcal{R}T}\right),$$

$$\mu =10^{\left(\frac{{\mu }_1}{T}+{\mu }_2\right)}\mathrm{Pa},$$

$$\eta ={\eta }_0\exp \left(\frac{E}{\mathcal{R}T}\right)\exp \left(-B\varphi \right),$$

$$\eta =10^{-7}\exp \left(\frac{40000\mathrm{K}}{T}\right){\left(1.35\varphi -0.35\right)}^{-5/2}\mathrm{Pa}\mathrm{s}.$$

$$\eta =10^{-7}\exp \left(\frac{40000\mathrm{K}}{T}\right)\mathrm{Pa}\mathrm{s}.$$

$$q_{\mathrm{BL}}=k_{\mathrm{therm}}\frac{T_{\mathrm{mantle}}-T_{\mathrm{surf}}}{\delta (T)},$$

$$\delta (T)=\frac{d}{2a_2}{\left(\frac{Ra}{R{a}_{\mathrm{c}}}\right)}^{-1/4}$$

$$Ra=\frac{\alpha g\rho d^4{q}_{\mathrm{BL}}}{\eta (T)\kappa k_{\mathrm{therm}}}.$$

$$T_{\mathrm{surf}}={\left(\frac{{\dot{E}}_{\mathrm{tidal}}}{4\pi R_{\mathrm{m}}^2\sigma }\right)}^{1/4},$$

$$C\frac{\partial T}{\partial t}-\frac{\partial }{\partial x}\left(D\left(1-x^2\right)\frac{\partial T}{\partial x}\right)=\left(S+S_{\mathrm{p}}\right)\left[1-A(T)\right]+\zeta -I(T),$$

$$\zeta =\frac{21}{38}\frac{{\rho }^2{R}^5{e}^2}{Q\mu }{\left(\frac{G{M}_{\mathrm{p}}}{a^3}\right)}^{5/2}.$$

$$\delta t<\frac{{\left(\Delta x\right)}^{2}C}{2D\left(1-x^2\right)}.$$

$$A(T)=0.525-0.245\tanh \left[\frac{T-268\mathrm{K}}{5\mathrm{K}}\right].$$

$$S=\frac{q_0}{\pi }\left(H\sin \lambda \sin \delta +\cos \lambda \cos \delta \sin H\right),$$

$$\cos H=-\tan \lambda \tan \delta ,$$

$$q_0=1.36\cdot 10^6{\left(\frac{M_{\star }}{M_{\odot }}\right)}^4\mathrm{erg}{\mathrm{s}}^{-1}\mathrm{c}{\mathrm{m}}^{-2},$$

Full text

1. 1 Preface

2. 2 Introduction

3. 2.1 Habitability of exoplanets

4. 2.1.1 Habitable zones

5. 2.1.2 Liquid water habitable zone

6. 2.2 Habitability of exomoons

7. 2.2.1 Circumplanetary habitability

8. 2.2.2 The importance of tidal heating

9. 2.2.3 Icy moons

10. 2.2.4 Solar System analogues

11. 3 Methods

12. 3.1 Habitable zone

13. 3.2 Tidal heating

14. 3.2.1 Viscoelastic model

15. 3.2.2 Internal structure and convection

16. 3.2.3 Climate models

17. 3.3 Albedo estimation

18. 3.3.1 Calculation of ‘plateau duration’ and ‘MO duration’

19. 3.3.2 Calculating the photon noise level of different instruments

20. 3.3.3 Orbital distances around M dwarfs

21. 4 Calculating the Liquid Water Habitable Zone

22. 4.1 Effective temperature – mass relation

23. 4.2 Luminosity – mass relation

24. 4.3 Habitable zone

25. 4.4 Conclusion

26. 5 Viscoelastic Models of Tidally Heated Exomoons

27. 5.1 Viscoelastic model

28. 5.2 Comparison to the fixed Q model

29. 5.2.1 Method

30. 5.2.2 Surface temperature

31. 5.2.3 Occurrence rate of ‘habitable’ moons

32. 5.2.4 The value of Qµ

33. 5.2.5 Scaling the Galilean satellite system

34. 5.3 Tidal heating with climate models

35. 5.4 Conclusions

36. 6 The possibility for albedo estimation of exomoons

37. 6.1 Solar-type stars

38. 6.2 Moons in systems of M dwarfs

39. 6.2.1 Plausibility of large icy moons

40. 6.2.2 Different wavelengths

41. 6.2.3 Next generation telescopes and the occultation of exomoons

42. 6.3 Discussion

43. 6.3.1 Accuracy

44. 6.3.2 Considering different parameters

45. 6.3.3 Importance of ice

46. 6.4 Conclusions

47. 7 Summary

48. 8 Short description

49. 9 Kivonat

50. 10 Acknowledgements

1.   Preface

The aim of my dissertation is to investigate habitability in extra-Solar Systems. Most of the time, only planets are considered as possible places where extraterrestrial life can emerge and evolve, however, their moons could be inhabited, too. I present a comprehensive study, which considers habitability not only on planets, but on satellites, as well.

My research focuses on three closely related topics. The first one is the circumstellar habitable zone, which is usually used as a first proxy for determining the habitability of a planet around the host star. The word habitability is used in the sense that liquid water, which is essential for life as we know it, may be present on the planetary surface.

Whether the planet is habitable or not, its moon might have a suitable surface temperature for holding water reservoirs, providing that tidal heating is in action. Tidal heating is generated inside the satellite and its source is the strong gravitational force of the nearby planet. The second topic of my research explores tidal heating and the habitability of extra-solar moons with and without stellar radiation and other related energy sources.

Life is possible to form even on icy planetary bodies, inside tidally heated subsurface oceans. The third topic probes the possibility of identifying an ice-covered satellite from photometric observations. A strong indication of surface ice is the high reflectance of the body, which may be measured when the moon disappears behind the host star, so its reflected light is blocked out by the star.

Before submerging in these topics, I give an overview in Chapter 2 about the current state of research regarding the habitability of planets and their moons outside the Solar System. I present the used methods for my investigations for all three major topics in Chapter 3. My work regarding the habitable zone is described in Chapter 4, the tidal heating calculations are shown in Chapter 5 and the proposed method for selecting icy moons by observations is presented in Chapter 6. Finally, I summarise the results in Chapter 7.

2.   Introduction

2.1 Habitability of exoplanets

More than three thousand extra-solar planets (exoplanets) are known as of today (source: exoplanet.eu, July 2016) and even more planet candidates wait for confirmation. The first planets discovered were hot Jupiters that are of super-Jupiter size and orbit the central star at a very close distance. We do not expect that these planets could harbour life as we know it, but with current technology we are capable of detecting Earth-size planets, as well.

2.1.1 Habitable zones

The usual method for searching life is to find those planets which are located in their star’s habitable zone, hence they may be suitable for supporting life. Depending on the inspected conditions, several habitable zones can be distinguished. These are the Galactic Habitable Zone (for details see Lineweaver et al., 2004), the Liquid Water Habitable Zone (Kasting et al., 1993), the Ultraviolet Habitable Zone (Buccino et al., 2006; Dobos et al., 2010), the Tidal Habitable Zone (Barnes et al., 2009), the Photosynthesis-Sustaining Habitable Zone (Franck et al., 2000; von Bloh et al., 2009) and the Life Supporting Zone (Neubauer et al., 2012).

The emergence of life needs long time, given that on Earth it took about a billion years after the end of the heavy bombardment (Jones, 2004, Chapter 3). The orbit of the planet needs to be stable during this long time interval to stay in the central star’s habitable zone.

The so-called Continuous Habitable Zones describe whether a planet is located continuously inside the habitable zone for a long time. As the star ages, its radiation increases, hence the habitable zone moves outward from the star (see for example Kasting et al., 1993; Underwood et al., 2003; Buccino et al., 2006; Guo et al., 2009). This means that a planet which was initially inside the habitable zone gets stronger stellar radiation with time, which can be harmful for living organisms. On the other hand a farther planet, that was too far from the star and did not receive sufficient energy may become habitable.

In the followings, the Liquid Water Habitable Zone will be discussed in details.

2.1.2 Liquid water habitable zone

The most studied habitable zone is the liquid water habitable zone (HZ) which is usually described as a region around a star where a hypothetical Earth-like planet could support liquid water on its surface (Kasting et al., 1993). A schematic figure of the HZ in the Solar System can be seen in Fig. 2.1. In order to calculate the boundaries of the HZ, the stellar parameters are needed. Other environmental elements, like the presence of an atmosphere, its thickness, composition, atmospheric pressure, or cloud formation can influence the surface temperature, i.e. the state of the water, as well. Kasting et al. (1993) used a climate model based on carbonate-silicate cycle, atmospheric composition and stellar radiation. This model is widely used and is considered to be the best tool to estimate the boundaries of the liquid water habitable zone.

On the surface of the Earth, the water reservoirs wouldn’t boil at 100 ∘C, because if the atmosphere were at this temperature, than the atmospheric pressure would be 2 bars (Sullivan & Baross, 2007, Chapter 4). However, water molecules can move to the atmosphere approximately at 60 ∘C, enhancing the greenhouse effect, which leads to a warmer temperature. Because of the higher temperature even more water molecules evaporate, resulting in a positive feedback process. In the upper atmosphere dissociation occurs, meaning that ultraviolet photons break up the molecules to atoms, and finally hydrogen atoms escape to space, preventing new water molecules to form in the future. This process is called the runaway greenhouse effect (Kasting et al., 1993). Usually, this runaway greenhouse limit is used as the inner boundary of the liquid water habitable zone.

Water vapour is responsible for about the two-third of the total greenhouse effect on Earth. The other major greenhouse gas is carbon-dioxide, which is essential in the carbonate-silicate cycle, too (Sullivan & Baross, 2007, Chapter 4). The importance of this cycle is that the circulation of CO2 stabilizes the climate on Earth over a long timescale. Atmospheric carbon-dioxide dissolves in rainwater and forms carbonic acid (H2CO3, see Fig. 2.2). Silicate rocks are dissolved by this weak acid and the products of this silicate weathering (e.g. bicarbonate and dissolved silica) are transported by rivers to oceans. In the ocean these are built into shells in the form of calcium carbonate (CaCO3). When these organisms die, they fall to the bottom of the ocean, where the shells either redissolve, or become buried in sediments on the seafloor. At plate boundaries, when the oceanic plate slides beneath the continental plate, the carbonate sediments sink to high temperature and high pressure depths. In the process of volcanic eruption, CO2 is released into the atmosphere, which is called carbonate metamorphism. This is a full cycle of carbone-dioxide, which has a timescale of approximately 200 million years (Kasting et al., 1993; Sullivan & Baross, 2007).

The climate regulation is due to the negative feedback mechanism of carbonate-silicate cycle. If the temperature is higher, then precipitation rate will be higher, too, removing more CO2 from the atmosphere. Consequently, the greenhouse effect will lessen, resulting in a lower surface temperature. On the other hand, decreasing temperature leads to slower weathering rates, leading to the accumulation of CO2 in the atmosphere, which will strengthen the greenhouse effect. This process can even bring a planet back from a snowball state, providing that the temperature is not low enough to permit the formation of CO2 ice clouds, which would further cool the surface because of the release of latent heat, and the increased albedo (Caldeira & Kasting, 1992).

Since carbonate-silicate cycle can reverse the frozen state of a planet, it can push the outer boundary of the habitable zone relatively far out from the star. For this reason, the carbonate-silicate cycle is usually taken into account when calculating the boundaries of the HZ. Mostly, calculations are made for Earth-like bodies, since the climate models are very complex, and different planetary sizes, atmospheric compositions, cloud coverage, or continent-to-ocean ratios would fundamentally change the required models. The HZ in the Solar System is located between 0.97 and 1.67 AU distances from the Sun (these are the runaway greenhouse and the maximum greenhouse limits, respectively, Kopparapu et al., 2013a, b).

The HZ of hotter stars is located farther than for cooler stars, because their stellar radiation is stronger (see Fig. 2.3). Also, the HZ is wider for hotter stars, for example it is approximately 18 AU wide for a B spectral type main sequence star, but only 0.5 AU wide for an M class main sequence dwarf.

The habitable zone is a useful tool for selecting those planets where life could potentially emerge, however, even if a planet can be called ‘habitable’ based on liquid water, it does not necessary mean that it meets every requirements for life. A habitable planet, as described above, is one that can maintain liquid water reservoirs on the surface for a long period of time.

In general, F, G and K class main sequence stars are in the main focus in connection with habitability, since they are similar to our Sun, and spend sufficiently long time on the main sequence for the emergence of life. The habitable zone of the less massive M class stars is much closer to the star. M type red dwarfs are among the most active stars, with erosive flares, which can be harmful to life on close-in planets. A strong magnetosphere can protect the planet from cosmic rays, however, planets in the HZ of M dwarfs are probably tidally locked, which may result in a small magnetic moment (Grießmeier et al., 2005). Even for a magnetosphereless body, flares may not induce a direct hazard on the planetary surface, if the planet has a protective, oxygen-rich atmosphere (Segura et al., 2010). In addition, the abundance and long lifetime on the main sequence make M stars preferable targets in searching for life.

The climate model developed by Kasting et al. (1993) is complex, since it takes into account the carbonate-silicate cycle, the atmospheric composition of the planet and the stellar properties. For faster estimations of the location of the HZ, Underwood et al. (2003) fitted parabolic equations to the results of Kasting et al. (1993) and used these equations for calculating the boundaries of the liquid water habitable zone. Later one of the fitted equations was modified (Jones et al., 2006; Jones & Sleep, 2010). Selsis et al. (2007) proposed similar equations by extrapolating the results of Kasting et al. (1993). Catling & Kasting (2016) also made a parabolic fit to the same results. Kopparapu et al. (2013a) proposed a quartic function for calculating the stellar flux from the effective temperature. Note, that the correct coefficients can be found in Kopparapu et al. (2013b). To calculate the boundaries of the HZ with any of these equations, only the stellar temperature and luminosity are needed. Dobos et al. (2013) gave empirical formulae to calculate these parameters as functions of the stellar mass for F, G, K and M class main sequence stars.

It is also possible, that a planet is habitable, while it is outside the HZ. For example, the planet’s atmospheric composition and thickness can strongly influence its surface temperature, as well as large amount of clouds, atmospheric pressure, salinity of water, high eccentricity of the orbit or tidal heating. These all influence habitability, although these parameters are not considered in the traditionally called habitable zone.

2.2 Habitability of exomoons

Beside exoplanets, their moons (exomoons) can be habitable, too. No exomoons have been discovered yet, but these measurements are expected in the next decade. Bennett et al. (2014) present a candidate, which has been detected via the MOA-2011-BLG-262 microlensing event. The best-fit solution for the data implies the presence of an exoplanet hosting a sub-Earth mass moon. This measurement however needs confirmation, since an alternate solution is also presented. Nevertheless, this measurement indicates that the era of exomoon detection is about to begin.

The most favourable method for exomoon discoveries is photometry. An exoplanetary transit may reveal the presence of a moon in the light curve. Details of this method are thoroughly discussed in the literature (Simon et al., 2007; Kipping, 2009a, b; Simon et al., 2010; Kipping et al., 2012; Simon et al., 2012).

In addition, habitability of exomoons is under examination as well (see e.g. Kaltenegger, 2010; Heller & Barnes, 2013; Heller et al., 2014). Hinkel & Kane (2013) investigated the influence of eccentric planetary orbits on moons, and concluded that a moon with sufficient atmospheric heat redistribution may sustain suitable temperature for life on the surface even if the moon orbits a planet that moves temporarily outside of the circumstellar HZ at each orbital period.

Although moons in the Solar System are small compared to planets (the largest moon, Ganymede has a size of about 0.4 Earth-radius and a mass of 0.025 Earth-mass), it is widely believed that much larger moons can exist in other planetary systems, as well. One of the reasons is that much larger planets (hot Jupiters) exist, too, indicating that at the planet formation period they could have had much more massive protoplanetary discs, from which more massive moons might have been able to form. According to Crida & Charnoz (2012), the mass of a moon is limited by the mass of its host planet. Also, the mass of satellite systems in the Solar System is proportional to the mass of their host planet. Canup & Ward (2006) showed that this might be the case for extra-solar satellite systems as well, giving an upper limit for the mass ratio at around 10$-4$. Thus 10 Jupiter-mass planets may have 0.3 Earth-mass satellites. However, if the moon is not originated from the circumplanetary disk, but from collision, like in the case of the Earth’s Moon, then even larger satellites might exist. When studying the habitability of exomoons, it is very common to consider Earth-size moons, as they can have significant atmosphere, and they could be similar to Earth.

M dwarfs could have specific characteristics regarding the exomoon identification and analysis. Although most search projects focus on solar-like stars, M dwarfs as targets should also be considered, because more than 70% of the stars are of M spectral type. In theory it is possible that no large exomoons are present around M dwarfs, as their smaller stellar mass allow the formation of only smaller exoplanets and exomoons than in the Solar System. But based on the discovered exoplanets, 68 M dwarfs are known to have 96 exoplanets, from which the heaviest planet’s mass is 62 $M_{\mathrm{J}}$, 20 other exoplanets have at least 10 Jupiter-masses, and there are 13 other planets more massive than 1 $M_{\mathrm{J}}$ (exoplanet.eu, January 2016). The average mass of the exoplanets of the 68 M dwarfs is 4.8 $M_{\mathrm{J}}$.

2.2.1 Circumplanetary habitability

When considering the habitability of exomoons, it is usually assumed that the host planet orbits the central star inside the circumstellar habitable zone (the blue region in Fig. 2.4), and the moon’s orbit is between the so-called habitable edge and the Hill radius, or half of the Hill radius of the planet (this is the green circumplanetary region in Fig. 2.4). If the satellite is too close to the planet, then the reflected stellar light from the planet, the thermal radiation and tidal heating will be strong enough to make the moon’s surface too warm to be habitable. The habitable edge is defined by the runaway greenhouse limit, hence in orbits closer to the planet the moon will loose its water reservoirs, providing that the atmospheric conditions are similar to that of the Earth (Heller & Barnes, 2013). The Hill radius surrounds the gravitational sphere of influence for a body. This is the outermost distance where a moon can orbit a planet. However, according to Domingos et al. (2006) and Donnison (2010), the orbits of direct orbiting satellites can stay stable only inside the inner half of the Hill sphere.

Another outer limit can be defined by climate models. Forgan & Kipping (2013) developed a model for exomoons in order to investigate the energy balance of satellites. This climate model combined with eclipses and the ice-albedo feedback is capable of predicting the orbital parameters of a moon that will result in a snowball state, defining an outer limit for the circumplanetary habitable zone (Forgan & Yotov, 2014; Forgan & Dobos, 2016). The ice-albedo feedback is a positive loop where cooling results in larger ice coverage, which leads to a higher albedo that means that more stellar light will be reflected back, and less energy will be absorbed, so the cooling will be stronger.

A comprehensive model should include all the energy sources and energy sinks for the exomoon. The energy sources for an exomoon are: stellar insolation, atmospheric retention, radiogenic heat, planetary illumination, thermal heat from the planet and tidal heating. The possible energy sink sources are: infrared cooling, planetary eclipses and high albedo caused by oceans, ice or clouds. Eclipses can be frequent (in case of low orbital inclination relative to the planetary orbit) and prolonged on a satellite of a giant planet, hence they can reduce the incident flux by as much as 6% (Heller & Barnes, 2013). In calculations the albedo of exomoons are usually just estimated, or the albedo of Earth is used.

2.2.2 The importance of tidal heating

A moon can be habitable outside the circumplanetary, or even outside the circumstellar habitable zones, as well. Since tidal heating of the moon is caused by the presence of the planet, this energy source is independent from the star. Tidal heating is an inner source of energy that is caused by the gravitational forces of the nearby planet. Since the planet is much more massive than the moon, the arising forces deform the moon and cause friction inside the body that leads to heat dissipation. Hence, opposed to stellar radiation, tidal heating comes from the inside. For this reason, the surface temperature of a moon can be suitable for liquid water for several billion years even in the case when the planet-moon pair does not have a central star.

There are a few factors that favour the detection of tidally heated exomoons: 1. they can be far more luminous than their host planet, 2. their luminosities are independent of their separations from the star, and 3. they are visible around nearby (thus bright) stars (Peters & Turner, 2013).

The tidal heat rate of a moon is usually calculated by the following expression (e.g. Reynolds et al., 1987; Meyer & Wisdom, 2007):

[TABLE]

where $G$ is the gravitational constant, $M_{p}$ is the mass of the planet, $R_{\mathrm{m}}$, $n$, $e$ and $a$ are the radius, mean motion, eccentricity and semi-major axis of the moon, respectively. $Q$ is the tidal dissipation factor and $k_{2}$ is the second order Love number:

[TABLE]

where $\mu$ is the rigidity, $\rho$ is the density and $g$ is the surface gravity of the satellite. This calculation method is called the fixed $Q$ model, because $Q$, $\mu$ and $k_{2}$ are considered to be constants.

It can be seen from Eq. 2.1 that in case the moon orbits the planet in a circular orbit, then the tidal heat rate will be zero. Since tidal forces tend to circularize the orbit of satellites, tidal heating attenuates quickly, unless some other force continuously maintains the eccentricity. If other moons exist which orbit the same planet, then it is possible that their orbits will be in mean motion resonance, which means that the ratio of their orbital periods can be expressed with small integers. Due to this resonance, they will experience regular, periodic gravitational pull that excites their eccentricity. Therefore mean motion resonance between the orbits of two or more satellites maintains their eccentricities in long timescales, which continuously sustain tidal heating.

The fixed $Q$ model is broadly used in tidal calculations, but highly underestimates the tidal heat of the body (Ross & Schubert, 1988; Meyer & Wisdom, 2007). Moreover, both $Q$ and $\mu$ are very difficult to determine, and vary on a large scale for different bodies: from a few to hundreds for rocky planets, and tens or hundreds of thousands for giants (see e.g. Goldreich & Soter, 1966). In addition, these parameters should not be constants, since they strongly depend on the temperature (Fischer & Spohn, 1990; Moore, 2003; Henning et al., 2009; Shoji & Kurita, 2014). As a consequence, tidal heat flux has a temperature dependency, as well: it reaches a maximum at a critical temperature ($T_{\mathrm{c}}$) as can be seen in Fig. 2.5. Between the solidus and the liquidus temperatures ($T_{\mathrm{s}}$ and $T_{\mathrm{l}}$, respectively) the material partially melts. Above the breakdown temperature ($T_{\mathrm{b}}$) the mixture behaves as a suspension of particles. The dashed curve represents the convective heat loss of the body. Circles indicate equilibria, for example, the solid circle between $T_{\mathrm{c}}$ and $T_{\mathrm{b}}$ is a stable equilibrium point. If the temperature increases, convective cooling will be stronger than the heat flux, resulting in a cooler temperature. In the case of decreasing temperature, the tidal heat flux will be the stronger, hence the temperature increases, returning the system to the stable point. The stable equilibrium between the tidal heat and convection is not necessarily located between $T_{\mathrm{c}}$ and $T_{\mathrm{b}}$, in fact, there are cases, when the two curves do not have intersection at all (see Henning et al., 2009, Fig. 6). In these cases tidal heating is not strong enough to induce convection inside the body.

In contrast to the fixed $Q$ model, viscoelastic models take into account the temperature dependency of the body, hence are more realistic. The first application of a viscoelastic tidal heating model to exomoons was introduced by Dobos & Turner (2015). This model did not include any other energy sources than tidal heating. For this reason it can only be used in cases when stellar radiation and other energy sources are negligible.

As mentioned above, comprehensive models exist which take into account all of the energy sources affecting moons (see for example Heller & Barnes, 2013; Forgan & Yotov, 2014). These models however, use fixed $Q$ models for calculating the tidal heating of the moon. Forgan & Dobos (2016) merged a comprehensive climate model with viscoelastic tidal heating calculation and compared the results with fixed $Q$ calculations.

2.2.3 Icy moons

Among the satellites in the Solar System, Europa (moon of Jupiter) and Enceladus (moon of Saturn) are the best candidates for supporting subsurface liquid water and related potential habitability and hence these are the most actively studied satellites. Both moons have icy crusts, that are at least a few tens of kilometres thick. Underneath the ice layer, a huge ocean resides, which is believed to be global in the case of Europa, and probably on Enceladus, too (Kivelson et al., 2000; Thomas et al., 2016). Both bodies are tidally heated which maintains the liquid phase state of the oceans. At the boundary of the ocean and the silicate layers, where the temperature is warm enough, the environment can be suitable for the appearance of prebiotic processes thus potentially further development toward simple forms of life.

Ice is very reflective in the visual wavelengths and for this reason, the albedo of an icy celestial body is usually very high. The Bond albedo (which is often called as albedo for short) is defined as the ratio of the reflected and the incident flux at all phase angles. The phase angle is the angle subtended at the observer between the star and a celestial body. In other words, the Bond albedo is the reflected radiation combined in all directions. The albedo of a fully reflective body, that reflects all incident radiation is one, and the lowest possible value for fully absorbing bodies is zero.

The geometric albedo however, can exceed one in some cases. The geometric albedo is the ratio of brightness at zero phase angle compared to an idealized, flat, fully reflecting, diffusely scattering disk (which is called Lambertian reflectance). Zero phase angle means that the observer measures the brightness of the body from the direction of the stellar light, so the visible hemisphere is illuminated completely. In most cases it is impossible to do such measurements in practice, since the Earth or the detector blocks out some part of the light. Usually the geometric albedo is measured in a phase angle very close to zero. This difference in the angle can be very significant, because in some cases the reflected light increases drastically when the phase angle approaches zero. This is due to the so-called opposition effect. Fig. 2.6 shows that Solar System planetary bodies reflect much more light when the phase angle is low. The decrease for higher phase angles is partly caused by the decreased area that is lit, but note that the Moon and Mars have the steepest phase function curve at low phase angles in Fig. 2.6 which is caused by the opposition effect.

There are two explanations for the opposition effect, one is called the shadow-hiding opposition effect (SHOE) and the other is the coherent backscatter opposition effect (CBOE). The SHOE is caused by the disappearing shadow at zero phase angle which is particularly important in fine powders with grain sizes less than 20 µm, because the electrostatic and van der Waals forces exceed the gravitational force of the planetary body (Hapke, 2012). Consequently, the grain particles form towers and bridges which cast long shadows except at zero phase angle. The Moon is known to have such fine, porous powders on the surface that cause strong SHOE.

The CBOE is basically an interference phenomenon: an incident wave front on a particulate medium is scattered twice or more before exiting the medium, and from the same wave front light will travel in the same path but in the opposite direction, as well. At zero phase angle the phase difference of the emerging wavelets will be zero and the amplitudes add coherently as they interfere with each other, resulting in higher intensity up to a factor of two (Hapke, 2012; Verbiscer et al., 2013). Europa and Enceladus show significant CBOE (Verbiscer et al., 2013).

The high albedo of planetary bodies favours occultation measurements, because higher reflectivity means higher flux difference in the light curve. As can be seen in Fig. 2.7, the stellar flux can be measured alone when the planet disappears behind the star (occultation), and just before and after this occultation, the total flux is higher because of the reflected stellar light arriving from the planet. Conversely, if we can measure the flux difference in occultation, the albedo of the planet can be obtained. Rowe et al. (2008) estimated the geometric albedo of a planet using this method. Besides planets, this method can also be used to moons in theory. Dobos et al. (2016) investigated the possibilities of such measurements for exomoons with next generation telescopes.

2.2.4 Solar System analogues

When considering the habitability of exomoons, Solar System satellites can serve as useful analogues, since similar moons could exist in other planetary systems, as well. Tidal heating is present in several moons in the Solar System, and lots of ice satellites orbit the giant planets, hence it may be useful to study them in more details before moving forward to the investigation of extra-solar moons.

Io, one of the Galilean moons of Jupiter, is the most volcanically active satellite in the Solar System. Voyager-1 took the first photo of a volcanic eruption on another planetary body than the Earth (see Fig. 2.8). Today several volcanoes are known on Io. According to Spencer et al. (2000), the mean surface temperature of the satellite is approximately 90-95 K, but at volcanoes it can be around 200-400 K. The main source of this internal energy is tidal heating. Io, Europa and Ganymede are in a 1:2:4 mean motion resonance, respectively, which maintains their orbital eccentricities. This eccentricity excites continuous tidal heating inside the bodies.

The temperature map of Enceladus shows warmer areas, the so-called tiger stripes (see Fig. 2.9) at the south pole (Porco et al., 2006). The stripes coincide with the location of plumes that are geyser-like jets (see Fig. 2.10). The erupting water originates from the subsurface ocean and is ejected to such heights that it basically feeds the E ring of Saturn (Juhász et al., 2007). Kite & Rubin (2016) showed that the eruptions are sustained by tidal heating on million year timescales.

On Europa, tidal friction leads to heat dissipation mainly in the silicate layer that is the underlying layer of the ocean. Part of the heat is dissipated probably through hydrothermal vents on the seafloor that heats the ocean locally. The warmer water moves upward to the ice layer, melting it and creating the so-called chaos terrains. At these regions the ice sheet of Europa becomes approximately 2 km thin (Kereszturi & Keszthelyi, 2013). It is possible, that similarly to Enceladus, Europa has plume activity, too (Roth et al., 2014).

The satellites of the giant planets hold substantial mass of water ice that is often located on the surface (Brown & Cruikshank, 1997; Verbiscer et al., 1990). Although there could be differences between satellite systems, the icy surface, including water ice, is abundant on them (Stevenson, 1985; Sasaki et al., 2010). Analysing the optical properties of the Solar System satellites, it can be seen that the high albedo in general indicates water ice. Non-ice covered satellites usually show low geometric albedo values between 0.04 and 0.15, while ice covered satellites show higher values, although different ingredients (mostly silicate grains) could decrease the albedo if they are embedded in the ice, or the ice is old and strongly irradiated/solar gardened (Grundy & Stansberry, 2000; Sheppard et al., 2005). These satellites show moderate geometric albedo between 0.21 for Umbriel and 1.00 for Dione (Veverka et al., 1991; Verbiscer et al., 2013). Beside water ice, nitrogen ice is also a strong reflector, producing albedo around 0.76 on Triton (Hicks & Buratti, 2004). On this satellite the bright surface can partly be explained by its relatively young, 6-10 Myr old age. (Schenk & Zahnle, 2007). The most important analogues for high albedo in the Solar System could be Europa with 1.02 and Enceladus with 1.38 geometric albedo (Verbiscer et al., 2013). In the case of Europa the relatively young (about 50 million year old) surface accompanies with clean and bright ice, while on Enceladus freshly fallen ice crystals increase the albedo significantly (Pappalardo et al., 1999; Porco et al., 2006).

The opposition effect (CBOE) is prominent on the rings (see Fig. 2.11) and nearby satellites of Saturn which results in high geometric albedos that in some cases exceed unity. Verbiscer et al. (2007) showed that there is a correlation between the reflectance profile of the E ring and the geometric albedo of nearby satellites. The plumes of Enceladus give material not only to the E ring, but eventually, the particles fall to the surface of the satellites, too. The clean ice coats them, enhancing their reflectance in opposition. The affected satellites are Mimas, Tethys, Dione and Rhea, and they have 0.96, 1.23, 1.00 and 0.95 geometric albedos, respectively.

3.   Methods

The used calculation methods for the circumstellar and circumplanetary habitability, and also for the albedo estimation are presented in this section. The description of the habitable zone calculation (section 3.1) is based on the work of Dobos et al. (2013), the tidal heating and the comprehensive climate models (section 3.2) are based on Dobos & Turner (2015) and Forgan & Dobos (2016) and the method for albedo estimation (section 3.3) is based on the work of Dobos et al. (2016).

3.1 Habitable zone

To calculate the boundaries of the liquid water habitable zone as a function of the stellar mass, the values of the stellar luminosity and effective temperature are needed. Many of the equations for calculating these parameters from the stellar mass have different terms for different stellar masses. The terms for F, G, K and M class main sequence stars are chosen and the equations are arranged to similar forms to provide an easier comparison.

The distance of the inner and the outer boundaries from the central star is given by

[TABLE]

where $r$ is the distance from the star in astronomical units, $L$ is the luminosity of the star in solar units and $S$ is the stellar flux at the boundaries in units of the solar constant (Kasting et al., 1993).

The value of $S_{\mathrm{inner}}$ and $S_{\mathrm{outer}}$ can be calculated by the fitted equations of Jones et al. (2006):

[TABLE]

where $T$ is the effective temperature of the star. Similar equations were given by Catling & Kasting (2016). These equations are equivalent.

The stellar temperature and luminosity can be calculated for main sequence stars from their mass as given by Zaninetti (2008):

[TABLE]

[TABLE]

Here $M$ is the mass of the star in solar units and $M_{\odot}$ represents one solar mass.

Other possible solutions are presented by the simple power law (SPL):

[TABLE]

[TABLE]

or by the broken power law (BPL) of Razzaque et al. (2009):

[TABLE]

[TABLE]

With these formulae the location of the liquid water habitable zone can be determined easily for different stellar masses between 0.3 and 2 $M_{\odot}$ which is the mass interval covered by all equations. Any of the equations above can be used for these stellar masses.

For numeric comparison of the different calculation methods the reduced chi-square ($\chi^{2}$) values were calculated for all cases. The following equation was used:

[TABLE]

where $N$ is the number of data points, $y_{i}$ is the measured value, $f(x_{i})$ is the calculated value and $\sigma_{i}$ is the uncertainty of the measurement.

3.2 Tidal heating

3.2.1 Viscoelastic model

In viscoelastic tidal heating models the term $k_{2}/Q$ in Eq. 2.1 is replaced by the imaginary part of the complex Love number Im($k_{2}$), which describes structure and rheology in the satellite (Segatz et al., 1988):

[TABLE]

Note that in this expression the mass of the planet and the semi-major axis of the moon are eliminated by the mean motion ($n=\sqrt{GM_{p}/a^{3}}$).

Henning et al. (2009) give the value of Im($k_{2}$) for four different models (see Table 1 in their paper). In this work we use the Maxwell model:

[TABLE]

where $\eta$ is the viscosity, $\omega$ is the orbital frequency and $\mu$ is the shear modulus of the satellite.

The viscosity and the shear modulus of the body strongly depend on the temperature. Below the $T_{\mathrm{s}}$ the shear modulus is constant: $\mu=50\,\mathrm{GPa}$ and the viscosity follows an exponential function:

[TABLE]

where $\eta_{0}=1.6\cdot 10^{5}\mathrm{Pa\,s}$ (reference viscosity), $E$ is the activation energy, $\mathcal{R}$ is the universal gas constant and $T$ is the temperature of the material (Fischer & Spohn, 1990).

Between $T_{\mathrm{s}}$ and $T_{\mathrm{b}}$ the body starts to melt. The shear modulus changes by

[TABLE]

where $\mu_{1}=8.2\cdot 10^{4}\mathrm{K}$ and $\mu_{2}=-40.6$ (Fischer & Spohn, 1990). The viscosity can be expressed by

[TABLE]

where $\phi$ is the melt fraction which increases linearly with the temperature between $T_{\mathrm{s}}$ and $T_{\mathrm{l}}$ ($0\leq\phi\leq 1$) and $B$ is the melt fraction coefficient ($10\leq B\leq 40$) (Moore, 2003).

At $T_{\mathrm{b}}$ the grains disaggregate, leading to a sudden drop in both the shear modulus and the viscosity. Above this temperature the shear modulus is set to a constant value: $\mu=10^{-7}\mathrm{Pa}$. The viscosity follows the Roscoe-Einstein relationship so long as it reaches the liquidus temperature (where $\phi=1$) (Moore, 2003):

[TABLE]

Above $T_{\mathrm{l}}$ the shear modulus stays at $10^{-7}\mathrm{Pa}$, and the viscosity is described by (Moore, 2003)

[TABLE]

In our calculations rocky bodies are considered as satellites, and for this reason we follow the melting temperatures of Henning et al. (2009), namely: $T_{s}=1600$ K, $T_{l}=2000$ K. We assume that disaggregation occurs at 50% melt fraction, hence the breakdown temperature will be $T_{b}=1800$ K.

3.2.2 Internal structure and convection

The structure of the moon in the model is the following: the body consists of an inner, homogeneous part, which is convective, and an outer, conductive layer. If the tidal forces are weak, the induced temperature will be low, resulting in a smaller convective region and a deeper conductive layer. But in case of strong tidal forces, the temperature will be higher, hence the convective zone will be larger with a thinner conductive layer.

For calculating the convective heat loss, we use the iterative method described by Henning et al. (2009). The convective heat flux can be obtained from

[TABLE]

where $k_{\mathrm{therm}}$ is the thermal conductivity ($\sim 2\mathrm{W/mK}$), $T_{\mathrm{mantle}}$ and $T_{\mathrm{surf}}$ are the temperature in the mantle and on the surface, respectively, and $\delta(T)$ is the thickness of the conductive layer. We use $\delta(T)=30\,\mathrm{km}$ as a first approximation, and then for the iteration

[TABLE]

is used, where $d$ is the mantle thickness ($\sim 3000$ km), $a_{2}$ is the flow geometry constant ($\sim 1$), $Ra_{\mathrm{c}}$ is the critical Rayleigh number ($\sim 1100$) and $Ra$ is the Rayleigh number which can be expressed by

[TABLE]

Here $\alpha$ is the thermal expansivity ($\sim 10^{-4}$) and $\kappa$ is the thermal diffusivity: $\kappa=k_{\mathrm{therm}}/(\rho\,C_{\mathrm{p}})$ with $C_{\mathrm{p}}=1260\,\mathrm{J/(kg\,K)}$. For detailed description see the clear explanation of Henning et al. (2009).

Because of the viscosity of the material the thickness of the boundary layer and convection in the underlying zone changes strongly with temperature. The weaker temperature dependencies of density and thermal expansivity are neglected in the calculations. The iteration of the convective heat flux lasts until the difference of the last two values is higher than $10^{-10}\mathrm{W/m^{2}}$.

Calculations of tidal heat flux and convection are made for a fixed radius, density, eccentricity and orbital period of the moon. We assume that with time, the moon reaches the equilibrium state. Henning et al. (2009) showed that planets with significant tidal heating reach equilibrium with convection in a few million years. However, change in the eccentricity can shift, or destroy stable equilibria. After finding the stable equilibrium temperature, the tidal heat flux is calculated, from which the surface temperature can be obtained using the Stefan-Boltzmann law:

[TABLE]

where $\sigma$ is the Stefan-Boltzmann constant. This is the first time of using a viscoelastic model for obtaining the tidal heat induced surface temperature on exomoons.

3.2.3 Climate models

We adopt initial conditions essentially identical to those of Forgan & Kipping (2013) and Forgan & Yotov (2014). The stellar mass is $M_{\star}=1M_{\odot}$, the mass of the host planet is $M_{p}=1M_{\mathrm{Jup}}$, and the mass of the moon is $M_{\mathrm{s}}=1M_{\oplus}$. This system has been demonstrated to be dynamically stable on time-scales comparable to the Solar System lifetime (Barnes & O’Brien, 2002).

The planet’s orbit is given by its semi-major axis, $a_{\mathrm{p}}$, and eccentricity, $e_{\mathrm{p}}$ , and the moon’s orbit by $a$ and $e$, respectively. We assume that the planet resides at the barycentre of the moon-planet system, which is satisfactory given the relatively large planet-to-moon mass ratio. The inclination of the planet relative to the stellar equator is $i_{\mathrm{p}}=0$ (i.e. the planet revolves 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_{\mathrm{m}}$, is zero unless stated otherwise. The orbital longitudes of the planet and moon are defined such that they equal to zero corresponding 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.

The 1D latitude energy balance model employed in this work solves the following diffusion equation:

[TABLE]

where $T=T(x,t)$ is the temperature at time $t,x\equiv\sin\lambda$, and $\lambda$ is the latitude (between $-90^{\circ}$ and $90^{\circ}$). This equation is evolved with the boundary condition $\frac{\mathrm{d}T}{\mathrm{d}x}=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 atmospheric heat capacity, the diffusion coefficient $D$ controls latitudinal heat redistribution, $S$ and $S_{\mathrm{p}}$ are the stellar and planetary insolation, respectively, $\zeta$ is the surface heating generated by tides in the moon’s interior, $I$ is the atmospheric infrared cooling and $A$ is the albedo. $D$ and $I$ are calculated as described in Forgan & Dobos (2016), and the surface heat flow can be calculated as (Scharf, 2006)

[TABLE]

The diffusion equation is solved using a simple explicit forward time, centre space finite difference algorithm. A global time-step was adopted, with constraint

[TABLE]

This time-step constraint ensures that the first term on the left-hand side of equation 3.21 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 generally true, as the star is the principal insolation source, and the moon rotates relative to the star on time-scales of a few days.

The albedo function is

[TABLE]

This produces a rapid shift from low albedo ($\approx 0.3$) to high albedo ($\approx 0.75$) as the temperature drops below the freezing point of water, producing highly reflective ice sheets. This shift in albedo affects the potential for global energy balance, and for planets in circular orbits, two stable climate solutions arise, one ice-free, and one ice-covered. Spiegel et al. (2008) show that such a function is sufficient to reproduce the annual mean latitudinal temperature distribution on the Earth. Note that we do not consider clouds in this model, which could 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 albedo, which in truth is not likely to be the case.

The stellar insolation flux $S$ is a function of both season and latitude. The total diurnal insolation can be obtained from

[TABLE]

where $q_{0}$ is the bolometric flux received from the star at a distance of 1 AU, $H$ is the radian half-day length at a given $\lambda$ latitude and $\delta$ is the solar declination. The radian half-day length can be calculated from

[TABLE]

and the bolometric flux is given by

[TABLE]

We implement planetary illumination and eclipses of the moon in the same manner as Forgan & Yotov (2014) (see also Heller & Barnes, 2013). While Heller & Barnes (2013) allow for the planet to be in synchronous rotation and have a significant temperature difference between the dayside and nightside (expressed in the free parameter $\mathrm{d}T_{\mathrm{planet}}$), we assume the planetary orbit is not synchronous, and we fix $\mathrm{d}T_{\mathrm{planet}}=0$. The planetary albedo is fixed at 0.3.

The habitability function $\xi$ is

[TABLE]

We then average this over latitude to calculate the fraction of habitable surface at any time-step: and the bolometric flux is given by

[TABLE]

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 $\xi(t)$, $\bar{\xi}$, and the sample standard deviation, $\sigma_{\xi}$. We use these two parameters to classify each simulations as follows.

Habitable moons – these moons possess a time-averaged $\bar{\xi}>0.1$, and $\sigma_{\xi}<0.1\,\bar{\xi}$, i.e. the fluctuation in habitable surface is less than 10 percent of the mean. 2. 2.

Hot moons – these moons have average temperatures above 373 K during all seasons, and are therefore conventionally uninhabitable, and $\bar{\xi}<0.1$. 3. 3.

Snowball moons – these moons have undergone a snowball transition to a state where the entire moon is frozen, and are therefore conventionally not habitable. As with hot moons, we require $\bar{\xi}<0.1$ for the moon to be classified as a snowball, but given the nature of the snowball transition as it is modelled here, these worlds typically have $\bar{\xi}=0$. 4. 4.

Transient moons – these moons possess a time-averaged $\bar{\xi}>0.1$, and $\sigma_{\xi}>0.1\,\bar{\xi}$, i.e. the fluctuation in habitable surface is larger than 10 percent of the mean.

3.3 Albedo estimation

In this section the method of albedo calculations from occultation light curve is presented. Because of the wide range of possible conditions (distance of the moon from its host planet and host star, orbital elements, size, albedo, etc.) only general and simple cases are used here. With this approach we get insight into the most important parameters without analysing all possibilities. The potential significance of parameters not used in this model are discussed in section 6.3.2.

Three different configurations are possible just before and after any occultation:

the exomoon is behind the planet (both the planet and the moon are in the same line of sight), 2. 2.

the exomoon is in front of the planet (both the planet and the moon are in  the same line of sight), 3. 3.

the exomoon and the planet are separately ‘visible’ (although spatially not resolved).

It can be generally assumed that the configuration does not change much in most cases during the occultation. The reason is that the time duration of the occultation varies in a couple of hours time scale, while the orbital period of the exomoon is assumed to be longer (Cabrera & Schneider, 2007). A reverse scenario is unlikely, but also discussed by Simon et al. (2010) and Sato & Asada (2010).

From the three possible configurations, the third one, which favours for the identification of the exomoon, since in the other cases the moon’s presence does not change the shape of the light curve (although the depth of the flux minimum can be different). The third configuration is also the most probable of all, based on geometric considerations (Heller, 2014). Apart from a few special cases, the moon spends just a small fraction of its orbital time in front of, or behind the planet (Heller, 2012). In this (third) case two sub-configurations are possible: the planet occults before the moon (a1 in Fig. 3.1), or the moon occults before the planet (a2 in Fig. 3.1). These cases can be seen in Figure 3.1 along with schematic light curves. Note that unlike at transits, stellar limb darkening does not affect the light curve at occultations.

Simple formulae are used for the calculations, because this is a basic phenomenon of exomoons that has not been analysed before. Simplified cases are considered because of the large number of possibilities and also to see the big picture without being lost in the details.

For the albedo calculation the following formulae and configuration are considered. Fig. 3.2 shows a schematic light curve of an occultation, assuming that first the planet moves behind the star, followed by the moon (case a1 in Fig. 3.1). The light curve has its maximum at $y_{1}$ and the minimum of the occultation is at $y_{4}$. The latter corresponds to the stellar flux (the planet and the moon are behind the star), denoted by $F_{\star}$ which is normalized so that $F_{\star}=1$. Reflected stellar light from the planet and from the moon give contributions to the observed flux. Values of $y_{3}$ and $y_{2}$ refer to the measured flux before the occultation of the moon and after the occultation of the planet (the moon is still occulting), respectively. Measuring both $y_{1}-y_{2}$ and $y_{3}-y_{4}$ can be used to refine their value (which should be equal in theory), assuming that the apparent configuration of the system does not significantly change during the occultation. The accuracy of this value is very important, because it plays a key role in calculating the moon’s albedo.

For the difference of the measured fluxes:

[TABLE]

where $F_{\mathrm{m}}/F_{\star}$ is the relative flux contribution of the moon normalized to the stellar flux. The flux contribution of the planet is irrelevant for our study. We call the $y_{3}-y_{4}$ difference moon occultation (MO) depth. Adapting Eq. (2) from Rowe et al. (2008) to this case, the ratio of the fluxes can also be written as

[TABLE]

where $A_{\mathrm{g}}$ is the geometric albedo, $R_{\mathrm{m}}$ is the radius and $a_{\mathrm{p}}$ is the distance of the exomoon from its host star, which is approximated by the semi-major axis of the planet. The radius of the moon is assumed to be already known from previous transit observations. The geometric albedo is the ratio of brightness at zero phase angle compared to an idealized, flat, fully reflecting, diffusely scattering disk. The Bond albedo on the other hand, is the reflectivity at all phase angles. The value of the Bond albedo is restricted between 0 and 1. The geometric albedo can occasionally be higher than 1, especially for atmosphereless bodies with icy surface (e. g. Enceladus, Tethys) because of the strong backscattering and opposition effect (Verbiscer et al., 2007). Because of the coherent backscatter opposition effect (CBOE), the phase difference of the light waves in the exact backscatter direction will be zero and the amplitudes add coherently as they interfere with each other, resulting in higher intensity up to a factor of two (Hapke, 2012; Verbiscer et al., 2013).

From Eqs. 3.30 and 3.31 one can get:

[TABLE]

The geometric albedo of an exomoon can be calculated by obtaining the MO depth from observations during occultation, and by using Eq. 3.32. However, very precise measurements are required, hence it is reasonable to consider the MO depth as an upper limit (Rowe et al., 2006). Several observations are needed for a phase-folded light curve that show the presence of the exomoon because of the photometric orbital sampling effect (Heller, 2014).

Using Eq. 3.32 we calculated the MO depth values for different bodies, in order to investigate the required precision of future observations. These photometric measurements may lead to obtain the albedo of an exomoon with the use of the described method.

3.3.1 Calculation of ‘plateau duration’ and ‘MO duration’

In order to achieve the detection of small flux differences, sufficiently long detector integration time is required. An important limiting factor for detection is the duration while the exomoon is ’visible’ alone and the exoplanet is behind the star (the time interval $t_{3}-t_{2}$ in Fig. 3.2), or when only the exoplanet is visible without the moon. We call this phase ‘plateau duration’, because it causes a plateau in the light curve.

The plateau duration is calculated for ideal cases, when the exomoon is located at its maximum elongation from the exoplanet (where it has only radial velocity component). Circular orbits were assumed both for the planet and for the moon. This is a crude simplification, but in the case of captured or impact-ejected moons (like Triton around Neptune or the Moon around the Earth) tidal interactions can strongly reduce the satellite’s eccentricity, almost circularizing the orbit. With this simplified model the plateau duration can easily be calculated by using the speed of the planet-moon pair around the central star:

[TABLE]

where $v_{\mathrm{p}}$ is the orbital speed of the planet, $G$ is the gravitational constant, $M_{\star}$ and $M_{\mathrm{p}}$ are the masses of the star and the planet, respectively. Since circular orbits are considered, $a_{\mathrm{p}}$ indicates the planet’s orbital distance. In most cases the apparent location of the moon does not change significantly during the occultation. For this reason the plateau duration can be estimated from the orbital speed of the planet, if the orbital distance of the moon ($a$) is known (this is the distance that the moon needs to travel after the planet disappears behind the star). If the moon was not in its maximal elongation, then the distance (and hence the plateau duration) would be shorter.

However, the moon can be fast enough on its own orbit around the planet to influence the plateau duration. In case the moon passes 25% of its orbit (starting from its maximal elongation) we terminate the plateau duration, because the moon must be in front of or behind the planet. These cases were taken into account as well. The orbital speed of the moon ($v_{\mathrm{m}}$) was calculated by

[TABLE]

Both time intervals were calculated: $t_{\mathrm{a}}=a/v_{\mathrm{p}}$ and $t_{25}=0.25\left(2a\pi\right)/v_{\mathrm{m}}$, and the shorter was considered as the plateau duration ($T_{\mathrm{event}}$) in each case.

The plateau duration was calculated for different mass host planets. It can be shown that having 1, 5, or 10 $M_{\mathrm{J}}$ mass planets ($M_{\mathrm{J}}$ denotes the mass of Jupiter) does not change the order of magnitude of the results. Table 3.1 shows such cases for two moons. The stellar distance is 3 AU in each case, and a Sun-like star was used in the calculations. For the calculation, the radius of the planet was calculated from the mean density (1326 kg/m3) and original radius (71,492 km) of Jupiter in each case. For the orbital distance of the two moons the real semi-major axes were used.

The plateau duration depends on many factors, including the exomoon’s apparent distance from its host planet, its orbital velocity and the size of the host star. Larger stellar disk, lower velocity and larger planetary distance of the moon increase the length of the plateau duration. If the apparent planet-moon distance is larger than the diameter of the star, then the ‘plateau’ will not appear in the light curve, because the planet and the moon will occult separately, making two separate drops in the light curve (i.e. the planet’s occultation ends before the moon’s occultation starts or the moon’s occultation ends before the planet’s one starts). In other words the occultations of the moon and the planet will not overlap in time. In this case the plateau duration is replaced by the duration of the moon’s occultation. We call it ‘moon occultation duration’, or MO duration for short.

In these cases, when the size of the star limits the MO duration, the following simple calculation was used. The velocity of the moon is already known from previous calculations, and the distance that the moon travels during the MO duration approximately equals to the diameter of the star. Since it is assumed that the impact parameter of the moon and the inclination of both the planet and the moon are zeros, the related time can be calculated from $t_{\mathrm{MO}}=2R_{\star}/v_{\mathrm{p}}$. If $t_{\mathrm{MO}}<t_{\mathrm{a}}$, then $t_{\mathrm{MO}}$ will be considered as the MO duration ($T_{\mathrm{event}}$).

The phase space of different stellar and planetary distances of a hypothetical moon was mapped. The plateau duration for different stellar and planetary distances can be seen in Fig. 3.3. The masses of the star and the planet are 1 solar mass and 1$M_{\mathrm{J}}$, respectively, for all cases. Black and white colours indicate the longest (41 hours) and the shortest (6 minutes) durations, respectively. Below the dashed grey line the plateau duration is terminated, because the moon reaches 25% of its orbit. It mainly occurs for higher stellar distances, which produces lower orbital velocity around the star. Because of this lower velocity $t_{25}$ will be smaller than $t_{\mathrm{a}}$. Above the solid horizontal white line the plateau duration is constant for each stellar distance. The growing semi-major axis of the moon results in larger $t_{\mathrm{a}}$, but at about 1,393,000 km ($2R_{\star}$) $t_{\mathrm{a}}=t_{\mathrm{MO}}$, or in other words, above this planetary distance the planet and the moon occults separately (the MO duration takes the place of the plateau duration). A few examples from the Solar System are plotted for comparison.

3.3.2 Calculating the photon noise level of different instruments

It is assumed that the star emits perfect black-body radiation. This is a good assumption if we intend to estimate the photon noise level for a given instrument. The star is characterized by its radius $R_{\star}$ and $T_{\star}$, the effective temperature of the stellar surface. According to Planck’s law, the power density on unit surface and unit solid angle as the function of the wavelength $\lambda$ can be written as

[TABLE]

where

[TABLE]

which is proportional to the location of the peak of the spectral energy distribution. The constants are the following: $k_{\rm B}$ is the Boltzmann constant, $h$ is the Planck constant and $c$ is the speed of light. The integral over the stellar surface $(4\pi R_{\star}^{2})$ as well as on the irradiated hemisphere according to Lambert’s cosine law yields a total power output at a certain wavelength which can be computed as

[TABLE]

The number of photons emitted at a certain wavelength can then be computed by dividing the power $P_{\lambda}(\lambda)$ by the photon energy $hc/\lambda$. The total number of photons detected by a detector having a gross quantum efficiency of $Q(\lambda)$ at unit time is then

[TABLE]

Here $A$ is the area of the entrance aperture on the telescope (on which the detector is mounted) and $d$ is the distance to the star. In practical applications, $Q(\lambda)$ is the product of the transparencies of various filters, lens and reflective elements forming the telescope optics and the quantum efficiency of the detector itself. In order to simplify our calculations, $Q(\lambda)$ can be written as a characteristic function having a constant value of $Q$ within an interval of $[\lambda_{\rm min},\lambda_{\rm max}]$ (and it is zero outside this interval). Let us define $\Delta\lambda$ as the effective bandwidth of the filter as $\Delta\lambda=\lambda_{\rm max}-\lambda_{\rm min}$ and the central wavelength as $\lambda_{\rm f}=(\lambda_{\rm min}+\lambda_{\rm max})/2$.

Using this assumption above, the total number of photons received at unit time can be written as

[TABLE]

By introducing the variable $x=\lambda_{\star}/\lambda$ and applying the respective measure of ${\rm d}\lambda=(\lambda_{\star}/x^{2}){\rm d}x$, we can write

[TABLE]

If the bandwidth of the filter is narrow, i.e. $\Delta\lambda\ll\lambda_{\rm f}$, the integral above can further be simplified by assuming the expression $x^{2}(e^{x}-1)^{-1}$ to be constant in the interval of $[\lambda_{\star}/\lambda_{\rm max},\lambda_{\star}/\lambda_{\rm min}]$. In this case, $\Delta x=x_{\rm max}-x_{\rm min}$ can be approximated as $\Delta x\approx\lambda_{\star}\Delta\lambda/\lambda_{\rm f}^{2}$ and hence the photon flow is

[TABLE]

Once the photon flow, i.e. the number of photons received by the detector is known, the photon noise of the instrument can be computed for any exposure time ($T$) as

[TABLE]

If the filter bandwidth is comparable to the central wavelength, the integral in Eq. (3.40) can be computed numerically. Since the integrand is a well-behaved function, even a few steps using Simpson’s rule can be an efficient numerical method for this computation. We note that in the case of Kepler, TESS or PLATO 2.0, where CCDs are used without any filters, the bandwidth of the quantum efficiency curve is comparable to its central wavelength, so this numerical computation is more accurate. However, the error in the accuracy introduced by this simplification is comparable to the discretisation of the $Q(\lambda)$ function (i.e. it is in the range of $10-20\%$ in total, depending on the actual values of $\lambda_{\star}$, $\lambda_{\rm f}$ and $\Delta\lambda$).

In order to apply Eq. (3.41) in practice, it is worth converting $R_{\star}$, $d$, $\lambda_{\rm f}$ and $\lambda_{\star}$ in astronomically relevant units (such as solar radius, parsec and microns) instead of SI units. After these substitutions, we can write the approximation

[TABLE]

This approximation is accurate within a few percents w.r.t. Eq. (3.41). In the constant $10^{12}\,{\rm Hz}$ we included both the dimension conversions (solar radius, parsec, microns) as well as the other constants such as $2\pi$ and the speed of light.

If we intend to detect a signal $S$ over the (cumulative) time $T$, the corresponding signal-to-noise ratio is going to be

[TABLE]

The integration time is calculated as the product of the MO duration (or plateau duration) and the number of events (occultations): $T=T_{\mathrm{event}}\cdot N_{\mathrm{event}}$. For observatories, 30 events are assumed for the measurements, and for survey missions the number of events is calculated from the length of the mission campaign and the orbital period of the planet: $N_{\mathrm{event}}=T_{\mathrm{campaign}}/P_{p}$.

For all the calculations presented in this work, $d=50$ pc and $S/N=5$ is used, the latter corresponds to $5\,\sigma$ detection. In some cases (see Table 6.3 for specific cases), instead of the $Q$ quantum efficiency, the system throughput was used, which contains $Q$ and other optical attenuations as well. In order to make these different cases comparable, either the system throughputs, or 80 percent of the quantum efficiency ($0.8Q$) was used in the calculations.

3.3.3 Orbital distances around M dwarfs

Beside solar-like stars, I have applied the calculations to a range of small M dwarfs. The orbital distances of the planet and the moon are calculated as functions of the stellar mass. In the case of M dwarf host stars, the planet-moon pair is set to the snowline. Icy exomoons are expected to form beyond or around the so-called snowline, which is the distance where water ice condenses in the protoplanetary disk, supporting the formation of large mass planetesimals that finally evolve toward gravitationally collected gas giants (Ida & Lin, 2008; Qi et al., 2013).

The location of the snowline is calculated from the equilibrium temperature at the planet’s sub-stellar point ($T_{0}$) (Cowan & Agol, 2011):

[TABLE]

$T_{0}$ is assumed to be approximately 230 K at the ice condensation boundary. To obtain $R_{\star}$ and $T_{\star}$ from the stellar mass ($M_{\star}$), a parabolic equation was fitted to the effective temperature and radius values given by Kaltenegger & Traub (2009, Table 1). Only M4 – M9 stars were used for the fitting. The obtained equations are:

[TABLE]

[TABLE]

that can be used if $0.075<M_{\star}<0.2$.

The moon’s orbital distance is set to 0.45 Hill radius ($R_{\mathrm{H}}$), hence it also depends on the stellar mass. It was calculated from the following formula:

[TABLE]

We chose this distance in order to avoid orbital instability, since Domingos et al. (2006) and Donnison (2010) showed that beyond approximately half of the Hill radius, the orbit of a direct orbiting exomoon becomes unstable.

4.   Calculating the Liquid Water Habitable Zone

In this work the liquid water habitable zone is investigated as a function of the mass of the main sequence star. Section 3.1 gave an overview of different calculation methods for the boundaries of the habitable zone and for the effective stellar temperature and luminosity. These methods are compared to each other, discussed in this section, and new empirical formulae are proposed. The following description is based on the work of Dobos et al. (2013).

4.1 Effective temperature – mass relation

The boundaries of the HZ strongly depend on the stellar temperature, but different calculation methods give different results. In Fig. 4.1 the effective temperature can be seen as a function of the stellar mass. The red dots indicate measured parameters of main sequence stars with known planets (data from exoplanet.eu, 2012 November, Schneider et al., 2012). The curves labelled as Zaninetti, Razzaque SPL, and Razzaque BPL show the values calculated from Eqs. 3.3, 3.5, and 3.7, respectively. It can be clearly seen that the broken power law gives the worst estimation compared to the measured data points, because the calculated values are too low for lower mass stars (0.3 – 0.7 $M_{\odot}$) and too high for higher mass stars (1.1 – 1.4 $M_{\odot}$). The other two equations give much better results comparing to the data points, but still give a bit higher temperature for stars with masses of 1.2 – 1.4 $M_{\odot}$.

For this reason, two equations were fitted to the data points. The first equation (labelled as fit1 on Fig. 4.1) is a second order polynomial, and the other equation (fit2) has a similar expression than the ones from the previous models (Eqs. 3.3, 3.5, and 3.7). The fitted equations for fit1 and fit2 are

[TABLE]

respectively, where $f_{i}$ is $T_{i}(M)$ or $L_{i}(M)$, i.e. these equations are used for calculating both the stellar temperature and the luminosity. The fitted values can be seen in Table 4.1. These equations can be used in the 0.3 – 1.4 $M_{\odot}$ interval. For higher mass values $T_{1}(M)$, and for smaller masses $T_{2}(M)$ have large deviations.

For calculating $\chi^{2}$, $200$ K was considered as uncertainty in the temperature measurements. The two fitted curves reproduce the measured parameters more precisely, than the curves of the other models. For numeric comparison, the reduced chi-square values were calculated for all five cases and the results can be seen in the second column of Table 4.2. From the three previously known models the Zaninetti method, and from all cases the fit1 have the lowest $\chi^{2}$ values. The fit2 model has nearly as low $\chi^{2}$ value as the fit1 model. The fitted $p$ and $q$ parameters are in a good agreement with the parameters in Eqs. 3.3, 3.5, and 3.7, but have slightly lower values.

4.2 Luminosity – mass relation

Similar calculations were made for determining the stellar luminosity. All five methods were used as described in Eqs. 3.4, 3.6, 3.8, 4.1a and 4.1b. The results are shown in Fig. 4.2. Here the values represented by red dots are calculated from measured stellar parameters given by Jones et al. (2006). These parameters are the distance ($d$), the apparent visual brightness ($V$) and the bolometric correction ($BC$). The luminosity is calculated from the following expression (Jones et al., 2006)

[TABLE]

When calculating the stellar luminosity values, the following uncertainties were considered: $d\pm 1\mathrm{pc}$, $V\pm 0.1^{\mathrm{m}}$, $BC\pm 0.1^{\mathrm{m}}$. These uncertainties propagate into the luminosity and are shown in Fig. 4.2 as error bars. The same uncertentaintes were taken into account when fitting the luminosity functions with Eqs. 4.1a and 4.1b. If the uncertainties were not considered in the fitting, then the parameters would significantly differ from the ones that are shown in the second row of Table 4.1.

The $\chi^{2}$ values were calculated again for all cases, and they show that the formula of fit2 reproduces the stellar luminosity more precisely than the other formulae (see the $\chi^{2}_{L}$ column of Table 4.2). The value of $q$ in the $L_{2}(M)$ equation is between the corresponding parameters of the SPL and BPL cases (3.6 and 4.8, respectively), and like in the Zaninetti model, a multiplier ($p$) is also used, which improves the fitting. The fit1 formula, however, gives negative value for the stellar luminosity between approximately $0.38$ and $0.41\,M_{\odot}$, which obviously cannot be a good estimation for any star.

4.3 Habitable zone

For calculating the boundaries of the habitable zone, Eqs. 3.1a and 3.1b were used. As the effective temperature and the luminosity can be calculated from different equations, a comparison was made between the different methods. The used equations for each model can be seen in Table 4.3. The first column contains the names of the models, and the second, third and fourth columns show the number of the used equations for calculating the effective temperature, the stellar luminosity and the boundaries of the HZ, respectively. As the $L_{1}(M)$ equation gave negative values for the luminosity for some stellar masses, $L_{2}(M)$ (Eq. 4.1b) is used instead.

The results of the calculations are shown in Fig. 4.3 and in the last two columns of Table 4.2. In Fig. 4.3 the red dots indicate the location of the inner, and the light blue diamonds indicate the outer boundary of the HZ, calculated by Jones et al. (2006) using Eq. 3.2a and 3.2b. These values were calculated from measured stellar parameters, hence we consider them as good references for different calculation models. The error bars shown are the propagated uncertainties originating from the measurements of the distance, stellar temperature, visual brightness and bolometric correction. The curves in Fig. 4.3 indicate the different models as described in Table 4.3. (The names of the methods in Fig. 4.3 are the same as in the first column of the table.) With these calculation methods both the inner and the outer boundaries were calculated for the same stellar masses that are presented by the data points (red dots and light blue diamonds). The fit1 and fit2 curves shown in the top and bottom panels, respectively, are nearly identical.

The solid curves surrounding the fit1 and the fit2 curves indicate the uncertainties (variance) of the fitted models. They were estimated for each stellar mass from the following expression:

[TABLE]

where $\Delta L$ and $\Delta S$ are the estimated uncertainties of the luminosity and the stellar flux, respectively. The values of $\Delta L$ and $\Delta S$ were calculated with the same method (but with their own equations, as shown in Table 4.3), using the uncertainties of the parameters presented in Table 4.1.

It can be clearly seen in Fig. 4.3 that the uncertainty of the fitted models is larger for more massive stars, than for smaller stars. It is probably caused by the larger deviation and uncertainty of the luminosity observed for more massive stars. The uncertainty region at 1.4 stellar masses is wider in the case of the fit1 model (about 0.32 and 0.58 AU at the inner and the outer boundaries, respectively), than for the fit2 model (about 0.17 and 0.35 AU, respectively).

From the reduced chi-square values (see the last two columns of Table 4.2) it can be seen that the fit1 and fit2 models provide nearly the same precision and they have the lowest $\chi^{2}$ values among the different models. It means that these two models give the closest results to the calculation made by Jones et al. (2006) in the 0.3 – 1.4 $M_{\odot}$ interval.

4.4 Conclusion

Different calculation methods of the liquid water habitable zone were compared to each other. Three previously defined and two new fitted models were used for the calculations. The first fitted equation (fit1) is a parabolic function, and the other (fit2) has a similar form than the other models. The equations of the fit2 model are:

[TABLE]

These models were used for calculating the effective temperature and stellar luminosity. From the three already known models the one called Zaninetti gave the best results for the effective temperature, and for the luminosity the Razzaque BPL model gave the best reproduction to measured stellar parametes. The fit2 model gave even better results for both the stellar temperature and luminosity values. Using the fit2 method, better estimations can be obtained for the boundaries of the habitable zone of F, G, K and M class (0.3 – 1.4 $M_{\odot}$) main sequence stars.

Since the observed deviation of stellar parameters is large for similar stellar masses, it is possible, that even more precise models could be obtained, if another parameter than the stellar mass (for example the effective temperature) would be considered in the fitting, as well.

5.   Viscoelastic Models of Tidally Heated Exomoons

The idea of a circumplanetary, tidally-heated habitable zone has emerged and was investigated by several authors (e.g. Reynolds et al., 1987; Scharf, 2006; Heller & Barnes, 2013). For the first time, a viscoelastic model was applied for studying tidal heat in exomoons. This work aims to give a detailed study of the circumplanetary Tidal Temperate Zone (TTZ), and discusses the differences to other models. For a complex investigation of habitability on exomoons, the viscoelastic model was also coupled to a climate model. The work described in sections 5.1 and 5.2 are based on the paper of Dobos & Turner (2015), section 5.3 is based on the work of Forgan & Dobos (2016) and section 5.4 presents the conclusions of both works.

5.1 Viscoelastic model

The satellite’s surface temperature is calculated as described in section 3.2, for different orbital periods and radii, at a fixed density and eccentricity. Stellar radiation and other heat sources are not considered, and have been neglected. The orbital period and the radius of the moon vary between 2 and 20 days, and between 250 km and 6550 km, respectively. It is common to consider Earth-mass moons in extra-Solar Systems when speaking of habitability, however, their existence is not proven. In the Solar System the largest moon, Ganymede has only 0.025 Earth mass. But the mass of satellite systems is proportional to the mass of their host planet. Canup & Ward (2006) showed that this might be the case for extra-solar satellite systems as well, giving an upper limit for the mass ratio at around $10^{-4}$. This means, that 10 Jupiter-mass planets may have 0.3 Earth-mass satellites. Besides accretion, large moons can also form from collisions, as in the case of the Earth’s Moon, and in such cases the mass ratio of the moon and the planet can be even larger than $10^{-4}$. Other possibility is the capturing of terrestrial-sized bodies through a close planetary encounter, as described by Williams (2013). For these reasons, we also take Earth-like moons into account.

The results of the calculations can be seen in Fig. 5.1, where the density of the moon is that of Io, and its eccentricity is set to 0.1. Different colours indicate different surface temperatures. In the white region there is no stable equilibrium between tidal heat and convective cooling. In other words, tidal heat is not strong enough to induce convection. For comparison, a few Solar System moons are plotted that have similar densities to Io’s. Yellow contour curves denote 0 and $100\,^{\circ}$C. The green area between these curves indicates that water may be liquid on the surface of the moon (atmospheric considerations were not applied). We define this territory as the Tidal Temperate Zone.

Interestingly, the location of the TTZ strongly depends on the orbital period, and less on the radius of the moon. Low radii are less relevant, since smaller bodies are less capable of maintaining significant atmospheres.

The dependency on the eccentricity can be seen by comparing Figs. 5.1 and 5.2. In the case of the latter figure the moon’s eccentricity is 0.01. For most of the orbital period–radius pairs there is no solution (white area). Due to this drastic difference, Europa analogues get out of equilibrium for smaller eccentricities, and the TTZ becomes narrower and shifts to shorter orbital periods. Note, that radiogenic heat is not considered in the model, which could push back the moon into equilibrium state, and would result in higher surface temperature.

Similar calculations were made for the density of the Earth and Titan (top and bottom panels of Fig. 5.3, respectively). The densities do not have high influence on the tidally induced surface temperature, however, the TTZ slightly shifts to lower orbital parameters for higher densities. (The density of Earth, Io and Titan are $5515\,\mathrm{kg/m^{3}}$, $3528\,\mathrm{kg/m^{3}}$ and $1880\,\mathrm{kg/m^{3}}$, respectively.)

In the left panel of Fig. 5.3 an example Earth-like moon is plotted inside the TTZ. This hypothetical body has the same mean surface temperature (288 K), radius (6370 km) and density as the Earth, hence its orbital period is 2.06 days. In the right panel a few Solar System satellites are plotted that have similar densities to that of Titan.

The stellar flux for moons with ambient temperatures of $\sim$100 K (which is similar to the case of the Galilean and Saturnian moons in the Solar System) is about one percent of the tidal flux in the TTZ. For this reason, stellar insolation may be safely ignored if the planet-moon system orbits the star at a far distance, or if they are free-floating. In systems in which the stellar irradiation alone is sufficient to heat the surface to levels of order the melting temperature or higher, the models presented here would need to be replaced by more complex hybrid ones to take into account both sources of heat and their very different spatial distributions on and within the moon.

5.2 Comparison to the fixed Q model

5.2.1 Method

It is clear from the results, that the viscoelastic model does not give solution in the case of small tidal forces. In other words, the amount of heat that is produced by tidal interactions is insufficient to induce convective movements inside the body, and for this reason there is no equilibrium between them. In contrast, the fixed $Q$ model provides solution in these cases, as well. However, the viscoelastic model describes the tidal heating of the body more realistically than the fixed $Q$ model, due to the temperature dependency of the $Q$ and $\mu$ parameters. How are the results of the two models related to each other?

For comparing the results of the two kinds of models, we use the expression of Eq. 7 from Peters & Turner (2013) for the fixed $Q$ calculation:

[TABLE]

where $T_{\mathrm{surf}}$ is the surface temperature of the moon induced by tidal heating, $G$ is the gravitational constant, $\sigma$ is the Stefan-Boltzmann constant, $R_{\mathrm{m}}$ is the radius, $\rho$ is the density, $\mu$ is the elastic rigidity, and $Q$ is the dissipation function of the moon, $e$ is the eccentricity of the moon’s orbit, and $\beta$ is expressed with the semi-major axis ($a$) and the mass of the planet ($M_{\mathrm{p}}$):

[TABLE]

where $a_{R}$ is the Roche radius of the host planet. These equations can be used for calculating the surface temperature of the moon heated solely by tidal forces.

The satellite’s mean motion can be expressed from $\beta$ by

[TABLE]

which makes the comparison of the two models easier.

5.2.2 Surface temperature

For comparison of the fixed $Q$ and the viscoelastic model, see Figs. 5.4, 5.5 and 5.6, which show the surface temperature of a moon, calculated with both the viscoelastic (red solid curve) and the fixed $Q$ model (green dashed curve) as functions of the eccentricity, radius and orbital period of the satellite. For the density of the moon we used 5515 $\mathrm{kg/m^{3}}$, which is the density of the Earth, and for the fixed $Q$ model we used $Q=280$ and $\mu=12\cdot 10^{10}\,\mathrm{kg/(m\,s^{2})}$ in each case (Peters & Turner, 2013, Table 1). The radius, orbital period and eccentricity of the satellite are set to that of the Earth, Io and 0.03, respectively, except that one of these parameters is varied in each figure (horizontal axes). The horizontal light blue, dashed lines indicate 0 and 100 ∘C (making the boundaries of the TTZ), and the solid blue lines denote the minimum and maximum temperatures ($-20$ and $60\,^{\circ}$C) that are probable limits of habitability on an Earth-like body (Sullivan & Baross, 2007, Chapter 4). In salty solutions the lower limit for microbial activity is around $-20\,^{\circ}$C, and the upper limit for complex eukaryotic life is $60\,^{\circ}$C. The latter temperature is also about the runaway greenhouse limit for Earth. These limits are only used for Earth-like bodies ($\rho=\rho_{\mathrm{Earth}}$ and $R_{\mathrm{m}}\approx R_{\mathrm{Earth}}$). Vertical lines show a few examples from the Solar System for different eccentricities, radii and orbital periods.

It is noticeable that the red curve is less steep than the green one, and larger portion of it is located inside the TTZ, especially in Figs. 5.4 and 5.5. It shows that the viscoelastic model stabilizes the surface temperature comparing to the fixed $Q$ model. These are just a few examples indicating that the viscoelastic model is less sensitive to these parameters, and that there are huge differences in the results of the models. In the next section the volume of the TTZ is investigated more thoroughly.

5.2.3 Occurrence rate of ‘habitable’ moons

Habitability on extraterrestrial bodies is an exciting, but complex question. Here we consider solely the tidally induced surface temperature of a hypothetical moon. We were curious about the occurrence rate of moons with suitable surface temperature for life. For this reason we mapped the phase space evenly with hypothetical moons that have different radii (between 250 and 6550 km) and eccentricities (between 0.001 and 0.1), and their densities are that of the Earth. We used both the viscoelastic and the fixed $Q$ model for calculating the surface temperature of these bodies, and then calculated the percentage of those that have suitable surface temperature, i.e. that are located inside the TTZ ($0\leq T_{\mathrm{surf}}\leq 100~{}^{\circ}$C). The calculation was made for different orbital periods between 0.1 to 3.5 days, and for each value there were 63100 hypothetical moons distributed in the radius-eccentricity phase space. The result can be seen in Fig. 5.7. Red solid and green dashed curves indicate the percentage of being inside the TTZ for the viscoelastic and for the fixed $Q$ model, respectively. The blue dotted curve shows the percentage of those cases that do not give result for the viscoelastic model. The top axis shows the $\beta$ parameter, which is the ratio of the moon’s semi-major axis and the planet’s Roche radius. It can be clearly seen that the red and green curves have a peak, which means that the probability of having suitable surface temperature has a maximum at a certain orbital period. The viscoelastic model predicts a much more efficient heating than the fixed $Q$ model, i.e. a much larger fraction of the hypothetical moons have their surface temperature between 0 and 100 ∘C. The ratio of the integral under the red to that under the green curve is 2.8, meaning that 2.8 times more exomoons are predicted in the TTZ with the viscoelastic model. For the viscoelastic model the maximum percentage appears around 1 day orbital period, and here the probability for the moon of being inside the TTZ is almost 80%. For higher orbital periods, this probability rapidly falls down, which is in contrast with the fixed $Q$ model. The latter has its peak around 1.5 days, and it has less than 20% chance for satellites being in the TTZ. Despite the high probabilities achieved by the viscoelastic model for small orbital periods, the fixed $Q$ model give more promising results for those moons that have their orbital periods at 2 days or more.

For detailed study, the 0 and 100 ∘C temperature contours were plotted in the radius-eccentricity plane for a few, specific orbital periods, namely $P=0.5$ day (top panel), $P=1$ day (middle panel) and $P=1.5$ days (bottom panel) (see Fig. 5.8). Again, red and green colours represent the viscoelastic and the fixed $Q$ model, respectively. Between the contour curves the region of the TTZ is filled with light red and light green colours. The result shows that the viscoelastic model mostly favours the small moons, especially at high eccentricities, but also some large moons at small eccentricities over the fixed $Q$ model. This suggests that the viscoelastic model is less sensitive to the parameters of the moon, and holds the temperature more steady than the fixed $Q$ model. This is due to the melting of the inner material of the moon that leads to a less elevated temperature, as discussed by Peters & Turner (2013). On the other hand, the lower temperature implies that the total irradiated flux of the moon will be also lower, hence making the detection of the moon more difficult.

We were also interested in ‘Earth-twins’ as satellites and in the probability of their ‘habitability’. For this reason we made similar calculations, but the radius and the density of the hypothetical moons were set to be close to that of the Earth: $R_{\mathrm{m}}=6378\,\mathrm{km}\,(\pm 5\%)$ and $\rho=5515\,\mathrm{kg/m^{3}}\,(\pm 5\%)$. The radius and density values were chosen randomly from these intervals. The eccentricity was altered similarly as in the previous case (uniformly between 0.001 and 0.1). Altogether 200000 cases were considered for each orbital period. The temperature limits were set to $-20$ and 60 ∘C, which are the probable limits for life on Earth. Fig. 5.9 shows the results of this calculation. Note, that the peaks of the red solid (viscoelastic model) and the green dashed (fixed $Q$ model) curves are shifted to higher orbital periods, comparing to Fig. 5.7. This is in part caused by the changed temperature limits, and in part by the much shorter radius range. The maximum probabilities are also higher, that is especially visible in the case of the fixed $Q$ model, which reaches more than 40% at the curve’s peak (in the previous case it was less that 20%). As one would expect, it suggests that larger moons are more probable of maintaining warm surfaces. The ratio of the areas under the red and the green curves is 2.3.

In general, it can be concluded that the viscoelastic model is not just more realistic than the fixed $Q$ model, but also gives more promising results for exomoons, since much larger fraction of the hypothetical satellites have been found in the TTZ. In those cases when the viscoelastic model does not give solution for the equilibrium temperature, one can use the fixed $Q$ model instead, however, the values of $Q$ and $\mu$ are highly uncertain.

5.2.4 The value of Qµ

With the product of $Q$ and $\mu$, one can easily calculate the tidally induced surface temperature of a moon without using a complex viscoelastic model. Using Eq. 5.1 is a fast way to obtain $T_{\mathrm{surf}}$, but a good approximation is needed for the $Q\mu$ value. For such calculations, in the following, we give the $Q\mu$ values for hypothetical moons. Because of the large number of possible variations in the physical and orbital parameters of the moons, only a few, Solar System-like bodies were considered. Since the $Q\mu$ varies several orders of magnitudes for different rocky bodies, a good estimation can serve almost as well as the exact value. The following examples can be used as a guideline for making such estimations. Note, that the used model can be applied to rocky bodies, but for icy satellites, such as Enceladus or Europa, the results may be misleading, because of the more complex structure and different behaviour of icy material.

From the surface temperature of the moon, that was calculated from the viscoelastic model, the $Q\mu$ value was determined using Eq. 5.1 for six orbital period–eccentricity pairs. In Table 5.1 the tidally induced surface temperatures and the logarithm of the $Q\mu$ values can be seen. The radius and density of the satellites are those of the corresponding Solar System bodies (see the first column), and the values are from Murray & Dermott (1999). The eccentricities are set to 0.01 and 0.1, and the orbital periods to 1, 2 and 3 days. ‘N.a.’ indicates that there was no solution (weak tidal forces).

5.2.5 Scaling the Galilean satellite system

Since no satellite has been discovered so far outside the Solar System, we used the Galilean system as a prototype for realistic calculations. Io, Europa and Ganymede are orbiting in a 1:2:4 mean motion resonance, that maintains their eccentricities, which play an essential role in forcing continuously their tidal heating. Ogihara & Ida (2012) investigated satellite formation in the circumplanetary disc of giant planets using $N$-body simulations including gravitational interactions with the circumplanetary gas disc. They have found that 2:1 mean motion resonances are almost inevitable in Galilean-like satellite systems, and based on their results they predict that mean motion resonances may be common in exoplanetary systems. For these reasons the Galilean satellite system seems to be a representative example for realistic calculations, since the moons are in resonance, and their scaled-up versions will probably stay in resonance, too.

The test systems consist of a planet (Jupiter) and the four Galilean moons. 91 cases are considered, one is the real Galilean system, and the others are the scaled-up versions: the masses of the planet and the moons were multiplied by the $scale$ factor ($scale=1.0,1.1,1.2,...\,10.0$), and the semi-major axes of the moons were altered with constant orbital periods for each $scale$ value.

[TABLE]

where $a$ is the semi-major axis of the moon, $M_{\mathrm{p}}$ and $M_{\mathrm{s}}$ are the masses of the planet and the moon, respectively. The fixed orbital periods guarantee that the satellites approximately stay in resonances. This calculation resulted in constant $\beta$ values for all $scale$ parameters.

Using both the fixed $Q$ and the viscoelastic models, the warmth of tidal heat was investigated in each case. The tidal heat induced surface temperature can be seen in Fig. 5.10, where the 91 cases are connected with solid curves for each satellite. 182 other cases were calculated, too, which are shown with dashed and dotted curves in the figure. These curves indicate that the densities (dashed curve) and the eccentricities (dotted curves) of the satellites are doubled compared to their original values in the Solar System. In the calculations $\mu$ and $Q$ were set for all satellites to that of Io, namely $10^{10}\mathrm{\,kg/(m\,s^{2}})$ and 36, respectively, except for Europa, which has the following parameters: $\mu=4\cdot 10^{9}\mathrm{\,kg/(m\,s^{2}})$ and $Q=100$ (Peters & Turner, 2013). Densities of the moons are from Lodders & Fegley (1998), and the reference for the semi-major axis, eccentricity and mass values is Murray & Dermott (1999, Appendix A).

For Io, in the $scale=1$ (Solar System) case the fixed $Q$ and the viscoelastic models give 60 K and 160 K, respectively. The observed surface heat flux induced by tidal heat on Io is around $2\mathrm{\,W/m^{2}}$, which is a lower limit (Spencer et al., 2000). In other words, tidal forces produce at least 77 K heat on the surface of Io. The fixed $Q$ model resulted in a lower value than this limit, but note that $Q$ and $\mu$ are very difficult to estimate. The viscoelastic model gave much higher temperature than the observation, but keep in mind, that the heat is concentrated in hotspots, and is not evenly distributed on the surface of Io. The temperature of the warmest volcano, Loki is higher than 300 K (Spencer et al., 2000).

The viscoelastic model gives solution only for Io (orange curves) and Europa (light blue curves), but not for all $scale$ values, as shown in Fig. 5.10. In those cases when the densities are twice than those in the Solar System (dashed curves), the surface temperatures are just slightly higher. In fact, in the viscoelastic model, higher densities result in less tidal heat because of the imaginary part of the second order Love number. Doubling the eccentricity instead of the density (dotted curves) makes the surface temperature higher in each case.

5.3 Tidal heating with climate models

For a complete habitability investigation of exomoons other energy sources are also needed to be considered besides tidal heating. Although it is possible that stellar radiation can be neglected (either because there is no central star, or because the planet-moon system orbits the star in a far distance), a general description must include the irradiated energy, as well. For this reason a complex climate model was used as described in section 3.2.3, and the viscoelastic tidal heating model was applied as a part of the whole calculation. From the viscoelastic model the $Q\mu$ value is obtained as described in section 5.2.4. From this value the surface heat flow ($\zeta$) was calculated and then inserted into Eq. 3.21. The stellar mass in the simulations is $1M_{\odot}$, the mass of the planet and the moon are $1M_{\mathrm{Jup}}$ and $1M_{\oplus}$, respectively. The planet revolves the star on a circular orbit at $a_{\mathrm{p}}=1$ AU distance.

Fig. 5.11 shows the results of the simulation data as functions of the semi-major axis and eccentricity of the moon’s orbit. For the data in the left-hand column the fixed $Q$ tidal heating was applied (calculated with Eq. 3.22), compared to viscoelastic tidal heating model in the right-hand column. Different rows show different orbital inclinations of the moon relative to the planet’s equator. The four colours correspond to the four possible states of the satellite as described in section 3.2.3 and also summarised in the caption of the figure.

The circumplanetary HZ is considerably wider when the viscoelastic tidal heating is used (right-hand column of Fig. 5.11) compared to the fixed $Q$ model. While the weaker tidal heating at low $e$ tends to pinch the HZ slightly, at high $e$ the HZ is much wider, since tidal heating prevents catastrophic heating of the moon. For zero inclination the HZ is at least as wide as 0.06 Hill radius, and likely wider than that is measured in this parameter space. Perhaps most notably, for $a\lesssim 0.07$ Hill radius, the outer edge is not present at any eccentricity. These moons experience rapid, brief eclipses which are compensated by thermal inertia and the carbonate-silicate cycle.

If the inclination of the moon is increased (middle and bottom rows) than the outer edge begins to disappear, even at relatively large moon semi-major axis. For the outer edge to become apparent, eclipses must be sufficiently frequent for their effects to be cumulative and induce a snowball state. Increasing inclination reduces the eclipse rate (see e.g. Heller, 2012) until it is close to zero. (For polar lunar orbits ($i=90^{\circ}$), the eclipse rate is not exactly zero – eclipses can still occur twice per year, when the lunar nodes align with the planet’s orbital position vector.) With both implementations of tidal heating, relatively large inclinations ($45^{\circ}$) erase the outer habitable edge completely in the investigated parameter space. At intermediate inclinations ($10^{\circ}$), we see a slightly mixed picture. The eclipse rate should be effectively zero, yet the fixed $Q$ tidal heating run shows that an edge still exists (albeit much further from the planet). This is an indication that here the ice-albedo feedback mechanism is activating due to the moon’s distance from the star rather than due to eclipses, which in turn may suggest that the albedo prescription we use may in fact be too simple for this scenario.

5.4 Conclusions

We have used, for the first time, a viscoelastic model for calculating the surface temperature of tidally heated exomoons. The viscoelastic model gives more reliable results than the widely used fixed $Q$ model, because it takes into account that the tidal dissipation factor ($Q$) and rigidity ($\mu$) strongly depend on the temperature. Besides, these values are poorly known even for planets and satellites in the Solar System. Using the viscoelastic model for exomoons helps to get a more realistic estimation of their surface temperature, and to determine a circumplanetary region, where liquid water may exist on them. It may help future missions in selecting targets for exomoon detections.

We have defined the Tidal Temperate Zone, which is the region around a planet where the surface temperature of the satellite is between $0~{}^{\circ}$C and $100~{}^{\circ}$C. No sources of heat were considered other than tidal forces. Assuming, that the planet-moon system orbits the star at a far distance, or the stellar radiation is low due to the spectral type, tidal heat can be the dominant heat source affecting the satellite. We have investigated such systems, and found that the TTZ strongly depends on the orbital period, and less on the radius of the moon. For higher densities or eccentricities of the moon, the location of the TTZ is slightly closer to the planet.

Comparing this model to the traditionally used fixed $Q$ model revealed that there are huge differences in the results. Generally, the viscoelastic model is less sensitive to moon radius than the uniform $Q$ model, keeping the surface temperature of the body more steady. The reason is that higher tidal forces induce higher melt fraction which results in a lower temperature than the fixed $Q$ model. The viscoelastic model demonstrates the way in which partially melting of a moon can act as a thermostat and tend to fix its temperature somewhere near its melting point over a wide range of physical and orbital parameters. As a consequence, the statistic volume of the TTZ is much larger in the viscoelastic case, which is favourable for life. But this lower temperature also means that the detectability of such moons is lower in the infrared. In addition, for low tidal forces there is no equilibrium with the convective cooling; hence, only the fixed $Q$ model provides solution. In these cases the challenge is to determine the values of $Q$ and $\mu$.

For a few characteristic cases the product of the tidal dissipation factor and rigidity was calculated from the viscoelastic model, in order to help in quick estimations of tidally heated exomoon surface temperatures. Since the viscoelastic model is more realistic because of the inner melting and the temperature dependency of the parameters, but the fixed $Q$ model is easier to use, these $Q\mu$ values (along with the surface temperature) are provided in Table 5.1. By inserting $Q\mu$ into Eq. 5.1, one can get the estimation of the tidally induced surface temperature of a moon. Connection between the quality factor ($Q$) and the viscoelastic parameters (viscosity and shear modulus) was given for the Maxwell model, too, by Remus et al. (2012).

Earth-like bodies were also investigated as satellites, and in these cases the $-20$ and $60~{}^{\circ}$C temperatures were used as limits of habitability. The results are similar, but the volume of this habitable zone is larger than that of the TTZ for wide range of satellite radii. This habitable zone includes atmospheric considerations of the moon, but stellar radiation was neglected in the calculations. In case of significant radiation from other sources, the surface temperature of the moon will be higher. Additional heat sources (such as stellar radiation, radiogenic processes, reflected stellar and emitted thermal radiation from the planet), and the effects of eclipses, or the obliquity of the satellite are thoroughly discussed by Heller & Barnes (2013).

For simulating realistic systems, the Galilean moons were used as a prototype. Their surface temperature was calculated with both models for different, scaled up masses. The mean motion resonance between the satellites helps to maintain their eccentricity, and consequently to maintain the tidal forces. By raising their masses, the temperatures of Io and Europa elevates less drastically in the viscoelastic model, than in the fixed $Q$ model (see Fig. 5.10). At $scale=5$ (masses are five times as in the Solar System case) the surface temperature of Europa is $\sim 150$ K calculated from the viscoelastic model. Assuming that its density does not change, its radius will be approximately 0.25 Earth radii. In case of an additional 100–120 K heat (e. g. from stellar radiation), the ice would melt, and this super-Europa would become an ‘ocean moon’, covered entirely by global water ocean. The used viscoelastic model might not be adequate, and can be oversimplified for such bodies that consist of rocky and icy layers, as well. Salty ice mixtures may also modify the results. The applied model ignores the structure, pressure and other effects, and applies melting for the whole body. However, it provides a global picture of the tidally heated moon. Even with a more detailed viscoelastic model, that describes Enceladus as a three layered body (rocky core, ocean and ice shell), Barr (2008) have found that tidal heat is $\sim 10$ times lower than that was observed by the Cassini Composite Infrared Spectrometer. Similarly, Moore (2003) concluded that observed heat flux on Io is about an order of magnitude higher than that can be explained with a multilayered, viscoelastic model. These results suggest that tidal heat can be much more relevant than what is predicted by models.

We have also used a 1D energy balance model of Earth-like exomoon climates, which contain stellar and planetary insolation, atmospheric circulation, infrared cooling, eclipses and tidal heating as the principal contributors to the moon’s radiative energy budget. The model also takes into account the carbonate-silicate cycle, a negative feedback mechanism to regulate planetary temperatures and the positive ice-albedo feedback system. We calculated tidal heating as part of this whole climate model with both the fixed $Q$ and the viscoelastic models.

We have found that the results show many more potentially habitable configurations when the viscoelastic model is used, even in our somewhat limited parameter space. The circumplanetary HZ is significantly wider with viscoelastic tidal heating, and extends much further from the planet (provided the moon’s orbital eccentricity is larger than about 0.02, which is slightly less than that of Titan). We have found that in case of zero inclination of the moon’s orbit there is a well-defined outer edge to the circumplanetary HZ, inside the orbital stability limit, due to a combination of eclipses and ice-albedo feedback. The outer edge requires the orbits of the planet and the moon to be quite closely aligned, so that eclipses are sufficiently frequent. We have shown that if the moon’s orbit is inclined so that eclipses are unlikely, the outer edge completely disappears for any eccentricity of the moon’s orbit, even for relatively large semi-major axes of the moon.

6.   The possibility for albedo estimation of exomoons

The aim of this work is to propose a new method to identify elevated albedo and thus the possibility of water ice on the surface of an exomoon using photometric measurements during occultations. We will also discuss the role of different orbital and physical parameters that may influence the observed data. The following work is based on the article by Dobos et al. (2016).

Surface ice cover might be more frequent for exomoons than for exoplanets, based on examples of the Solar System and theoretical argumentation. Although the photometric signal would be small, it is worth considering the possibility and potential results of such measurements, because of the following reasons: the findings might orient the instrumental development, and search programs to focus on given stellar type, and also give a hint on surface properties of exomoons, or provide information on their habitability. Analysing different configurations helps to see the best possibility to estimate exomoons’ albedo. We also aim to compare the scale and possible role of such configurations of exomoons that are realistic but not present in the Solar System. Just like in the case of exoplanets, unusual situations could be present and occasionally dominate because of observational selection effects (see for example the case of many hot Jupiters during the first decade of exoplanet discoveries). In this work we analyse these possibilities, pointing to the difficulties, possible ways and best conditions for such estimation.

Although such precise measurements are difficult to make, the method could be used to differentiate among exomoon surface properties. Despite no exomoon has been discovered yet, probably those are numerous in other planetary systems, and will play an important role in astrobiology research in the future (Kaltenegger, 2010; Heller & Barnes, 2013) giving rationality to elucidate the potential method outlined here. Because of the large number of possible configurations and different parameters of exomoon systems, specific cases are considered in this paper with the aim to give general insights about the factors that are important for albedo estimation. We also consider such possible conditions that are not present in the Solar System but could be observed at certain exomoons, e.g. large exomoons, objects on eccentric, inclined or retrograde orbits (section 6.3.2).

The basic idea behind this work is that the exomoons’ albedo could be estimated from the flux change during occultation. High albedo values might suggest water ice on the surface (Verbiscer et al., 2013). With the analysis of periodic changes in the transit time and its duration (Kipping, 2009a, b), the mass of the moon can be calculated, from which the bulk density could also be estimated. The density of the moon helps in determining whether or not the moon contains significant amount of ice. For the albedo calculations Solar System moons and hypothetical bodies were considered as exomoons, in order to determine their observability during the occultation. The possible usage, problems and uncertainties of this method are described in section 6.3.

6.1 Solar-type stars

First, the flux change caused by the presence of different hypothetical exomoons was calculated using Eq. 3.32. As there might be a great variety of sizes, orbital distances and stellar luminosities, we restricted our analysis to certain number of example cases. For stellar distance 3 AU was used, which is the approximate location of the snowline in the Solar System. In the following calculations the Sun was used as central object, because future surveys and monitoring programs will probably focus on Sun-like stars to search for Earth-like exomoons (Kaltenegger, 2010; Peters & Turner, 2013). Beside solar type stars we also calculated some example cases for M dwarfs (see section 6.2.3), as they came into the focus of exoplanet related research recently (see the MEarth project, Berta et al., 2013).

A few parameters of interesting examples are shown in Table 6.1. Beside Solar System bodies (Earth, Europa, Enceladus, Io), ‘icy Earth’ is introduced as a hypothetical moon. Icy Earth is considered to reflect as much light as Enceladus (almost 100% of the incident light) and its radius equals to that of the Earth. It is important to note that Europa and Enceladus have no relevant atmospheres and thus we also consider the hypothetical ‘icy Earth’ to satisfy this criterion as well. Io does not contain water ice on the surface, but still has moderately elevated albedo because of fresh sulphur containing material. Although there is no Earth-sized moon in the Solar System, a natural satellite with this size might be considered realistic, and can serve as a useful for comparison to other objects (Ogihara & Ida, 2012). As it can be seen from Table 6.1, the required photometric precision is many orders of magnitude beyond the capability of current technology and such measurements will not be feasible in the foreseeable future.

6.2 Moons in systems of M dwarfs

6.2.1 Plausibility of large icy moons

Moons are more likely to form by ejection around Neptune than around Jupiter-like giants because. Jupiter-like gas giants are composed of mostly liquid and gaseous hydrogen at their outer layers, thus the ejected satellite-forming material might be $\mathrm{H_{2}O}$ poor, composed mostly of gaseous $\mathrm{H_{2}}$ and could easily escape. But if this event happens before the dispersion of the circumplanetary disk, the co-accreted ice in the disk could provide icy surface for such satellite. Neptune-like planets seem to be better candidates. They are mainly composed of $\mathrm{H_{2}O}$, and the ejected satellite, too. The occurrence of impacts at giant planets are favoured by certain theories, as a contributor to the continuous gas accretion (Broeg & Benz, 2012), increasing the probability of impact-ejected icy satellites.

Besides ejection, capture could also produce a large moon. In the so-called binary-exchange captures, a terrestrial sized moon can be captured by a 5 Jupiter-radii planet around an M dwarf, if the planetary encounter is close enough (Williams, 2013). For Neptune-mass planets the capture is easier than for Jupiter-mass planets, because generally the encounter speeds are lower. The capturing is much more difficult at small stellar distances, but if it occurs before or during the planetary migration, then it is possible that such a planet-moon pair will orbit an M dwarf at a close distance (around the snowline). In the Solar System, Triton may be an example for such binary-exchange capture (Agnor & Hamilton, 2006). For these reasons we consider the existence of large exomoons (up to a few Earth-masses) of Neptune-mass planets around M dwarfs.

6.2.2 Different wavelengths

It is well known that M stars are much brighter in the infrared than in the visible spectrum. For this reason it seems obvious to measure the flux at longer wavelengths, however, water ice is most reflective in the visible. In this section we discuss which wavelengths are optimal to use for successful measurements of icy moons.

Beside water ice, methane ice is also very reflective and such bodies are also known in the Solar System that are covered at least partially with methane ice. For this reason, if the estimated albedo is high, it may be challenging to decide whether water ice or methane ice is present on the surface. To make the distinction easier, see the differences in the spectra of Enceladus and Eris in the bottom panel of Fig. 6.1 (the spectra are from Verbiscer et al. (2006) and Alvarez-Candal et al. (2011)). Eris is a good example for a small body covered by methane ice, and Enceladus is used as a representation of a moon with water ice on the surface. The spectra are normalized, but these normalized values must be in good correspondence with the geometrical albedo for both bodies: the spectra of Enceladus are normalized at 0.889 $\mu$m where the geometric albedo (at phase angle ${\approx}3^{\circ}$) of the leading and the trailing hemispheres are 1.02 and 1.06, respectively (Buratti et al., 1998), and the reflectance of Eris is normalized at 0.6 $\mu$m (Alvarez-Candal et al., 2011), and its geometric albedo in the V band is 0.96 (Verbiscer et al., 2013) which is very close to one. From their spectra, the albedo values at different photometric bands are estimated and presented in Table 6.2. All values in the table are rough estimations, which can differ for various exomoons. The albedo of Enceladus (water ice case in the table) in the V band is from Verbiscer et al. (2007).

To see the albedo variations respective to wavelength, it is best to measure the flux in all photometric bands indicated in Table 6.2, but it seems that the V and the J bands are the most relevant. Measuring in these two bands provide useful data, partly because of the apparent difference in the albedo for the two kinds of ice, which helps determining whether water or methane ice is present on the surface. Another reason for the V band is the maximum reflectance of water ice, and for the J band is the maximum spectral radiance of low-mass M dwarfs. The radiance of an M6 – M7 main-sequence star in the J band is about 5 times stronger than in the V band. Because of the larger number of photons that arrive to the detector, it is easier to observe the occultation in the J band, where the albedo of water ice is still relatively high. However, the shape of the light curve will not change in different bands.

Clouds are also very reflective and for this reason it may be challenging to determine whether we see an icy surface or clouds in the atmosphere. Since water is very abundant is the universe, we assume that clouds form from H2O. The upper panel of Fig. 6.1 shows the intensity of water cloud on Earth (Gottwald et al., 2006, Chapter 7). In the F140M filter band (indicated by a green vertical band in the figure) the intensity of the clouds are very low while the reflectance of water ice is moderately high (see the Enceladus’ spectrum in the bottom panel). This photometric band gives opportunity to separate clouds from ice.

6.2.3 Next generation telescopes and the occultation of exomoons

It seems easier to estimate the albedo of satellites around M dwarfs than around larger stars. This is because the fact that the lower mass the star has, the closer the snowline is located. It means that the planet-moon system can be closer to the central star without the sublimation and/or melting of the ice on the surface. If the moon is closer to the star, then the light drop caused by its presence will be more significant in the light curve. This effect favours detection. In other words, it may be easier to observe icy moons around lower mass stars in occultation, than around more massive stars. It also means, that smaller moons may also be detected in occultation around late type stars.

The capabilities of next generation telescopes are of fundamental importance in successful albedo estimations. Beside Kepler, the most relevant missions for exomoon detection are the CHaracterising ExOPlanets Satellite (CHEOPS, Broeg et al., 2013), the Transiting Exoplanet Survey Satellite (TESS, Ricker et al., 2010, 2014), the James Webb Space Telescope (JWST, Clanton et al., 2012), the PLAnetary Transit and Oscillations of stars (PLATO 2.0, Rauer et al., 2011, 2014) and the European Extremely Large Telescope (E-ELT, Ramsay et al., 2014). The photon noise levels of these instruments were calculated as described in section 3.3.2. The data used for the calculations for each instrument are shown in Table 6.3. All these numbers are subject to change, since most of these missions are in the planning and testing phase. The numbers in the last column indicate references as explained in the footnote1111: Ricker et al. (2014), 2: Auvergne et al. (2009), 3: Broeg et al. (2013), 4: The CHEOPS Study Team (2013), 5: Rauer et al. (2014), 6: Van Cleve & Caldwell (2009), 7: Avila et al. (2016), 8: Dressel (2016), 9: JWST webpage: http://www.stsci.edu/jwst/instruments/nircam/instrumentdesign/filters/, 10: Davies et al. (2010).

[FIGURE:]

Figs. 6.2 and 6.3 show the relation between the MO depth of icy moons and the mass of the host star in the V and J photometric bands, respectively. The objects covered in this figure are low-mass M dwarfs: M9 – M5 according to Kaltenegger & Traub (2009). The MO depth is in logarithmic scale. The radius of the moon indicated by black curves was calculated for all stellar mass and MO depth pairs. For this calculation Eq. 3.32 was used and the albedo was fixed to that of Enceladus in each case, i.e. 1.38 for Fig. 6.2 and 0.84 for Fig. 6.3. The distance from the star was set to the snowline, which was calculated from Eq. 3.45, and is indicated in the upper axis. This is the minimum stellar distance that is required in order to prevent the melting of the surface ice. It also means that our calculation is an optimistic approach: the minimum distance and the approximately minimum moon radius are determined for a successful detection. The vertical grey lines indicate the 60, 90 and 120 minute MO durations. In all investigated cases, the moon and the planet occults separately. For larger stellar masses the MO duration is longer, because of the larger stellar distance. The expected noise levels for the telescopes are indicated by coloured curves. Those instruments that measure in or near the J photometric band (HST WFC3 IR, JWST NIRCam and E-ELT MICADO) are shown in Fig. 6.3, all the others are in Fig. 6.2, however, because of their limited capabilities, TESS, and CoRoT do not appear in Fig. 6.2.

The highest radius that is shown in Figs. 6.2 and 6.3 is 2.5 Earth-radii (for comparison, the radius of Neptune is 3.9 Earth-radii). We choose this 2.5 Earth-radii as an upper limit for a rocky moon with surface ice to guarantee that its mass is not larger than 5 Earth-masses. (The radius of a body with the density of Enceladus and with a mass of 5 Earth-masses is 2.58 Earth-radii.) Such large body around a Neptune-mass planet could also be called a binary planet, which would probably be a different system than the Saturn–Enceladus planet–moon pair. For simplicity, and also because we use this 2.5 Earth-radii as an upper limit, we still call them moons in the rest of the paper.

In the visible spectrum, the photon noise of the telescopes is much higher than the required precision for detecting large, rocky moons in occultation. However, in the J photometric band (see Fig. 6.3) the E-ELT’s photon noise is just a few ppm greater than the MO depth of large icy exomoons. (Note that the MO depth is shown in a logarithmic scale.) For example, for a 2.5 Earth-radii moon around an 0.11 stellar mass star, the difference is about 3-4 ppms, and for a 1.5 Earth-radii moon at an 0.085 stellar mass star, it’s about 7 ppms. The accuracy of the measurements and additional parameters that can influence the success of the observations are discussed thoroughly in the next section.

6.3 Discussion

The findings on the observability of different exomoons, the estimation of their surface albedo, and other related considerations are summarised below.

• •

Smaller stellar distance: The depth of the flux drop caused by the moon in occultation as well as its observability increases with decreasing stellar distance. In the case of smaller stars, the surface of the moon can keep its icy composition even for small stellar distances, because the star’s luminosity is low. Small stellar distances (such as 0.01–0.1 AU in our calculations) are important regarding the long term stability of a planet-moon system, too. According to migration theories and perturbation effects, numerous criteria have to be satisfied in order to support the presence of an exomoon that survived a certain time period after it got close to the central star (Barnes & O’Brien, 2002). During this process they might be able to keep their satellite system protected against perturbations by the central star and other planetary objects (Mosqueira & Estrada, 2003). However, the migration is not relevant, if the moon is captured or impact-ejected (which seems more likely for large satellites) after the planet’s migration phase.

• •

Larger size: Among the analysed cases, the best possibility to determine the albedo is for Earth-sized, or even larger exomoons, because the reflecting surface is larger.

• •

Higher albedo: The higher the albedo that an exomoon has, the more stellar light it reflects, and the flux difference in the light curve is larger during the occultation. Extremely high geometric albedos (for example 1.375 in the case of Enceladus in the Solar System) are indicators for the presence of relatively young and fresh water ice since this is the only known material that could reflect such high percentage of light. The presence of water ice on the surface may imply astrobiological importance. Considering strongly reflecting atmospheric clouds, their presence could also cause high albedo. The separation of cloud and ice reflection with optical methods could be difficult, but in case of successful measurements with an F140M filter, it may be feasible, as discussed in section 6.2.2. However, our results show that such measurements will not be achievable with the instruments presented in section 6.2.3, meaning that next generation telescopes will not yet be able to separate atmospheric clouds from surface ice on exoplanetary bodies.

Beside the challenge in observation of the small flux contribution of the exomoon, analysing the measured data will also be difficult, as to fit a realistic model, because many possible configurations are needed to consider. To see the possible roles of these factors, and in order to help to plan and target future observations, in the following we evaluate how several basic parameters affect the observation and analysis.

6.3.1 Accuracy

For an adequate result from the light curve measurements, repeated observations are needed. It is not necessary that the moon is positioned exactly at the same apparent position every time when the occultation occurs. For this reason the small flux drop caused by the moon may be shorter or longer, and may not be seen at every observation. Simon et al. (2012) propose using an averaged light curve of several observations to deal with the scattered position of the flux change caused by the presence of the moon. They also present a four-step detection strategy for discovering exomoons by transit method. Occultation events with the same apparent positions of the moon and the planet might be rare, but in tidally locked situations, a given configuration during subsequent occultations could be frequent or regular and repeated at every revolution of the exoplanet, and the exomoon might be at the same position viewed from our line of sight.

Variable surface features of the host star, such as granules, micro-flares, etc. could influence the result, producing small flux fluctuations. Such variability is expected especially in the late type stellar atmospheres related to convective motion (Ludwig, 2006). From the moon’s transit measurements and from the occultation of the host planet, the expected time of the moon’s occultation can be estimated, which may help in distinguishing the MO signal from the stellar surface fluctuations, especially because the occultation will cause a periodic change in the light curve.

Beside the uncertainties in flux measurements, the accuracy of albedo estimation also depends on the accuracy of the exomoon’s size and the stellar distance. In our work both are assumed to be already known from transit detections, and determining their accuracy was already discussed by several authors (see e.g. Pál, 2008; Carter et al., 2008). However, the errors from transit measurements propagate into the determination of the geometric albedo. The total error of the albedo can be obtained from

[TABLE]

where the three terms of the right side of the equation can safely be considered to be independent of each other. The first term is $(1/5)^{2}$ since we set the measurements to $5\,\sigma$. The second term is typically around (1–5%$)^{2}$ (Pál, 2008). Assuming that $a_{m}/R_{\star}\geq 1$ (i.e. the moon and the planet occult separately, which was true in all cases described in section 6.2.3), and that the same number of occultations are detected for the moon as for the planet, and that these measurements were made with the same instrument, the third term can be expressed as

[TABLE]

where $Q$ is the quantum efficiency. This error can be calculated for given geometric configurations. In any case, the last term must be smaller than the first one in Eq. 6.1, because the S/N ratio is much larger for transits than for occultations, hence transit measurements are much more precise. See e.g. the case of CoRoT-1b which was measured in transit and occultation, as well (Snellen et al., 2009). The same effect is true for exomoons. It means that the first term of Eq. 6.1 is dominating in the error of the albedo estimation.

In the calculations shown in section 6.2.3 only photon noise was considered based on the quantum efficiency or on the system throughput. However, in near infrared measurements other noise sources become significant, as well: CCD read-out, CCD dark signal, Earth’s sky variations and zodiacal light. The produced noise will be comparable to the photon noise. The system throughput that was considered in the calculations (even when the quantum efficiency was used, since in these cases the quantum efficiency was multiplied by $0.8$) already contains some of these extra noise sources, but in future missions these are subject to change, since their precise values will be known only after the commissioning of the instruments. In addition, we assumed 30 occultation detections for observatories that would require at least 5 years for E-ELT because the orbital period of the planet at the snowline is approximately 31 days, and the occultation can only be measured at nights. Considering that not all noise sources were included in our estimations, and that measuring 30 occultations may be unrealistic in most cases, it can be concluded that our calculations are too optimistic, and the real noise level of the instruments will be higher than presented.

According to Simon et al. (2007), measuring the photometric transit timing variations (TTVp) leads to the estimation of the size of the exomoon. By analysing transit duration variations (TDVs) as well, rough mass estimation is also possible under favourable conditions (Kipping, 2009a, b). Using the derived mass and radius together, information on bulk density could be gained, and putting together this value with the presence of surface ice gives strong hint on the possible existence and ratio of water ice in the interior.

6.3.2 Considering different parameters

In the following we categorize parameters by their importance in influencing the results of the observations. Specific cases that might be relevant for certain exomoons, but are not usual in our Solar System are discussed below.

• •

Eccentricity: For the analysed, relatively large exomoons we do not expect highly eccentric orbits as the large satellites in the Solar System (even the captured ones) have low eccentricity. However, perturbations might elongate the orbit, or resonances with other moons may maintain a higher eccentricity. The variable orbital velocity around the host planet influences the apparent speed of the exomoon perpendicular to our line of sight. The changing velocity influences the length of the plateau/occultation duration. Shorter plateau/MO duration means less integration time for the telescope which results in larger photon noise. The moon-star distance periodically change with the orbital phase of the moon, influencing its brightness, but this effect is very small. Eccentricity of the moon may enhance this phenomenon. It could be substantial in cases when the ratio between the stellar distance and the exomoon’s semi-major axis is relatively small. Sato & Asada (2010) described the effect of the satellite’s eccentricity on the transit light curve in details, which can similarly be used to occultations with some changes.

• •

Inclination: Large tilt of the orbital plane of the exomoon relatively to the planet’s orbital plane around the host star could also influence the plateau / occultation duration, especially when the exoplanet orbits at large distance from the host star. In case of large inclinations, instead of occulting, the moon may pass beside the star. The length of the plateau/occultation duration also depends on the inclination and the ratio of the exomoon’s apparent distance from its host planet and the diameter of the host star. Larger stellar disk and smaller orbital distance increase the influencing effect of higher inclinations, however if the stellar disk is smaller, then the moon will not occult for high inclinations. For transits the effect of the exomoon’s inclination is discussed by Sato & Asada (2010).

• •

Reflected light from planet: In ideal cases the exomoon’s brightness might be elevated by reflected light from its host planet. The strongest reflected light flux reaches the exomoon during the ‘new moon’ phase (watching from the exoplanet). In different configurations mutual transits and mutual shadows affect the shape of the light curve (Cabrera & Schneider, 2007). In the albedo estimation this effect of the reflected light from the planet is not significant, as the exomoon shows almost full moon phase during the occultation (looking from the Earth). This reflected light flux is negligibly small comparing to the stellar irradiation (Hinkel & Kane, 2013). The illumination from the planet with other light sources is thoroughly discussed by Heller & Barnes (2013).

• •

Orbital direction: Retrograde orbits affect the plateau/occultation duration as stellar eclipses by the host planet are more frequent, than with a prograde orbit and can be prolonged because of the large size of giant planets (Forgan & Kipping, 2013). Such eclipses may occur just before or after the occultation with the star, changing the shape of the light curve. Using the phase-folded light curve of several observations may help in filtering out such cases. Beside eclipses, it does not influence the observability whether the exomoon is approaching to or receding from the host star viewing from the Earth.

• •

Effect of hemispherical asymmetry: An exomoon that faces with different sides toward the star (and the toward the Earth) might have different brightness at different occultations, if significant albedo difference exists between the two hemispheres. Calculating with a Solar System analogy, Iapetus shows huge difference in its Bond albedo with $\sim$0.31 at the trailing, and maximum 0.1 at the leading hemisphere (Howett et al., 2010). The difference between the observations of the two hemispheres produces such brightness difference that can not be neglected. This effect may influence the result of the albedo estimation, but using the phase-folded light curve such flux changes may be noticed, and from the average light curve the average albedo can be estimated.

Based on this subsection, different parameters are categorized by their importance in occultation observations. The results can be seen in Table 6.4. The most relevant factors influencing the observation are: 1. close stellar distance that is observationally favourable especially at small main sequence stars; 2. star/planet diameter ratio: the larger the ratio is the longer MO duration could be observed which helps in reaching smaller photon noise for the measurement (only relevant at small stars); 3. larger exomoons are easier targets, and are favourable to form by impact ejection or by capture; 4. larger planetary distance of the moon can increase the length of the plateau duration; 5. higher albedo means that the satellite is brighter, hence makes the observation easier.

6.3.3 Importance of ice

Icy moons with subsurface oceans might hold large volume of liquid water for geologic time-scales and thus may potentially be habitable or at least could be considered as favourable environments for life. In the Solar System such subsurface oceans exist in Europa, Titan, Enceladus, and possibly also in Ganymede and Callisto. The liquid phase state of these oceans are maintained by tidal heat production, freezing point depressing solved ingredients, and/or internal heat due to accretional energy conservation and radiogenic heat (Carr et al., 1998; Khurana et al., 1998; Kargel et al., 2000; Zimmer et al., 2000; McCord et al., 2001; Schenk, 2002; Collins & Goodman, 2007; Roberts & Nimmo, 2008; Lorenz et al., 2008; Postberg et al., 2009; Iess et al., 2014; Dobos & Turner, 2015). Having information on the existence of surface ice, using mass estimation from TTV and TDV might increase the probability to identify exomoons that are probable of having subsurface oceans.

The brightness of the ice could point to relatively young surface and active resurfacing processes. In the case of Europa the ocean gives possibility for active resurfacing thus for the existence of clean and bright ice there, although locations with non-ice ingredients can also be found on Europa, while on Enceladus the subsurface liquid water contributes to the geyser-like eruptions, and to the increase of surface albedo caused by fresh water ice crystals that fall back to the surface (Kadel et al., 2000; Fagents, 2003; Porco et al., 2006; Spencer et al., 2006; Verbiscer et al., 2007; Howett et al., 2011). Although other reasons could also produce elevated albedo, the very high albedo is still a strong indication of the water ice on the surface, and even of a liquid subsurface ocean.

Estimation of moons’ albedo is much easier at M type stars, if the planet and the moon are close to the snowline. However, strong, regular flares and coronal mass ejections are expected, that potentially make the orbiting rocky bodies uninhabitable, because of high UV irradiation and loss of the atmosphere (Scalo et al., 2007). Exomoons still might be habitable because the ice cover attenuate UV radiation, and serve as a physical screening mechanism (Cockell, 1998; Cockell et al., 2000). The loss of the atmosphere could also easily happen at small stellar distance, but does not influence directly the oceans beneath ice sheets.

6.4 Conclusions

In this work we proposed a method for the first time that hints on the surface albedo of an exomoon and thus might indicate the existence of water ice. The argumentation presented in this work can be useful to orient future research and to identify the best candidate systems for such observations. The methods and the results of our calculations are presented firstly, then secondly those parameters are listed which substantially influence such measurements and also those that are not relevant – as such information is useful in planning and targeting observations.

Moons of smaller stars seem to be easier to detect, if their stellar distance is smaller. In addition, the smaller the star, the closer the snowline located, meaning that rocky bodies can stay icy even if they are orbiting the star in a closer orbit (but still outside the snowline). The vicinity to the star makes the body brighter, thus an icy exomoon is easier to detect in occultation close to the snowline of an M dwarf, than at the snowline of a solar-like star. The largest flux difference is expected from large, $\sim$ 2.5 Earth-radii icy moons (about 5 Earth-mass body with similar density of Enceladus) orbiting small M dwarfs ($\sim$ 0.07–0.13 solar masses) close to the snowline. We have found that such moons cannot be observed in occultation with next generation space missions in the visual spectrum, because the instruments’ photon noise is far greater than the flux drop caused by the moon. However, in the near infrared (J band), the E-ELT’s photon noise is just about 5-6 ppm greater than the MO depth, despite the lower albedo of ice. This result, however, is too optimistic, because it does not take into account other noise sources, such as CCD read-out and dark signal, which become significant in the near infrared measurements. Considering that not all noise sources were included in our estimations, and that measuring 30 occultations may be unrealistic in most cases, since it would require years to achieve, it can be concluded that the real noise level of the instruments in the near infrared wavelengths will be higher than presented in section 6.2.3. Flux fluctuations caused by stellar activity makes the measurements even more difficult. We conclude that occultation measurements with next generation missions are far too challenging, even in the case of large, icy moons at the snowline of small M dwarfs.

We outlined the parameters that can influence the detection. Based on assumptions, the orbital and physical parameters of exomoons might be highly diverse, thus we evaluated a range in the space of several parameters. We categorized these parameters by their influence on the success of MO detections (see Table 6.4).

We also discussed the possible properties of the most suitable exomoons for characterization which may be useful once the instrumentation is available. Exomoons around Neptune-sized exoplanets of small stars seem to be the best targets for occultation-based albedo estimation opposite to Jupiter-like exoplanets, if the satellites, formed by impact ejection or capturing, are large enough. In the case of such satellites the surface ice cover might be co-accreted if the ejection happened before the clearing up of the circumplanetary disk, or formed by ejection from planets composed mostly of water ice, like Neptune and not of hydrogen, like Jupiter. Binary-exchange capture is also a possibility for a giant planet to have a large icy satellite (neptunes are better candidates for this process than jupiters).

Based on our work, the first albedo estimations of exomoons from occultations are expected around low mass M dwarfs at small stellar distances. These moons will probably be large, and have high albedo which may imply the presence of ice on the surface.

7.   Summary

In this work I investigated the habitability of exoplanets and their satellites. Three major topics were described in details: the habitable zone, tidally heated exomoons and the possibility of albedo estimation of icy satellites. These topics cover not only circumstellar and circumplanetary habitability, but also connect observation with potentially habitable moons.

The habitable zone is widely used to investigate the habitability of planets around different stars. This zone is determined as a region around a star in which an Earth-like planet could support liquid water on its surface for a long period of time. For calculating the location of habitable zones around different stars, usually complex climate models are used which take into account the atmospheric properties of the Earth and the carbonate-silicate cycle, which regulates the temperature of the planet. Beside this complex climate model, simple equations can be used, too. These simple formulae are fitted to the results of the more complex climate models. For this simple calculation the stellar temperature and luminosity are needed as input parameters, which can be calculated from the stellar mass. However, different calculation methods give different results.

We applied three previously defined models, and we also proposed new empirical formulae for calculating the temperature and luminosity from the stellar mass. We numerically compared the results to measured stellar parameters, by calculating the reduced chi-square values. The new equations result in a better reproduction of the boundaries of the HZ (calculated from measured stellar parameters) than the other models.

Beside exoplanets, their moons can also be habitable. Tidal heating may play a key role in their habitability, since the elevated temperature can melt the ice on the body and prevent the snowball state even without significant stellar radiation. For calculating the tidal heating rate in exomoons, usually the so-called fixed $Q$ models are used, which highly underestimate the tidal heat of the body. These models use $Q$ and $\mu$, the tidal dissipation and rigidity (respectively) as constants, however these parameters are functions of the temperature. In addition, the correct values of $Q$ and $\mu$ are extremely difficult to estimate even for Solar System bodies. Their values can vary a few orders of magnitude for similar bodies, which means that their estimation for yet unknown exomoons is not possible, one can only guess or assume some value for these parameters. For these reasons, we applied a viscoelastic tidal heating model for exomoons for the first time. This model is more realistic than the widely used, fixed $Q$ models, because it takes into account the temperature dependence of the tidal heat flux and the melting of the inner material. Because of the temperature dependency of the shear modulus and viscosity that are used in the viscoelastic model, the tidal heating will also be a function of the temperature. The tidally heated surface temperature of exomoons is calculated for the equilibrium temperature of the body, which corresponds to the stable equilibrium between tidal heating and convective cooling. In the calculations no stellar radiation, or atmosphere were considered.

Using the viscoelastic model, we introduced the Tidal Temperate Zone, which is a circumplanetary region where the tidally induced surface temperature of the satellite is between 273 and 373 K. We have found that the location of the TTZ strongly depends on the orbital period of the moon and less on its radius. We applied the calculation method to different densities and orbital eccentricities of the moon. In some cases there was no solution (no equilibrium temperature), because the tidal forces were so weak that they could not induce convective cooling inside the body. We compared the results with the fixed $Q$ model and investigated the statistical volume of the TTZ using both models. We have found that the viscoelastic model predicts 2.8 times more exomoons in the TTZ with orbital periods between 0.1 and 3.5 days than the fixed $Q$ model for plausible distributions of physical and orbital parameters. The viscoelastic model provides more promising results in terms of habitability because the inner melting of the body moderates the surface temperature, acting like a thermostat.

For Earth-like bodies other habitability limits than the 0 and $100\,^{\circ}$C were used, as well. The lower limit for microbial activity in salty solutions is around $-20\,^{\circ}$C, and the upper limit for complex eukaryotic life is approximately $60\,^{\circ}$C, which also corresponds to the runaway greenhouse limit on Earth. For most of the calculations these probable limits of habitability were also considered when applicable ($\rho=\rho_{\mathrm{Earth}}$ and $R_{\mathrm{m}}\approx R_{\mathrm{Earth}}$), because these are related to biological and atmospherical constraints. When investigating the volume of habitable orbits, we found that the results were similar to the ones of the TTZ, but the volume was larger for Earth-like moons.

We have also applied a 1D energy balance model of Earth-like exomoon climates, which contain stellar and planetary insolation, atmospheric circulation, infrared cooling, eclipses and tidal heating as the principal contributors to the moon’s radiative energy budget. The model also takes into account the carbonate-silicate cycle and the positive ice-albedo feedback system. We used both the viscoelastic and the fixed $Q$ models for calculating the tidal heating inside the moons. We have found that the circumplanetary HZ is significantly wider with viscoelastic tidal heating, and extends much further from the planet, than with the fixed $Q$ model. In the case of zero inclination of the moon’s orbit relative to the planet’s equator there is a well-defined outer edge to the circumplanetary HZ, due to a combination of eclipses and ice-albedo feedback. We have shown that if the moon’s orbit is inclined so that eclipses are unlikely, the outer edge completely disappears for any eccentricity of the moon’s orbit, even for relatively large semi-major axes of the moon.

Icy moons might have subsurface oceans, such as Europa or Enceladus. If the moon is tidally heated, then life might emerge at the bottom of the ocean, where the rocky seafloor interacts with water through hydrothermal vents. In order to separate icy satellites from other exomoons, we propose a new method for obtaining their albedo. The albedo may be estimated from photometric measurements during occultations, when the reflected light from the satellite is blocked out and a small flux drop occurs in the light curve. The depth of the flux drop is proportional to the albedo and the radius of the moon.

If the albedo is high, it implies that the surface is covered with ice, however methane ice is very reflective, too. In order to determine whether the surface ice is of water or methane, we suggest to make measurements in several photometric bands, especially in the J band, where the two kinds of ices show different spectral features, but they are still very reflective. Clouds can also cause high albedo on a planetary body. Since water is very abundant in the Universe, we assume that water clouds might be common, and we compared the spectrum of the Earth’s water cloud with ice spectrum. We found that the cloud’s reflectance significantly drops in the F140M photometric band, where the ice’s albedo is still very high. In order to distinguish clouds from surface ice, we propose to make measurements in this band, too, upon the instrumentation is available for such precise observations.

We applied simple calculations for different stellar masses in the V and J photometric bands, and compared the flux drop caused by the moon’s occultation and the estimated photon noise of next generation missions with $5\,\sigma$ confidence. We found that albedo estimation by this method is not feasible for moons of solar-like stars, but small M dwarfs are better candidates for such measurements, because their weaker stellar radiation brings the snowline closer to the star, and if the planet-moon system orbits the star at a closer distance, then they will reflect more light. Our calculations in the J photometric band show that E-ELT MICADO’s photon noise is just about 4 ppm larger than the flux difference caused by a 2 Earth-radii icy satellite in a circular orbit at the snowline of an 0.1 stellar mass star. Because of other noise sources that were not included in our calculations, we conclude that occultation measurements with next generation missions are far too challenging, even in the case of large, icy moons at the snowline of small M dwarfs.

We also discussed the role of different parameters (some of them were neglected in the calculations, e.g. inclination, eccentricity, orbiting direction of the moon), and categorised them by their influence on the measurements. The most substantial influencing factors are the stellar distance, the size of the exomoon, its albedo, the star-planet size ratio and the planetary distance. We also predict that the first albedo estimations of exomoons will probably be made for large icy moons around the snowline of M4 – M9 type main sequence stars.

8.   Short description

I investigated the habitability conditions of extra-solar planets and satellites.

We proposed new empirical equations for calculating the temperature and luminosity from the stellar mass. By using numeric investigations, we have found that these formulae reproduce the stellar parameters more precisely than other models that we applied. With the empirical equations, the boundaries of the circumstellar habitable zone fit well to the ones which were calculated from measured stellar parameters of F, G, K and M main sequence stars.

Tidal heating of exomoons may play a key role in their habitability, since the elevated temperature can prevent the body to enter a snowball state even without significant stellar radiation. We used a viscoelastic model for exomoons for the first time, which is more realistic than the widely used, so-called fixed $Q$ models. We introduced the circumplanetary Tidal Temperate Zone (TTZ), and compared the results to the fixed $Q$ model. We have found that the viscoelastic model predicts 2.8 times more exomoons in the TTZ, than the fixed $Q$ model, for plausible distributions of physical and orbital parameters. The viscoelastic model provides more promising results in terms of habitability because the inner melting of the body moderates the surface temperature, acting like a thermostat. We have also used a 1D energy balance model of Earth-like exomoon climates, and found that the circumplanetary HZ is significantly wider with viscoelastic tidal heating, and extends much farther from the planet, than with the fixed $Q$ model. We have shown that if the moon’s orbit is inclined, then the outer edge completely disappears for any eccentricity of the moon’s orbit.

Icy moons might also be habitable, if they have tidally heated subsurface oceans. Occultation light curves of exomoons may give information on their albedo and hence indicate the presence of ice cover on the surface. We compared the flux drop caused by the moon’s occultation with the estimated photon noise of next generation missions with $5\,\sigma$ confidence. We found that albedo estimation by this method is not feasible for moons of solar-like stars, but small M dwarfs are better candidates for such measurements. Our calculations in the J photometric band show that E-ELT MICADO’s photon noise is just about 4 ppm greater than the flux difference caused by a 2 Earth-radii icy satellite in a circular orbit at the snowline of an 0.1 stellar mass star. We conclude that occultation measurements with next generation missions are far too challenging. We predict that the first albedo estimations of exomoons will probably be made for large icy moons around the snowline of M4 – M9 type main sequence stars.

9.   Kivonat

Disszertációmban a Naprendszeren kívüli bolygók és holdak lakhatóságát vizsgáltam.

Új, empirikus formulák használatat javaslom a hőmérséklet és a luminozitás csillagtömegből történő meghatározására. Numerikus vizsgálatok alapján arra a következtetésre jutottam, hogy ezek a képletek pontosabban visszaadják a mért csillagparamétereket, mint a többi vizsgált modell. Az új egyenletek használatával a lakható zóna határaira kapott eredmények jól illeszkednek a mért paraméterek alapján számolt értékekhez F, G, K és M típusú fősorozati csillagokra.

Az exoholdak árapályfűtése kulcsszerepet játszhat a holdak lakhatóságában, hiszen az árapályerők elhanyagolható csillagsugárzás mellett is megemelik az égitest hőmérsékletét, ami megakadályozhatja a hólabda állapot kialakulását. Első alkalommal használtam viszkoelasztikus modellt az exoholdak árapályfűtésének leírására, ami realisztikusabb a széles körben használt, ún. rögzített $Q$ modelleknél. Bevezettem a bolygó körüli mérsékelt árapályfűtési zóna (MÁZ) fogalmát, és az eredményeket összehasonlítottam a rögzített $Q$ modellével. Arra jutottam, hogy észszerű fizikai és keringési paraméterek használata mellett a viszkoelasztikus modell 2,8-szer több exoholdat jósol a MÁZ-ban, így a lakhatóság szempontjából kedvezőbb eredményeket ad, mert a hold belsejének megolvadása mérsékli a felszíni hőmérsékletet. A földihez hasonló légkört feltételezve az 1D energiaegyenleg modell alkalmazásával arra jutottunk, hogy a bolygó körüli lakható zóna sokkal szélesebb a viszkoelasztikus modell használata esetén, és a bolygótól nagyobb távolságokig kiterjed, mint a rögzített $Q$ modellel. Megmutattuk azt is, hogy ha a hold pályasíkja dőlt, akkor a külső határ bármely excentricitás értékre eltűnik.

Jeges exoholdak is lehetnek lakhatók, ha a felszín alatti óceán árapályfűtött. Az exoholdak okkultációs fénygörbéje információt adhat az albedóról, amivel a felszíni jég jelenlétére következtethetünk. Összehasonlítottam a hold okkultációjakor bekövetkező fluxuscsökkenést a következő generációs távcsövek fotonzajával $5\,\sigma$ szórás mellett. Megállapítottam, hogy a Naphoz hasonló csillagok esetén az albedó becslése nem valósítható meg ezzel a módszerrel, de a kis M törpék ígéretesebb célpontok lehetnek. Megmutattam, hogy a J fotometrikus sávban az E-ELT MICADO fotonzaja csak 4 ppm-mel haladja meg egy 2 földsugarú jeges hold fluxuscsökkenését egy 0,1 naptömegű csillag hóhatáránál húzódó körpályán keringve. Arra a következtetésre jutottam, hogy az okkultációs mérések túlságosan meghaladják a következő generációs távcsövek képességeit. Az első exoholdas albedóbecsléseket várhatóan a nagy, jeges holdakra fogják elvégezni, melyek M4 – M9 típusú fősorozati csillagok hóhatárának közelében keringenek.

10.   Acknowledgements

I would like to thank the work and help that I got from my supervisors during the years. I express my grateful thanks to Prof. László L. Kiss, who helped with useful discussions and ideas that improved my research, and funded my participation at conferences and summer schools through his projects, which led to new collaborations and research works, part of them presented in this thesis. I am also grateful to Prof. Edwin L. Turner for his invitation to work at Princeton University, for his hospitality, and the helpful research consultations in person and also online across continents. I thank the support of Dr. Judit Orgoványi who sacrificed her free time to help me along the way to my PhD. I thank Prof. József Vinkó the help he provided with the recalculation and improvement of the habitable zone boundaries. I thank Prof. László Szabados for reading my dissertation, providing useful comments and grammar checking. I also thank all of the co-authors whom I worked with, without them I couldn’t have written my thesis. I am indebted to my parents and family for helping me through my education, which finally, after following a long and curvy road, led me to becoming an astronomer and hopefully to getting a doctorate degree in the near future. I owe my husband a debt of gratitude for being patient and very supportive even during the most stressful times. Special thanks to my friends, especially the (extended) flight team, for the happy moments that we spent together, which helped me to ‘recharge my batteries’ so I could continue working. Finally, I would like to thank the Eurest cantine in the Ericsson building in the Infopark, the bonappétit food service, and the Genomics Café in the Icahn building for the lunches, that saved me a lot of time, since I didn’t have to cook on weekdays.

I have been supported by the Hungarian OTKA Grant K104607, the Hungarian National Research, Development and Innovation Office (NKFIH) grant K-115709, the Lendület-2009 Young Researchers Program of the Hungarian Academy of Sciences, the TÉT-14FR-1-2015-0012 project, and the ESA PECS Contract No. 4000110889/14/NL/NDe.