Resonant sub-Neptunes are puffier

A systematic, population-level discrepancy exists between the densities of exoplanets whose masses have been measured with transit timing variations (TTVs) versus those measured with radial velocities (RVs). Since the TTV planets are predominantly nearly resonant, it is still unclear whether the discrepancy is attributed to detection biases or to astrophysical di ff erences between the nearly resonant and non resonant planet populations. We defined a controlled, unbiased sample of 36 sub-Neptunes characterised by Ke-pler , TESS, HARPS


Introduction
Planets with radius in the 1-4R Earth range are estimated to exist within a hundred days of orbital period around 30 − 50% of all Sun-like stars (Lovis et al. 2009;Howard et al. 2012;Fressin et al. 2013).In order to understand the nature of these objects, it is crucial to constrain both their masses and radii, and, thus, their densities.The bulk of exoplanet discoveries is done by transit surveys such as Kepler/K2 and TESS, which measure the planets' radii.Then, the mass is typically estimated using the radial velocity method (RVs).For compact multi-planetary systems, the mass can also be estimated using transit timing variations (TTVs).In particular, when the period ratio of two planets is close to commensurability, namely, P out /P in ≈ (k + q)/k, the planets can exhibit TTVs due to their proximity to the mean motion resonance (MMR).Despite their relative rarity (Fabrycky et al. 2014), (nearly) resonant systems are over-represented in the population of planets with both mass and radius measurements because the four-year baseline of the Kepler mission allowed estimations of their masses through TTVs at no additional cost.
Over the last decade, numerous studies (e.g.Wu & Lithwick 2013;Weiss & Marcy 2014;Steffen 2016;Mills & Mazeh 2017;Hadden & Lithwick 2017;Cubillos et al. 2017;Millholland 2019;Leleu et al. 2023;Adibekyan et al. 2024) have noted (and discussed), the apparent discrepancy in density between the planets characterised by TTVs and radial velocities RVs.However, the origin of this discrepancy remains unclear: it could be due to sensitivity biases inherent to each method, with photom-etry biased towards larger planets and radial velocity biased towards more massive planets.Recent results also showed that part of the TTV-characterised population had underestimated densities due to the difficulty of extracting transit timings for lowsignal-to-noise ratio (low-S/N) transits (Leleu et al. 2023).Hadden & Lithwick (2017) put forward a selection bias as possible explanation, since TTVs tend to allow the characterisation of small planets on larger orbital periods (hence, cooler orbits) than the bulk of the RV characterisation.It has also been proposed that the systems characterised by RVs and TTVs formed in different environments, such as with different disk metallicity (Adibekyan et al. 2024).However, the differences in physical properties could be due to the orbital configuration in which the planets are embedded (e.g.Weiss & Marcy 2014;Mills & Mazeh 2017;Goyal et al. 2023), since TTVs mainly characterise sub-Neptunes that are near mean motion resonances (MMRs), while the RV-characterised planets are more representative of the bulk of known exoplanets.In this paper, we explore the possibility that there is an intrinsic connection between the densities of sub-Neptunes and their resonant orbital configurations.with their uncertainties.The host properties are required in order to check for possible correlations between these parameters and the densities of the planets.Restricting this population to closein systems with periods in the 5-60 days range and radii between 2 and 4 R Earth has reduced this number to 133.The lower limit of the period range is chosen to avoid the lower part of the Neptunian desert, which is partly shaped by photoevaporation (Owen & Lai 2018).We define the (nearly) resonant population (in orange) as planets whose period ratio with an inner or outer planet satisfies ∆ MMR <0.05, where for q = 1 and k ∈ [1, 2, 3, 4, 5] or q = 2 and k ∈ [3, 5].While the non-resonant population (in blue) is defined ∆ MMR >0.05.This limit is set by the edge of the clump of nearly-resonant system found in Kepler (Fabrycky et al. 2014).The (nearly) resonant population, in orange, appears to be composed of lower-density planets than the non-resonant population.However, these populations could be affected by numerous biases.Notably, the (nearly-)resonant population is mainly characterised by TTVs, while the non-resonant population is mainly characterised by RVs.
Regarding the RV-characterised masses, a possible bias could come through the follow-up observation process; some planets could be dropped out after few RV points were taken, if the RV signature of the planets did not seem large enough.This selection process would lead to a bias in the literature towards higher masses and higher densities for RV-characterised planets.Using all sub-Neptunes in the 5-60 days range that were followed-up by HARPS or ESPRESSO, we show in Appendix A that such bias is absent for planets whose radius is above 2.7R Earth , since 93% (26 out of 28) of the followed-up planets in that radius range have published masses.These 20 RVcharacterised planets are the first part of our controlled sample.
For TTV-characterised planets, the mass-radius relation can strongly be affected by mass-eccentricity degeneracies Lithwick et al. (2012) and the manner in which the TTVs are extracted from the light curves (Leleu et al. 2023).To address the first point, we only use planets whose mass estimations have been shown to be robust against mass and eccentricity degeneracy (e.g.Hadden & Lithwick 2017;Leleu et al. 2023).Regarding the second point, Leleu et al. (2023) showed that TTVs are typically correctly estimated by usual methods if the signal to noise ratio of individual transits (S/N i ) is high enough.As shown in Fig.
A.2, we checked that the planets with radius above 2.7 had an S/N i > 3.5, which ensures that their individual transit timing can robustly be recovered (Leleu et al. 2023).Here also, it translates to the selection of stars that are bright enough and not overly active.If individual transits of planets can reliably be observed, there is no reason why we could not characterise denser planets, as their TTV signals would either be larger (Lithwick et al. 2012) or faster (Nesvorný & Vokrouhlický 2016), depending on whether the pair is near or inside a MMR.
We therefore define our controlled sample as the planets in the 2.7 to 4 R Earth range in the mass-radius diagram (shown in Fig. 2).This sample, detailed in  ity2 to be drawn from the same underlying population.However, the probability that their masses have been drawn from the same underlying population is p value =0.002 +0.010 −0.001 .We then checked whether this discrepancy could be attributed to different stellar metallicities, equilibrium temperatures (e.g.Hadden & Lithwick 2017), or stellar effective temperatures.With respect to all of these quantities, the two populations are similar, with p value of 0.62 +0.06 −0.13 , 0.25 +0.31 −0.18 , and 0.62 +0.12 −0.18 , respectively.Exploring possible 2D relations between these parameters, we performed 2D Kolmogorov-Smirnov tests (Peacock 1983).These results are shown in Appendix A. Across all our tests, the p value involving the mass are lower than the rest by two orders of magnitudes.We therefore conclude that the proximity to MMR is the main factor in the discrepancy between the mass, hence, the density, of the two sub-populations.

