Survival of exomoons around exoplanets

Despite numerous attempts, no exomoon has firmly been confirmed to date. New missions like CHEOPS aim to characterize previously detected exoplanets, and potentially to discover exomoons. In order to optimize search strategies, we need to determine those planets which are the most likely to host moons. We investigate the tidal evolution of hypothetical moon orbits in systems consisting of a star, one planet and one test moon. We study a few specific cases with ten billion years integration time where the evolution of moon orbits follows one of these three scenarios: (1)"locking", in which the moon has a stable orbit on a long time scale ($\gtrsim$ 10$^9$ years); (2)"escape scenario"where the moon leaves the planet's gravitational domain; and (3)"disruption scenario", in which the moon migrates inwards until it reaches the Roche lobe and becomes disrupted by strong tidal forces. Applying the model to real cases from an exoplanet catalogue, we study the long-term stability of moon orbits around known exoplanets. We calculate the survival rate which is the fraction of the investigated cases when the moon survived around the planet for the full integration time (which is the age of the star, or if not known, then the age of the Sun).The most important factor determining the long term survival of an exomoon is the orbital period of the planet. For the majority of the close-in planets (<10 days orbital periods) there is no stable orbit for moons. Between 10 and 300 days we find a transition in survival rate from about zero to 70\%. Our results give a possible explanation to the lack of successful exomoon discoveries for close-in planets. Tidal instability causes moons to escape or being tidally disrupted around close-in planets which are mostly favoured by current detection techniques.


Introduction
Regular moons are predicted to form generally in planetary systems, as a direct outcome of planet formation. During the core-accretion, gravitational perturbations between planet embryos imply a series of constructive impacts up to the formation of a fully grown planet (Nagasawa et al., 2007). This phase may happen between a few Myrs to a few 100 Myrs after the star formation (see e.g. Chambers and Wetherill, 1998), and in these processes, satellites may form around the growing planets. Such giant impacts are advocated for the formation of the Earth's Moon (Canup, 2004) and for the formation of Uranus and Neptune's satellites (Morbidelli et al., 2012). Satellites may also be formed in rings, which could be natural outcomes of giant impacts or tidal disruptions, either in the planet formation phase or later in a relaxed planetary system (see e.g. Charnoz et al., 2009;Canup, 2012;Crida and Charnoz, 2012).
An alternative scenario invokes the accretion of moons in the gaseous circumplanetary envelopes that surround the most massive giant planets during their growth in the gaseous protoplanetary disks (see e.g. Canup and Ward, 2006;Sasaki et al., 2010;Miguel and Ida, 2016). Inside the circumplanetary disk, the solid material is replenished by the surrounding protoplanetary disk. This solid material is thought to accrete in the form of large satellites in a way similar to planets. Then the young satellites migrate inward, and sometimes sink into the planet. When the circumplanetary disk disappears, the satellite system remains and then undergo long-term tidal evolution.
After the planets have been formed, the orbits of moons evolve due to the planet's tides. Tidal forces lead to the dissipation of energy. On the one hand, this usually results in the gradual expansion of the moon's orbits, if the evolution starts outside the planet's synchronous radius. (A synchronous orbit is where the moon's orbital period is equal to the planet's rotational period.) On the other hand, moons inside the synchronous orbit will migrate inwards and can reach the planet's Roche limit where the moons are disrupted by tidal forces.
The timescale of orbital evolution depends mainly on the intensity of dissipation in the planet's interior, which relies on the material composition and internal structure of the planet. The dissipation coefficient also varies with the changing planet radius (for gas-rich planets Recent papers attempting to compute the quality factor, Q, (also known as tidal dissipation) show that Q could even be frequency dependant, and thus may strongly depend on the orbital distance (see e.g. Ogilvie and Lin, 2004). However, in the current stage of our knowledge, the large uncertainties in the origin of the dissipation process and in the internal structure of detected exoplanets preclude any theoretical calculation of Q for a specific planet. For this reason, in the present paper, a classical approach is adopted where Q will be considered as a constant, and many values will be explored. The dissipation scales as Q −1 and we can assume it to be in steady-state when exploring those scenarios which lead to moons on the long-term stability. For rocky planets, Q is low and the dissipation is high (Q ∼10-100), conversely for giant planets Q is high (∼ 10 4 -10 5 ) and dissipation is low (Goldreich and Soter, 1966;Alvarado-Montes et al., 2017).
The quest for exomoon discovery is a challenging task in itself, and upon finding an exomoon we could empirically constrain the Q value of the host planet by considering the planetary distance and the age of the star. This will provide invaluable information about the exoplanet interior, too, since high/low dissipation rates may depend on the planet core and interface with mantle Remus et al., 2012).
Recently, Alvarado-Montes et al. (2017) investigated the orbital changes of exomoons around close-in giant planets due to tidal effects. In addition to the tidal evolution, they implemented the changes of the giant planet's physical parameters into the model. The radius of the planet decreases over time, and this change is most prominent in the first billion year after the formation of the planet. The value of k 2 /Q slowly increases, which also has an impact on the orbital evolution of the moon.
Probably the most intriguing question concerning the moons is the lack of a firmly confirmed moon in exoplanetary systems. Before novel observation projects are designed for exomoon detection, first we need to find the answer to the question where the exomoons would be. Did we use suboptimal signal detection strategies, or did we simply look for moons in systems where their existence is not too likely? This latter possibility is the central question of this paper. Barnes and O'Brien (2002) showed that exomoons are more stable at longer planetary orbital periods, and also that the orbit of exomoons can tidally evolve first outward and then inward (towards the planet) as the planetary rotation slows due to tidal interaction with the star. Heller (2012) studied the orbital stability of habitable moons using the method of Kipping (2009) and it was found that for stars of masses below ≈ 0.2M ⊙ , moons are not likely to stay on stable orbits around planets which are in the habitable zone. Later Zollinger et al. (2017) showed that for planets in the habitable zones of 0.5M ⊙ stars, the orbit of moons are not stable as they migrate towards the planet.
Here we investigate the evolution of hypothetical satellites orbiting the currently known exoplanets with numerical integration, under the influence of tidal interaction. We map a wide range of parameters, including planet periods, moon to planet mass ratios, and k 2 /Q. Instead of studying a few cases in large details and with high accuracy, which would be meaning less due to our lack of knowledge on both the physical origin of tidal dissipation, and on the planet's internal structure, we are more interested in the statistical survival rate of moons. In order to investigate a large range of configurations over billions of years of evolution, we used the average description of Barnes and O'Brien (2002) which are formulated in Eqs. 1-3 (Section 2). We apply this method to identify different evolution pathways of moons (Section 3) and then we calculate survival rates of hypothetical moons around known exoplanets with the aim of pointing out to specific planets as observation targets which might have satellites (Section 4 and Appendix A).
Independently of our work, Sucerquia et al. (2019) used similar methods to investigate the orbit stability of exomoons around close-in planets, and estimated their detectability through TDV and TTV effects. Guimarães and Valio (2018) had a similar goal with the focus on the photometric detectability of hypothetical moons around Kepler planets. They selected 54 planets which could possibly host observable moons. Although we also list planets with survival rates for moons, we focus our investigations on the survival statistics.