Full sample
We go on to consider the full sample shown in Fig. 1, taking the robust masses from Hadden & Lithwick (2017) and Leleu et al. (2023) when available.For this sample, the probability that the mass of the (nearly) resonant and non-resonant population is drawn form the same underlying population drops to p value =2.3e − 06 +2.7e−05 −2.1e−06 , while the rest of the explored parameters are consistent between the two sub-population.The full 2D comparison is in Appendix B.1.The 2D comparison was also performed for the full sample restricted to the 2.7-4R Earth (see Fig.

B.2).
To compare the relative densities of planets with different radii, we defined ρ rel = ρ/ρ re f , which is the ratio between the density measured for a planet and a reference density for a planet of the same radius, using the M-R relation from Parc et al. (2024) (black line in Fig. 1).To further illustrate the effect of the proximity to MMR on the density of planets, in Fig. 3 we show ρ rel as a function of ∆ MMR .Formally resonant planets (including resonant chains) are shown on the left, while single planets are shown on the right 3 .In the figure, we estimated the local mean and scatter of ρ rel by fitting a Gaussian distribution in a box sliding over log 10 (∆ MMR ) with a width of 1.Most planets for which ∆ MMR < 0.05 have a relative density below 1, and this is even more so for the formally resonant systems, which have a mean relative density estimated at 0.63 ± 0.05.On the contrary, planets for which ∆ MMR > 0.1 are more uniformly spread around ρ rel = 1 and single planets are on average denser, with a mean relative density of 1.38 ± 0.10.The right panel shows the envelope of the 0.16 − 0.84 quantiles of the Gaussian distributions of the relative density as a function of the mass measurement method used and the distance to the resonance.Only the nearly resonant population (0.001 < ∆ MMR < 0.05) had enough measurements to be estimated by both methods.From this analysis, we have drawn three observations.First, the correlation between the distance to MMR and the relative density of sub-Neptunes is visible when using RV-characterised planets alone (grey, green, and dark green Gaussians).Second, RV-and TTV-characterised planets have similar relative density where they overlap in the nearly-resonant population (dark green and dark purple Gaussians).Third, the distribution of relative densities increase both in mean value and in spread away from the resonance, with resonant planets having similar relative densities, while single planets have a larger spread.

Discussion
If the difference in density between (nearly) resonant and nonresonant planets is not due to observational or selection biases, nor to the type of star that the planets orbit, it must be a result of different formation and evolution pathways.For example, the (nearly) resonant planets could be puffier as a result of atmo-spheric inflation due to tidal heating (Millholland 2019;Millholland et al. 2020).This would arise if the (nearly) resonant have systematically larger eccentricities or planetary obliquities.This is supported by previous studies, which suggested that planets captured in mean-motion resonances might have their spinaxes tilted as a result of secular spin-orbit resonance capture during the orbital migration process (Millholland 2019;Millholland et al. 2024).
Another hypothesis is that resonant and non-resonant planets were formed in different locations with diverse formation conditions.For example, Lee & Chiang (2016) posited that most super-Earths and sub-Neptunes formed in situ in gas-poor discs towards the end of the disc lifetime.However, they suggested that the rarer class of low-density "super-puff" planets formed outside ∼ 1 AU and accreted thicker gaseous atmospheres due to more efficient cooling.Because they formed further out, they would have entered MMRs upon inward migration.
Alternatively, lower densities for planets in MMRs is expected by the model known as 'breaking the chains' (see section 4 of Bean et al. 2021).In that model, close-in systems of sub-Neptunes form in resonant chains due to the migration of planets in the protoplanetary discs and the positive torque at the inner edge of the disc (e.g.Goldreich & Tremaine 1979;Weidenschilling & Davis 1985;Masset et al. 2006;Terquem & Papaloizou 2007).Resonant chains can then become dynamically unstable after the gaseous disc dissipates (Terquem & Papaloizou 2007;Ogihara & Ida 2009;Cossou et al. 2014).This model reproduces the observed period ratio and multiplicity distribution of the close-in sub-Neptune population if ≈ 95% of the chains become unstable after the disc dispersal (Izidoro et al. 2017).Some of these instabilities lead to giant impacts which can eject part, or all, the primordial H/He atmospheres of the planets (Biersteker & Schlichting 2019a), resulting in nonresonant planets that are on average denser than their resonant counterpart.
This scenario leads to a second observable: known resonant chains tend to be remarkably coplanar (e.g.Agol et al. 2020; Leleu et al. 2021a), while planet-planet scattering is expected to increase the mutual inclination between planets.In Fig. 4 we show the minimum mutual inclination (i.e.assuming that the ascending nodes of all the planets are aligned in the sky) of all pairs of planets that were confirmed in the exoplanet archive.Here, we can also see that pairs further away from MMRs have a larger scatter in the mutual inclination, in agreement with the 'breaking the chains' model.We note that the actual trend might be stronger, since here we only measure the minimal mutual inclination, the impact parameter of (nearly) resonant systems could be biased by unaccounted TTVs (García-Melendo & López-Morales 2011) and misaligned, spread-out systems have a lower transit probability.
The 'breaking the chains' mechanism (e.g.Izidoro et al. 2017Izidoro et al. , 2022) ) is also observed in different planet formation models (e.g.NGPPS Emsenhuber et al. 2021a;Burn et al. 2024).In Fig. 5, we show synthetic systems from Izidoro et al. (2021Izidoro et al. ( , 2022) ) and NGPPS (Emsenhuber et al. 2021a;Burn et al. 2024).These populations (described in Appendix C.1 and C.2, respectively) simulate the formation of planetary systems from planet embryos in the proto-planetary discs up to ∼50 millions of years after the disc dispersal.As can be seen in the middle panels, both populations harbour lower-density planets for the (nearly) resonant sub-population, while the non-resonant populations have larger and more diverse relative densities.The right panels show that the mutual inclination between each planet pairs tend to be larger for larger distance to MMRs.In both populations, these features are linked with post-disc instabilities and collision, shown with star markers in Fig. 5.These instabilities can be due to the configuration of the chain itself, but also to the ex-istence of an outer more massive planet (see Appendix C.2 and also Schlecker et al. 2021;Izidoro et al. 2022).In addition, for the NGPPS population, the larger density of the non-resonant population is partially due to the accretion of more rocky embryos from the inner system during the giant impact stage after gas disk dissipation.