Orbital evolution of satellites
We calculate the orbital evolution of the planet and the moon from the following equations (Barnes and O'Brien, 2002;Sasaki et al., 2012).
where n m and n p are the mean motions of the moon and the planet, k 2p is the second order Love number of the planet which relates to the dilatation due to tidal response, Q p is the quality factor of the planet which describes the dissipation of energy in the body, G is the gravitational constant, R p is the radius of the planet, M m , M p and M ⋆ are the masses of the moon, the planet and the host star, respectively, Ω p is the rotation frequency of the planet, and α is the prefactor of the moment of inertia (αM p R 2 ). Circular orbits are assumed for both the planets and the satellites. Since Equations (1) and (2) have the form ofẏ = C y p (where y is either n m or n p , p is the exponent of 16/3 while C contains all of the astrophysical constants) and Equation (3) can directly be integrated once the time evolution of n m and n p are known, this set of equations can be solved analytically on the domains until n m , n p or Ω p do not surpass each other; in this case, the presence of the sgn() function would not allow the fully analytical solution. We therefore implemented the solution in a piecewise-analytical manner, monitoring the sign changes in the differences between n m , n p or Ω p while in the distinct pieces, the exact soluion is evaluated. Thus, besides the vicinity of the points where the signs commute, the solution does not even depend on this stepsize, allowing us to highly speed-up the computations.
For calculating the survival rate of moons around known exoplanets, we let the orbits of the planet and the moon to evolve until the current age of the star, or if it is not know then until the age of the Sun (4.57 Gyr). We solve the equations analytically with fixed time steps to monitor whether the moon reaches the Roche lobe of the planet or the critical distance (half of the Hill sphere). In these cases the integration ends because the moon is considered to be lost either because of tidal disruption or escape from the planet's gravitational bond (for the latter case see Domingos et al., 2006). The critical distance is calculated by R crit = 0.5 R Hill = 0.5 a p M p /(3M ⋆ ) 1/3 . If at the end of the calculation the moon is still orbiting the planet, then it is counted as a "survival" case.

Parameters used for the integrations
Stellar mass (M ⋆ ), planet mass (M p ), semi-major axis (a p ) and radius of the planet (R p ) were taken from the The Extrasolar Planets Encyclopaedia (http://exoplanet.eu, data obtained on March 25, 2021). Only those planets are considered where the M p or the M p · sin i or the R p value is known (at least one of them). Planets with higher masses than 13 Jupiter mass are ignored as well as those with M ⋆ < 0.08M ⊙ . In addition, the selected planets also had to have either the orbital period (P p ) or the semi-major axis data. After this selection 4064 planets remained.
For both the radius and mass (or M p ·sin i if M p is not given) of the planet, the uncertainties were taken into account wherever possible. If there was no uncertainty in the data, then simply the given mass or radius was used. If there was a symmetric uncertainty given, then the mass or radius values were generated by choosing a random value with a Gauss distribution. The symmetry was determined with the same method as described in Chen and Kipping (2017, Section 2.2) where the difference between the upper and lower errors must be below 10%. For these symmetric cases the mean of these two values was taken as the standard deviation of the data. If the difference was higher than 10%, then the errors were ignored and simply the mean value was used. We applied exactly the same method for those planets where only the M p · sin i value was known, not the real mass. In these cases the error of the M p · sin i was taken into account wherever possible. This way we could generate realistic data for all planets, taken into account the uncertainties in the measurements.
There are 1007 cases (out of 4064) where both the radius and the mass (M p or M p ·sin i) of the planet are known. For all the other planets, the Forecaster model was used for generating realistic values for the missing parameter (Chen and Kipping, 2017). This tool predicts the mass (or radius) of the planet based on the radius (or mass). As an input, the generated (abovementioned Gauss distributed) mass (or radius) values were used as posteriors. This way, for each run the planet had a different mass or radius, or both (if the original data had symmetric errors).
See the following description of the defined ranges and distribution of the chosen random values for the rest of the (initial) parameters.
• Mass of the moon: The mass of the hypothetical moons is selected randomly with a uniform distribution between 1% and 10% of the host planet's mass, for each planet.
• Semi-major axis of the moon: The initial semi-major axis of the moon is chosen randomly with a uniform distribution on a logarithmic scale between 2 R p and R crit .
• Radius of the moon: Estimating the radius of exomoons is very challenging, because only the Solar System serves as an example. We choose to define it through their densities (ρ m ) using their mass which is already chosen. Based on the big moons (radius >200 km) in the Solar system, we define three categories depending on the locations with respect to the snowline (a snow ). This is because beyond the snowline there are more icy, hence less dense bodies in the Solar system. The snowline is calculated as a snow = T 2 eff ·R ⋆ /T 2 0 where T eff is the effective temperature of the star and T 0 ≈ 230 K is the equilibrium temperature at the planet's sub-stellar point (Cowan and Agol, 2011). A fourth category is also given for those cases when T eff or R ⋆ is not known, hence a snow cannot be calculated. The ρ m values are selected randomly with a Gauss distribution with the following parameters.
-If a snow is not known: the mean value of ρ m is 2.5 g/cm 3 , σ = 1/2 g/cm 3 .
• Rotation period of the planet: In the Solar system the planet that has the slowest spin around its axis is Jupiter with 9.9 hours. Scholz et al. (2018) defines an empirical spinmass relation for planets and brown dwarfs. According to their work, the spin velocity scales as v ∼ √ M, and from this relation it can be concluded that the spin period does not exceed 5 days even for brown dwarfs. (Spin-orbit resonances are not taken into account.) Based on these information and considering possible deviations, we set the spin period to a random value with uniform distribution between 10 hours and 5 days.
• Quality factor of the planet: For the quality factor (Q p ) three categories are definied.
-Rocky planets if R p < 2R ⊕ : 10 < Q p < 500 (Goldreich and Soter, 1966). Since lower values are more probable, the distribution of the randomly chosen Q p value is uniform on a logarithmic scale. -Ice/gas giants if R p ≥ 2R ⊕ and P p > 10 days: 10 3 < Q p < 10 6 (Goldreich and Soter, 1966;Ogilvie and Lin, 2004;Lainey et al., 2012Lainey et al., , 2017. For these planets intermediate values are the most probable in this range, hence Q p values are chosen with a log-normal distribution with a mean of Q p = 10 4.5 and σ = 10 0.5 . -Hot jupiters if R p ≥ 2R ⊕ and P p ≤ 10 days: Q p ∼ 5 · 10 6 (Ogilvie and Lin, 2004).
We do not have much information on the possible Q p value of these planets, hence we use a Gauss distribution with Q p = 5 · 10 6 as the mean value and σ = 2 · 10 6 . Q p is a critical parameter, as it controls the tidal evolution of the planet, however its value is very poorly constrained, and, in the current state of our knowledge it could vary by several order of magnitudes even for similar planets . For the giant planets in the Solar System, Q is mostly constrained by the secular acceleration of a big moon orbiting relatively close to the planet and also by theoretical considerations. For these planets Q is estimated to be in the range of 10 4 -10 5 (see e.g. Goldreich and Soter, 1966;Ogilvie and Lin, 2004;Lainey et al., 2009). However, if Q is not treated as a constant but as a function of the tidal frequency, then we obtain much lower values for Saturn with Q about 2500 (Lainey et al., , 2017. Q is mostly unknown for ice giants like Uranus and Neptune, but theoretical considerations put them in the range > 10 4 (Goldreich and Soter, 1966). However this estimate is based on some unconstrained assumption on the moon formation scenarios and ages. There is also a more recent estimation for Neptune which puts the Q value between 9000 and 36000 (Zhang and Hamilton, 2008).
Conversely, the Q of terrestrial planets is very low, in comparison to gas giants, and is typically between ten and a few hundreds (see e.g. MacDonald, 1964;Goldreich and Soter, 1966;Smith and Born, 1976;Ray et al., 2001;Rainey and Aharonson, 2006;Lainey et al., 2007). The Q for these planets is better constrained because of the possibility of laboratory tests of rocks (where possible). While Q might go up to 500 for rocky planets, it is typically below 200 in the Solar System.
It is highly challenging to draw conclusions about the possible Q value of exoplanets. It seems reasonable to assume that rocky planets have similar values as in the Solar System. We can also speculate that maybe the Q of giant planets is in the same order of magnitude as their Solar System counterparts, however, theoretically, with low-density ice-rich cores, Q in the range a few hundreds would also be possible for a Saturn-like planet (Lainey et al., 2017). In hot-jupiters, because of the assumed synchronous spin and circularized orbit, Q might be around 5 · 10 6 (Ogilvie and Lin, 2004).
The Love number k 2 varies modestly from one planet to another in the Solar System but is between 0.1 and 0.6 (Gavrilov and Zharkov, 1977). We choose k 2p = 0.5. However k 2 never appears alone, but only in k 2 /Q. So we fix the value of k 2 and vary the value of Q to save a free degree of liberty that would lengthen the computation time.
For the same reason, the parameter used for calculating the gyration radius is also fixed to α = 0.3. It is expected that planets with different compositions and structures have different gyration radii, but by keeping α fixed, a free parameter is removed.
For each planet ten thousand runs are made to get a wide sample of the different initial parameters. What we call the "survival rate" is simply the ratio of surviving moons (those that stayed in orbit around the planet until the end of the integration) and the total number of tested configurations.
For a few example cases, the parameters of the planet and the moon are set to specific values which are described in details in Section 3.

Example cases
In the example cases (see Figs. 1 and 2), the system configuration consists of a Solar-mass star with an Earth-like planet and a moon which has different mass in the different cases. The quality factor of the planet is Q p = 10 in all cases. The initial parameters used for the different runs are summarized in Table 1. Note that the time shown in the horizontal axis is significantly shorter in Fig. 1 (where it ends at 5 · 10 6 years) than in Fig. 2 (where it ends at 10 10 years).

Evolution outcomes
In the integrations, the moon has been observed to follow three different families of evolutionary tracks, and depending on the actual parameter set, all of these three tracks can evolve rapidly, resulting in an unstable system, or the evolution can be slow enough to lead to a surviving moon. Both the morphology and the stability depends on the position of the Hill sphere, thus eventually, on the masses and orbital periods of both the planet and the moon, the stellar mass, and Q p .
3.1.1. Disruption scenario When the planet is close enough to the star, the tidal torque on the planet from the star very effectively spins down the planet. The tidal bulge is shifted  Fig. 2  slightly toward the trailing side of the moon's orbit, the moon quickly looses energy, spirals inward and is disrupted at the Roche lobe. The time scale of this migration is very short (in the order of 1-200 million years, mostly depending on the stellar distance and on the moonplanet mass ratio) compared to the 10 Gyr integration time, therefore we consider that these moons are unstable. This scenario is described in details in Alvarado-Montes et al. (2017), and our results are completely compatible with their findings (Fig. 1, top panel).

Escape scenario
If the planet's orbit is close to the star, then the Hill radius is smaller. On any orbit inside the Hill sphere, the moon suffers strong tidal forces, especially if the planet is relatively massive, which makes the outward migration very fast. As a consequence, if the moon formed close to the critical distance from the planet (half of the Hill sphere), then it will quickly escape, before synchronization could happen ( Fig. 1, lower panel).

Tidal locking
In connection with tidal locking, we define three different scenarios.
Moon-planet locking: the Selenoid scenario. The evolution track is somewhat similar to the Earth-Moon case. First, the moon migrates outwards, until a spin-orbit synchronisation is reached between the orbital period of the moon and the rotation of the planet. In the later evolution, the moon slowly turns back, the typical time scale before colliding with the planet was in the order of 10 billion years in our integrations. (More specifically, the end state of the moon is reaching the Roche radius of the planet, where the moon will be disrupted). Therefore, we consider this scenario as stable (Fig. 2, top panel). We note that the similarity with the Earth-Moon case is limited to the synchronization between the two bodies, but there is a qualitative difference: here the planet spin is synchronizing with the orbit of the moon, while in the Earth-Moon case the spin of the Moon is synchronized with the orbit of the Earth.
Star-planet locking. This scenario is similar to the Selenoid, but rather than synchronising with the moon, the planet is synchronised to the central star. We did find this possibility in several test runs, and seems to be typical for certain configurations (closein systems with lightweight moon). In these cases, the moon's migration is slow and has a similar time scale than in the Selenoid case. The path is qualitatively identical: first the moon migrates outwards and then later inward (Fig. 2, middle panel). This configuration is also stable on the long term, and the moon can be observed (with difficulties coming from its low  Table 1. Blue dashed horizontal line at 6 · 10 −6 s −1 represents the location of the critical distance from the planet (half of the Hill radius of the planet) beyond which the moon escapes. Red dashed horizontal line shows the Roche radius of the planet expressed in the unit of mean motion. Upper panel (disruption scenario): the moon quickly synchronises with the planet, but when migrating inwards to the planet, it reaches the Roche lobe and becomes tidally disrupted. Lower panel (escape scenario): the moon slowly migrates outwards meanwhile the planet spin evolves. When the moon reaches the critical distance (after about 4 million years), it escapes from the planet. mass and small size). Temporarily locked scenario. We found another tidal evolution scenario with rather complex evolution for which we have not found earlier discussion in the literature. At the beginning it has essentially the same evolution as the star-planet synchronisation, but since the moon migrates inwards, here the tidal forces from the moon is increasing, and from a point they majorate over the tidal forces from the star. At this point, the planet rotation will be decoupled from the mean motion, and resynchronizes with the moon. The inward migration of the moon is enhanced from that point, and the moon soon reaches the Roche limit of the planet. In the scenario shown at the bottom panel of Fig. 2, the planet's synchronisation with the star happens after approximately 2 · 10 9 years, where there is the turning point of the moon's orbit, too, and the break-up of the star-planet synchronisation occurs at ∼ 7 · 10 9 years.  Table 1. Blue and red dashed horizontal lines represent the location of the critical distance for moon escape and the location of the Roche limit, respectively, expressed in the unit of mean motion. Upper panel (moon-planet locking scenario): the moon synchronises with the planet after ∼ 3 · 10 9 years; up to this point, the moon slowly migrates outwards from the planet, but after synchronisation the moon's orbital period slowly decreases. Middle panel (star-planet locking scenario): for smaller M m /M p mass ratio the planet synchronises with the star and becomes tidally locked. Bottom panel (temporarily locked scenario): first the planet synchronises with the star (at around 2·10 9 years) but because of the moon's inward migration, later it synchronises with the moon instead (∼ 7 · 10 9 years), but after that, the moon shortly reaches the Roche limit of the planet.
Because of the large time scales, this configuration is also stable on the long term, and moons in such systems could be found.

Time scales
The time scales of the discussed scenarios depend on the migration rate, hence indirectly on the distance of the Hill-sphere which can limit the evolution path very close to the star. Close-in planets always evolve very rapidly, because (i) the Hill-sphere is close to the planet, and the moon reaches the critical distance from the planet very soon, or turns back at a very close distance to the planet; and (ii) since the moon is always very close to the planet, very prominent tidal forces emerge and the orbit evolution will be very fast. For most close-in configurations, similarly to the hottest sub-Jupiters known, the tidal escape time scale of the moon can be even less than 1000 years. The only possible scenarios for those cases when the planet is too close to the star are the disruption and the escape scenarios ( Fig. 1 upper and lower panels).
If the planet orbits the star at a larger distance, the path to the the critical distance from the planet is longer and the evolution is slower, hence the time scale of stability can reach billions of years, or even the Hubble time scale (these are the tidal locking scenarios, see Fig. 2). For selenoids, this means that the fallback takes several billion years, and for runaway moons, the runaway time takes several billion years. In other words, with increasing planet orbital period first we get quickly falling-back selenoids, then stable selenoids ( Fig. 1 and 2 upper panels).

Known exoplanets
We calculated the evolution of possible model moons in the known planetary systems with the aim of pre-selecting those exoplanets which can physically have a satellite. These planets might worth a detailed follow-up observation for a moon survey. Planets and their model moons were simulated according to the recipe described in Sect. 2. The results are listed in Table A1 (in Appendix A), ordered by decreasing survival rate. (Only those cases are shown in the table which have non-zero survival rates.) The table also contains the mean of the integration times for each planet. If the survival rate is a hundred percent, then the mean of the integration times will be equal to the age of the star (or of the Sun, if the star's age is not known, as described in Section 2). If in some cases the integration stops earlier, then the mean of the integration times will be shorter. This parameter together with the survival rate can be considered as a measure for the stability of moon orbits. Fig. 3 shows the survival rate for test exomoons as function of the orbital period of the planet. Because the critical escaping distance is small below 10 days orbital period, only a few exomoons stay in a stable orbit around the planet. Above ∼ 300 days, the survival rate is about 70-90 percent, as a consequence of the large critical distance. Between about 5 and 300 days we find a "transition zone" where the increasing size of the critical escaping distance results in higher and higher survival rates. Survival rate [%] planet with measured radius and mass planet with estimated radius planet with estimated mass Figure 3. Survival rate of moons around known exoplanets as function of the planet orbital period. Black coloured dots: planets with known radius and mass data. Teal coloured dots: planets without known radius data, for these, their radius was estimated based on their measured mass (or M p sin i). Yellow coloured dots: Planets without known mass data, for these the masses were estimated from their measured radius.
The saturation at ∼85% survival rate for large orbital periods is caused by the chosen initial semi-major axes of the moons which in some cases is below the synchronous orbit. In these cases the moons spiral inwards to the planet until they reach the Roche sphere.
Only a fraction of known exoplanets are promising exomoon hosts according to this work. This is primarily because the orbital period behaves as a moon "censor": not much satellites are expected (survival rate mostly below 40%) for planets in shorter period orbits than 10 days. This is not surprising, since planets this close to the star have small Hill radii, hence stable orbits for a moon are limited to a very narrow space. As the orbit of the moon evolves, sooner or later it will reach the instability limit and escape from the planet's gravitational bond.
Interestingly, the planets with known mass and radius data (see the black dots in Fig. 3) are showing smaller deviation above ∼100 days orbital period, compared to those for which the radii or masses are estimated from their other parameter (teal and yellow coloured dots, respectively). For the planets with known mass and radius value, the survival rate reaches 50% at around 100 day orbital period of the planet, while for the other planets it is much more diffused in between approximately a hundred and a thousand days. Choosing a mass or radius value for the planet is an extra free parameter which naturally causes higher uncertainty in the results. This might be the reason of the shift of the transition zone from 5-100 days to 30-300 days orbital period for the two kinds of planets. For those where we have to add random mass or radius values, we will find more unstable cases, hence the survival rate starts increasing only at higher orbital periods where the critical escape distance from the planet is larger. Another explanation for the shift could be an observation bias, since the planets for which both mass and radius data are measured are those for which both transit and radial velocity observations were made. This means that they are orbiting close to the star and in the meantime they are also massive. This might also explain why the transition zone is closer to the star for these planets.
For known exoplanets with long orbital periods, the host stars are quite faint, and the S/N ratio required for a successful moon detection is not achievable by even as precise instruments as the Kepler telescope. This finding is the resolution of the seeming paradox of missing exomoons in Kepler data. Moons are more likely to have a stable orbit around long orbital period planets, but these cases are difficult to observe. This result is also an alert for later transit finding surveys to focus on planets around bright stars with 10-100 days orbital period.

Discussion
The successful detection of a moon around an exoplanet relies basically on two independent factors: the moons must exist (occurrence), and those moons must be observable. The occurrence is a product of two factors: the formation rate and the stability time scale. Formation studies of moons (and planets or stars) are beyond the scope of this paper, however, with the described formalism, we estimate the survival rate of moons around currently known planets.
The survival rate of model moons in known exoplanetary systems shows that satellites can only stay in orbit around the host planet if the orbital period of the planet is larger than five days. In our investigations we found a transition zone between 5 and 300 days orbital period where the survival rate is growing with larger orbital periods. Above 300 days orbital period the moon survival rate reaches 70-90% for the investigated planets. These specific numbers however depend on the range of chosen initial parameters. Also, increasing the integration time could result in decreasing the survival rates in some cases. This would be only a quantitative change as it would not alter the relation between the survival rate of moons and the orbital period of the planets. Our investigation serves statistical purposes, and it shows a clear trend in survival rates that increase with planetary orbital period.
This result is in good agreement with the findings of Alvarado-Montes et al. (2017), who developed a more detailed model including the changes of the planet's physical parameters over time as an additional factor influencing the orbital evolution of satellites. Their results showed that moons of planets with lower than ∼ 60-day orbital periods collide with the planet, but for planets with larger stellar distances, the satellites slowly migrate outwards and stay inside the Hill-sphere of the planet during the whole 4 Gyrs integration time ('realistic' case, Fig. 9 in their article). Including the change in the planet radius and k 2 /Q tidal parameter over time in their model increased the chance of moon survival. This implies that neglecting these factors in our study might underestimate the moon survival rates.
Our results are in line with those of Sucerquia et al. (2019), too, who concluded that long-period planets are more favourable moon hosts than short-period ones. They used a very similar calculation method with the aim of estimating the TTV and TDV signals for detection. They conclude that in their investigated phase space no system produces TDV signals strong enough to be observed, but there are a very few cases where the TTV signal would be observable with TESS and Kepler. Larger moon-to-planet ratios and stellar distances could help observations. With a similar aim, but using a different method, Guimarães and Valio (2018) also investigated known planets which could host moons. Beside checking the possible semimajoraxes of moons, they also took into account their detectability, providing a list of 54 Kepler planets and candidates which could have detectable moons.
Recently, Tokadjian and Piro (2020) made a similar work to ours in which they investigated the orbit evolution of moons around known exoplanets using simplified tidal lag models. They identified 36 habitable zone planets around which a moon might be stable for long time scales (at least for 1 Gyr). We found that from these planets 29 have at least 50% survival rate for exomoons (typically between 50 and 65%.). There are also four planets with survival rates higher than 65%: these are Kepler-459 b (69%), Kepler-456 b (72%), Kepler-1654 b (70%) and Kepler-1647 b (69%). They also found that for 21 planets the time scale before loosing the moon is larger than the Hubble time in all investigated cases. Our results show moon survival rates typically between 55 and 70% for these planets with mean integration times (t mean ) between 2 and 3 Gyrs (see Table A1). There is one outstanding case (Kepler-1628 b) which resulted in only 9% survival rate in our study. The explanation of our lower stability levels in general is that a much bigger phase space was explored for the initial parameters compared to the study by Tokadjian and Piro (2020). Also, since the age of these planets is typically lower than 5 Gyrs (and in some cases it is not known), our calculation stopped much earlier.
Martínez-Rodríguez et al. (2019) investigated the orbital stability and habitability of exomoons around known planets of M dwarfs. They found 4 planets which were in the habitable zone and for which the orbit of the moon is stable for a long time scale (the migration time of the moon is larger than the Hubble time). These four planets are HIP 12961 b, HIP 57050 b, GJ 876 b and c (in their paper under the names of CD-23 1056 b, Ross 1003 b, IL Aqr b and c, respectively). Interestingly, for these planets our model gave 0%, 0%, 6% and 0% survival rates, respectively. One of these four planets, however, HIP 57050 b was thoroughly investigated also by Trifonov et al. (2020) (planet name: GJ 1148 b in their paper) where they were using secular theory and direct N-body integrations to study the possibility of a moon. Similarly to our results, they have found that it is unlikely that a moon would stay in stable orbit around the planet.
Zollinger et al. (2017) investigated the tidal evolution and habitability of exomoons assuming possible configurations, and found that moons in the habitable zone of M dwarfs experience a tidal heating rate comparable to the tidal heat flux of Jupiter's moon, Io (∼ 2W/m 2 ). For stars with masses 0.2M ⊙ , because of the strong gravitational effect of the nearby star, the tidal heating rate of a Mars-like moon around a Jupiter-like planet is much stronger than that of Io's, implying that these moons might not be habitable. The maximum limit of tidal heating for habitability, however, is not known, and atmospheric effects can influence the habitability of bodies, as well (Barnes et al., 2013).
Recent hydrodynamical simulations showed that icy moons can form from circumplanetary disks of ice giants such as Uranus and Neptune, indicating that there might be a large population of moons around exoplanets beyond the snowline (Szulágyi et al., 2018). According to our findings, these moons, if formed, can stay in stable orbits for a long period of time.
Older stars host more evolved systems, in which there is a higher possibility that most satellites have disappeared due to either disruption (at the Roche radius of the planet) or escape. The presence of a massive exomoon may prevent synchronization of the planet with the star and promotes synchronization with the moon. In this case, satellites also evolve into a 1:1 spin-orbit resonance, getting to a mutually synchronised planet-moon system. Less dissipative planets (gas giants and ice giants) are more prone to keep satellites because of slower orbital migration (see e.g. Sasaki et al., 2012). If an exomoon is observed around a planet, it can provide a strong constraint on the value of the quality factor, Q, which describes the dissipation in the planet, and gives hint on its interior.
Since younger stars host less evolved systems, there is a higher chance for planets to still have satellites (which will escape or be destroyed later). But satellites on less distant orbits are more difficult to observe (Szabó et al., 2006;Simon et al., 2012). If more than one moon is found, than their orbital architecture may give information on their formation, and provide a strong constrain on Q. It may also be possible to infer the presence of rings.

Conclusion
There is no confirmed exomoon detection to date, but the quest for finding moons in Kepler data is still in progress. Teachey et al. (2018) reported an exomoon candidate, Kepler-1625 b I (see also Teachey and Kipping, 2018), but this discovery is under debate (Kreidberg et al., 2019;Teachey and Kipping, 2019). Regardless of the existence of this moon, we found stable orbits for satellites around this planet. Recently, the presence of an exomoon candidate was suggested based on observations of sodium and potassium lines around planet WASP-49 b (Oza et al., 2019). Further observations are needed for confirmation.
As a result of our calculations, we conclude that planets with long orbital periods support high survival rates for moons, and close-in planets, which are easier to observe, are less likely to host moons. This can serve as an explanation for not having a confirmed exomoon discovery so far.
Thanks to high accuracy photometry and long-term follow-up missions, CHEOPS is measuring the radius of numerous known exoplanets through transits with unprecedented accuracy. This paves the way for characterizing the close environment of exoplanets, and can unveil the presence of exomoons and exorings. Moons and rings are natural companions of planets. Neptune sized planets are especially promising targets owing to their (putative) lower dissipation and large mass compared to Earth-like planets. Such targets will be in the reach of the CHEOPS mission. Following the CHEOPS mission, ESA's PLATO 2.0 mission, capable of discovering and following exoplanet systems and host stars, and the ARIEL mission will continue the task of collecting precise exoplanet light curves. Whereas exomoons have never been firmly detected up to date, they should exist and their successful detections will provide invaluable information on the planet's interior and formation process (Crida and Charnoz, 2012). Detecting the first exomoons will be a major discovery as moons could provide habitable environments, in particular if they orbit volatile rich planets (like giant planets) and are located inside, or close to, the habitable zone (Forgan and Dobos, 2016;Dobos et al., 2017).  Table A1 for a list of exoplanets in a decreasing order of moon survival rate (see Section 4 for details). A horizontal line at 50% moon survival rate separates the cases where moons are more probable to stay in a stable orbit around the planet. Keep in mind that, as discussed in Section 5, the calculations served only statistical purposes, and the actual survival rates listed below might change for different initial conditions and different restrictions.  Table A1: List of planets ordered by decreasing moon survival rate. Only those cases are listed where the survival rate (surv. rate) is non-zero. Note also the mean of the integration time (t mean ) which is the mean value of the lengths of the simulations. In some cases the age of the system is believed to be very short which results in high survival rates but where t mean is very short.