Conclusion
Our results support the idea that the apparent discrepancy between TTV-and RV-characterised planets is astrophysical and due to different formation and evolution pathways of the characterised populations, rather than a method-related bias.The significance of the controlled sample (p value of 0.002 +0.010 −0.001 ) can be improved by a homogeneous analysis of the RV-characterised population, and a larger completion in the publication of mass measurement for smaller planets.For the TTVs, the population needs to be systematically checked for the mass or eccentricity degeneracies (Hadden & Lithwick 2017) and the robustness of the TTV extraction or the photo-dynamical analysis (Leleu et al. 2023), as well as ensuring that small dense resonant planets are not missed due to large TTVs (Leleu et al. 2021b(Leleu et al. , 2022)).In this study, we were also able to show that TTV and RV characterised planets had a similar relative density distribution for the nearly resonant population.A next step would be to get enough RV characterised resonant systems of sub-Neptunes, such as resonant chains, to check whether these results hold for that population.Finally, PLATO will enable the discovery and characterisation of systems by both TTVs and RVs, which should help to further alleviate the biases inherent to each method.In Izidoro et al. (2021), planetary seeds grow through pebble accretion and mutual collisions.Pebbles beyond the snowline are presumed to contain 50% water ice mass.As these pebbles migrate inward and cross the water snowline, they sublimate, losing their water component and releasing silicate grains.Collisions are modelled as perfect merging events conserving mass and momentum.
The planet formation simulations of Izidoro et al. (2021) provide planetary mass and composition but not planet size or radius.To compare the simulations outcome with observational trends involving planet sizes such as the exoplanet radius valley and the peas-in-a-pod trend, Izidoro et al. (2022) employ mass-radius relationships (Zeng et al. 2019) to convert mass into planetary radius.In addition, in this model, giant impacts (with projectile-to-target mass ratios greater than 0.1) occurring after gas disk dispersal strip primordial atmospheres, leaving behind either bare rocky or water-rich cores (Biersteker & Schlichting 2019b).According to Izidoro et al. (2022), approximately 80-90% of late impacts qualify as giant impacts in their model.Their model does not account for the formation of secondary/outgassed atmospheres but for only for primordial atmospheres accreted during the gaseous disk phase (H/He rich).
The stability of a primordial atmosphere for a planet that did not experience giant impacts after gas dispersal is estimated us- ing an energy-limited escape prescription, considering stellar Xray and ultraviolet radiation (Owen & Wu 2017).The criterion by Misener & Schlichting (2021) compares atmospheric binding energy to the energy received by the planet from 100 million to 1 billion years.An energy ratio smaller than unit indicates sufficient energy for atmosphere photo-evaporation.Planet sizes are computed following different planet models, which account, or not, for the presence of a primordial atmosphere following the atmospheric instability criterion of Misener & Schlichting (2021).

C.2. New generation planetary population synthesis (NGPPS)
The second set of theoretical planetary system calculations was obtained from the planetary population synthesis exercise conducted by Emsenhuber et al. (2021b) and updated by Burn et al. (2024).This set of simulations use the Bern model of global planet formation and evolution (Emsenhuber et al. 2021a).The starting point of the simulations are protoplanetary dust and solid disks around Solar-type stars with observation-informed distributions of initial conditions (Tychoniec et al. 2018, see also the discussion in Emsenhuber et al. 2023).The gas disk viscous surface density evolution equation (Pringle 1981) is solved assum-ing an α-viscosity of 2 × 10 −3 , photoevaporative mass loss, and a consistent temperature structure using opacities from Bell & Lin (1994), viscous and irradiation heating.Apart from 1% dust used as opacity source, the solid disk is assumed to be present in the form of planetesimals that are modelled as fluid and accreted by 100 growing seed embryos randomly distributed in the disk at initialisation.The composition of gas, embryos and planetesimals, initialised following Marboeuf et al. (2014), is tracked during planet growth.Gas accretion onto the planets proceeds by cooling and contraction of the previously accreted gas which is modelled in one dimension assuming hydrostatic equilibrium and energy release at the boundary between the solid core and the gaseous envelope.The simulations explicitly model gravitational interactions between the embryos using the mercury code (Chambers 1999) where an additional force is added due to the interaction of the planets with the gaseous disk.The force is derived from prescriptions of type I (Paardekooper et al. 2010(Paardekooper et al. , 2011) ) and II (Dittkrist et al. 2014) migration timescales, as well as eccentricity and inclination damping (Coleman & Nelson 2014).This will commonly lead to an inwards movement at moderate excitation of the planets during the gas disk stage.If two embryos collide, we assume a perfect merging of the solid core (including volatile species) and immediate ejection of the  2023), as function of the S/N of individual transits (S /N i ).The vertical grey dashed line shows the S /N i = 3.5 threshold above which we consider that large TTVs do not prevent the detection of the planet (Leleu et al. 2023).The horizontal grey dashed line represent the radius cut-off above which we consider that the published mass-radius population is not affected by selection biases.
hydrogen-helium envelope of the smaller of the two planets.The impact energy of the impacting core is then added as a luminosity source over a smoothing timescale on the core-envelope boundary of the larger target, which can lead to radius inflation and gas loss, namely, impact stripping.For a full description of the technical implementation, we refer to Emsenhuber et al. (2021a).Here, we note that the perfect merging of volatiles, such as water, is a current model shortcoming and impact stripping of volatiles need to be included in future versions of the model.Nevertheless, the stripping of heavier water or carbon-bearing species is less efficient than that of the lighter hydrogen and helium (Burger et al. 2020).
The evolution of the 1000 NGPPS planetary systems, following 100 Myr of N-body integration with the aforementioned model, was re-calculated by Burn et al. (2024) considering an improved equation of state for water (Haldemann et al. 2020) for each simulated planet individually.If water or other ices were present from the formation modelling, they are mixed with any H/He left at this stage.For simplicity, all volatile species are modelled as water.The revised internal structure modelling results in significantly different, larger, radii of volatile-rich planet due to the lower density of supercritical water.Over typically 0.1 Gyr timescales, the loss of the mixture is calculated using a mass-weighted mass loss rate of the two constituents due to high-energy irradiation by the star (Kubyshkina et al. 2018;Kubyshkina & Fossati 2021;Johnstone 2020).The fraction of elements in the envelope different from hydrogen and helium, Z env , is kept constant motivated by numerical results in the regime of efficient mass loss (Johnstone 2020).Here, we show the planetary population at an evolutionary age of 5 Gyr after star formation; however, we caution that due to computational limitations, the dynamical state represents a 100 Myr-old system and, in some cases, slightly varying masses (prior to atmospheric mass loss).
Finally, we briefly discuss the origin of the trends with distance to mean-motion resonance.To make a comparison with the observed systems, we selected similar-sized, close-in planets, and display their properties in Fig. 5 as well as in Figs.C.1 and C.2. From the latter figures, we can see that more collisions occur after the end of the disk lifetime for the non-resonant population and that their volatile fraction is on average lower.In addition, they also have larger masses and (in order to meet the radius selection criterion) they have higher bulk densities.Article number, page 11 of 14

Fig. 1 :
Fig. 1: Full sample of 133 sub-Neptunes used in this study.The top-left panel shows the mass-radius relation.The (nearly)-resonant population is defined as ∆ MMR <0.05 (see Eq. 1), while the non-resonant population is defined as ∆ MMR >0.05.The black line is the sub-Neptune mass-radius relation from Parc et al. (2024).Other panels show the cumulative distributions for parameters of the planets or their host star.The p value given for each parameter is the probability that the distribution of that parameter is drawn form the same underlying distribution for the (nearly) resonant and non-resonant populations.The potential 2D correlation between parameters is explored in Fig. B.1.

Fig. 2 :
Fig. 2: Same as Fig. 1, for the controlled sample.The potential 2D correlation between parameters are explored in Fig. A.3.

Fig. 3 :
Fig. 3: Relative density of planets as a function of the distance to the closest MMR (left).Resonant chains and single planets are arbitrarily set at ∆ MMR = 0.001 and 100, respectively.The colors indicate the method used to obtain the mass.The dark grey area indicates the 1σ confidence interval for the local mean value of ρ rel , while the dotted lines show a local estimation of its scatter.Distribution of relative densities assuming gaussian distributions, binning by method and distance to the resonance (right): resonant, nearly resonant (0.001 < ∆ MMR < 0.05), non resonant (0.05 < ∆ MMR ), and single planets.

Fig. 4 :
Fig. 4: Minimal mutual inclination of transiting confirmed planets.The dark grey area indicates the 1σ confidence interval for the local mean value of ∆i min , while the lighter grey area shows a local estimation of its scatter.
Fig. A.1: Number of points taken by Mars 2022 by HARPS and ESPRESSO facilities for planets in the 2 to 4 R Earth range.Each ESPRESSO measurements are counted as five HARPS measurements to account for photon noise.The horizontal grey dashed line represents the radius cut-off above which we consider that the published mass-radius population is not affected by selection biases.
Fig. A.2: Radius of the Kepler TTV characterised planets the full sample that came either fromHadden & Lithwick (2017) orLeleu et al. (2023), as function of the S/N of individual transits (S /N i ).The vertical grey dashed line shows the S /N i = 3.5 threshold above which we consider that large TTVs do not prevent the detection of the planet(Leleu et al. 2023).The horizontal grey dashed line represent the radius cut-off above which we consider that the published mass-radius population is not affected by selection biases.

Fig
Fig. A.3: Controlled sample for planets in the 2.7 to 4 R Earth range.(Nearly-)resonant planets are shown by orange points, while non-resonant planets are shown by blue points.The p value shown on the diagonal are 1D Kolmogorov-Smirnov test, while the p value on the bottom-left triangle are 2D Kolmogorov-Smirnov tests(Peacock 1983).Median and uncertainties on the p value are estimated by drawing 1000 samples assuming a Gaussian distribution for the radius of each planet, then computing the 0.16, 0.5, and 0.84th quantiles of the resulting p value distribution.

Fig
Fig. B.1: Same as Fig A.3 but for the full sample shown in Fig. 1, for radius between 2 and 4 Earth radii.

Table A
Notes.RV planets with no marker have either been only observed by HARPS, or by simultaneously (i.e.overlapping points withing few tens of days) by HARPS and another telescope.K2-136 was first observed by HARPS-N then ESPRESSO, WASP-47 was first observed by Coralie, then by ESPRESSO, HD 73583 was first observed by Coralie, then by HARPS, and TOI-1231 was first observed by PFS.