© The Authors 2024

1 Introduction

Since the 1995 discovery of the first exoplanet orbiting a main sequence star (Queloz & Mayor 1995), there has been explosive growth in the number of detections of these objects. The results from the Kepler mission provided the first statistical picture of the properties of planetary systems (Borucki et al. 2010), and we can compare the Solar System with other star systems in our galaxy. The transit method used in this mission, and its extension, the NASA K2 mission (Howell et al. 2014), contributed to the determination of fundamental planet parameters such as radius and orbital period for more than 3500 planets. Follow-up radial velocity (RV) surveys and observed transit timing variations (TTVs) have helped to determine planet masses in an ever-growing number of systems.

In the early days of exoplanet exploration, the main challenge was to detect the gravitational acceleration a small planet induced on its host stars or other planets in the system with the available ground or orbital instruments. Therefore, the idea of the theoretical mass–radius relation (MRR) was appealing both for planet mass prediction and the development of the first theories of observed exoplanet interiors and atmospheres. The MRRs are typically based on exoplanet data collected in the main exoplanet catalogues, such as observed mass and radius distributions, and are designed to explain these observations. Recent studies (Chen & Kipping 2016; Weiss et al. 2013; Weiss & Marcy 2014; Wolfgang et al. 2016; Bashi et al. 2017; Otegi et al. 2020) used different approaches, such as power-law fitting to the whole known population, or suggesting transition regions in the populations.

Initial analyses of available data indicated a rich diversity in observed exoplanets (Santerne et al. 2019; Armstrong et al. 2020; Hoeijmakers et al. 2018; Winn et al. 2018; Espinoza et al. 2020; Demory et al. 2016; Evans et al. 2016). It was soon established that the MRR diagram presents a large scatter, making it difficult to disentangle subsets in the population. The power-law fitting and other simplified methods of producing MRRs began to face increasing difficulties in explaining the observed diversity. A different approach to suggesting a specific MRR is to take into account thermal equilibrium calculations and to predict the likely occurrence of particular species and mineralogy within the interior of a planet, along with their equations of state (EOS). This approach was implemented by Zeng et al. (2016, 2019), who calculated the MRR curves for prescribed planet bulk compositions, including compression. Both the Kepler and ongoing TESS (Ricker et al. 2016) mission discovered a large population of super-Earths and mini-Neptunes, sometimes at close-in orbital distances. The lack of a planet of analogous size (and possible composition) in the Solar System means that we cannot fully rely on knowledge of internal planetary structures in our neighbourhood to interpret the majority of known exoplanets.

While radius is the most straightforward to measure, the density and interior structures of exoplanets hold key information about their physical nature and origins. Carter et al. (2012), who present an example of a Kepler-36 system with the smallest fractional separation between orbits of two planets and the largest density contrasts between them, raised interesting questions. If the two planets formed in situ, the compositions and densities of the planets could have changed with time, for example because of the erosion of the atmosphere of the closest planet by instellation. On the other hand, formation at two widely separated locations in the protoplanetary disc, for example in the volatile-poor and volatile-rich regions, would suggest a ‘migration’ mechanism capable of drawing two planets together. Luque & Pallé (2022) proposed a density gap that separates rocky from water-rich exo-planets and their results are inconsistent with the previously known bimodal radius distribution arising from atmospheric loss of a hydrogen and helium envelope. These considerations call for a focus on planet density and interior studies.

To understand the interiors and constituents of exoplanets, it would be helpful not only to base our theories on mass and radius determinations of planets but also to gather insights from a planetary system as a whole. Exoplanet observations indicate a wide variety in the architecture of multi-planet systems (Lissauer et al. 2011; Fabrycky et al. 2014; Winn & Fabrycky 2015), but most known systems have planets at separations of less than the orbit of Mercury that are simultaneously much larger and heavier than our innermost planet. The Solar System was once again found to be quite different from the others. Furthermore, the observations hint at the existence of several intrasystem correlations (Millholland & Winn 2021; Otegi et al. 2022) also not observed in our system. Exoplanets within an exo-planetary system tend to be of similar size and mass and evenly spaced, and are often ordered in size and mass, as was first suggested by Weiss et al. (2018), Wang (2017), and Millholland et al. (2017). Similarity implies that P_i+1/P_i ≈ 1, where P is a chosen planet parameter, which can be mass, radius in pairs of adjacent planets, or period ratios for the second and first pair in triples. However, the authors of the several papers mentioned above, and of some mentioned below, formulated these trends differently using different restrictions posed on their samples, and sometimes looked only for similarity in one or two parameters. Weiss et al. (2018) found the correlation between the size of a planet compared to the next planet in the system in order of increasing orbital period for the sample, which they constructed solely from the Kepler photometry using both confirmed planets and planet candidates, but posing multiple restrictions on the accuracy of the determined planet radius in order to work with a high-purity planet sample. These latter authors also found similarity in the period ratio in the sample restricted to period ratios of smaller than 4. Wang (2017), working with a sample of 29 systems with at least three RV-detected planets, found 15 systems with planets with M_p<30 M_⊕ that show unusual similarity in mass. Most of their sample consists of non-transiting planets, which prevented similarity analyses in radius. Analysing a small sample of 89 planets with masses determined by TTVs only, Millholland et al. (2017) found that planets in a system have both masses and radii that are significantly more similar than in the mock systems assembled from randomly shuffling planets in their sample. Otegi et al. (2022) reported significant similarity in radius, mass, and density in planet pairs belonging to the same system, especially if the system hosts planets with M_p<100 M_⊕. These authors found that planetary systems tend to be more similar in radius than in mass, and that generally, the planets most similar in radius do not correspond to those similar in mass.

Together, all spotted correlations were collectively referred to as ‘peas in a pod’ trends in the architecture of planetary systems. These trends jointly suggest that planets similar in size can have similar bulk densities (and perhaps similar chemical compositions). The physical reality of these trends has been challenged and attributed to observational biases, for instance, arising from Kepler’s discovery efficiency at the detection threshold (Zhu 2020, but see Weiss & Petigura 2020). Murchikova & Tremaine (2020) showed that apparent patterns can arise in radius distribution models with simple assumptions, without requiring planet radii to be correlated between neighbours. On the other hand, some of these correlations have been reproduced in planet population simulations (He et al. 2019, 2020; Mishra et al. 2021). He et al. (2019) found that models that assume planets have independent periods and sizes are inconsistent with the observationally derived properties of multi-planet systems in the Kepler sample. These latter authors suggest a clustered model where planets in the system tend to be similar in radius within the cluster and periods are drawn conditionally on the period scale of the parent cluster. They restricted their analyses to modelling planets with R_p<10 R_⊕, and in order to adequately simulate mock systems, they drew planet masses using the MRR from the non-parametric model in Ning et al. (2018). He et al. (2019) did not analyse the masses for similarity. The same authors later showed that the observed patterns in the ordering of planet sizes and uniform spacings cannot be produced by the detection biases of the Kepler mission (He et al. 2020). Using a modern planet population-synthesis model and forward modelling of the detection biases and completeness of the Kepler survey, Mishra et al. (2021) revealed the peas in a pod trends present in their simulations. These authors produced samples of planets with masses of up to 10⁴ M_⊕ and showed that both their models – one that is unrestricted by observational biases and one that simulates the Kepler observational field – exhibit a high level of correlation in adjacent planet radii (R ≈ 0.64–0.75) and in mass (R ≈ 0.6− 0.76). The models also exhibit a weak (R = 0.25) to moderate (R = 0.55) correlation in period ratios for restricted and unrestricted simulated populations, respectively. They formulated these trends as such:

• Similarity in size: planets within a system tend to be either similar in size, in the sense described above, or ordered in size.

• Similarity in mass: planets within a system tend to be either similar or ordered in mass.

• Similarity in spacing: for systems with three or more planets, the spacing between a pair of adjacent planets is similar to the spacing between the next pair of adjacent planets.

The main challenge for evaluating the peas in a pod trend as an astrophysical feature could arise from the incompleteness of the observed exoplanet population. First of all, at present, and with sufficient precision, we can determine the radii of transiting planets only (Struve 1952). Further, the Kepler mission, which was for many years the most productive mission in providing fundamental exoplanet parameter data, had its observational limitations, especially for small planets, due to transit depth measurement constraints and difficulties in detecting exoplanets at orbital periods of greater than 200 days. As planet occurrence increases towards small size, and Earth-sized planets are very common (e.g. Petigura et al. 2013; Fulton et al. 2017; Burke et al. 2015; Mulders et al. 2018), it is likely that at least some of the known systems harbour additional undetected planets that could skew the observed peas in a pod trends.

The other important factor that could affect trends emerging from the known population is the accuracy of measurements of exoplanet parameters such as mass and radius. The determination of these parameters relies on measurements of the host star radius and mass and can be imprecise and inhomogeneous. Recently released results of the Gaia astrometric space mission (Gaia DR3), which provide the astrometry, photometry, and dispersed light for hundreds of millions of sources from the Gaia prism spectra and the spectrograph (Fouesneau et al. 2023), provide a homogeneous set of stellar radii. Initially, in this study, we worked with planet parameters retrieved from the exoplanet databases as described below in Sect. 2. We subsequently re-evaluated the radius trends using homogeneously derived stellar radii from the Gaia DR3.

Intrasystem correlations in planet radii and mass can reflect the degree of uniformity in the internal structures and compositions of planets within the same system. One might expect planets born in the same system to have a similar composition, especially if they are formed near one another. After all, the terrestrial planets in our Solar System and the TRAPPIST-1 system are remarkably uniform (Agol et al. 2021). Indeed, Otegi et al. (2022) found a strong correlation in bulk density between adjacent planets in multi-planet systems, especially around M dwarfs. However, bulk density is a simplistic proxy for interior composition and does not fully capture mass-dependent self-compression and pressure-temperature-related phase transitions and material changes. In a recent study, Rodríguez Martínez et al. (2023) investigated the core and water mass fractions for multi-planet systems orbiting M dwarfs and found no evidence of intrasystem uniformity (and found a hint of significant non-uniformity in some cases). Though the sample studied by these latter authors is small, their work indicates that at least some stars host rocky and volatile-rich planets side by side. It would be interesting to extend this analysis to a larger sample of carefully vetted multi-planet systems.

In the present paper, we revisit the peas in a pod trends, as formulated by Mishra et al. (2021), in the exoplanetary system architecture, using recent exoplanet and stellar data. We aim to assess whether or not the similarity of a pair of planets in one parameter (e.g. radius) can be a good predictor of their similarity in other parameters (e.g. mass or density). The paper is organised as follows: in Sect. 2, we explain the process of assembling this study sample and our methods of analysing the uniformity of planetary systems in planet parameters. In Sect. 3, we then report our findings in terms of the patterns we see in mass, radius, density, and period spacing. Our conclusions are discussed in Sect. 4 and are summarised in Sect. 5.

2 Methodology

2.1 The sample

The data for this study were initially queried from The NASA Exoplanet Archive¹ (Akeson et al. 2013), which was accessed on July 18, 2022. We cross-checked the retrieved parameter values with the Extrasolar Planets Encyclopaedia² (Schneider et al. 2011). Additionally, the methods of mass detection for planets in the sample were retrieved from papers that reported this parameter. The reported masses were cross-validated between the two databases in order to exclude parameters obtained from theoretical MRR. The majority of this study’s sample came from the Kepler and TESS missions. In order to explore the architecture of the known planetary systems, and to study them as a family of objects, it is critical to assemble a sample consisting of multi-planet systems, where the configuration of the system is at least partially established.

We assembled an initial sample of 528 exoplanets in 175 multi-planet systems by selecting exoplanets from the systems where at least two planets have measurements of mass and radius. Planets discovered by the direct imaging method, which typically probes young A-type stars and very large planets on wide orbits, often have conflicting reported parameters, such as mass and radius. Supermassive planets usually come from young unstable, often multi-star systems, and their parameters are observed with large uncertainties (Galicher et al. 2016). HR 8799, PDS 70, and TYC 8998-760-1 systems formally made it to the sample but were eliminated for these reasons. Therefore, the mass and density samples include only planets with masses measured by RV or TTVs, and the radius sample includes planets with radii measured from transits. In order to analyse planetary systems architecture, we include in our sample all known parameters of the host star and all planets. We then extracted separate samples for adjacent inner and outer planet parameters, if two consecutively located planets have them reported. For the analyses of features in the systems architecture such as mass and size distribution, and spacing, the subsamples were further vetted to include only the systems containing at least two adjacent planets with the known parameters of choice (P) and which have the sum of uncertainties σ_Psum ≤ 2P (i.e. σ_Psum is equivalent of 2σ, where σ is standard deviation). For period ratio analyses, systems from our main sample with at least three planets were included. For studying an intrasystem parameter’s similarity, we include all planet pairs or triples belonging to the same system.

We started with forming separate subsamples for masses (184 pairs), radii (269 pairs), calculated bulk densities (167 pairs) and period ratios (167 triples) of adjacent planets by filtering the data to consider only exoplanets satisfying our uncertainty criteria. Additionally, planets with uncertainty flags³ were excluded. In our subsamples, if the planetary system has more than 2 planets, all formed planet pairs are included in the correlation calculation, representing multiple data points. To further address using uncertainties in parameters in our sample, slightly larger than usual, we compared our methods with ones used by Otegi et al. (2022), and only a handful of planets in our sample are above the threshold suggested by them (see Fig. A.1, left). As we demonstrate below in Sect. 3, our outcomes were largely not affected even when the sample was further restricted for uncertainties even smaller than that.

It should be noted that as with any other, this study catalogue suffers from observational biases, and database errors and is essentially incomplete. In Sect. 3, we explore the effect of inhomogeneity in stellar parameter measurements on the results by constructing a new sample of planet radii. General Stellar Parametrizer from Photometry (GSP-Phot; Recio-Blanco et al. 2023), which is designed to produce a homogeneous catalogue of parameters for hundreds of millions of single non-variable stars based on their astrometry, photometry, and low-resolution BP/RP spectra, provides stellar radii (Andrae et al. 2023). Alternatively, Final Luminosity Age Mass Estimator (FLAME) (Creevey et al. 2023) estimate the stellar radii based on the stellar astrophysical parameters from either GSP-Phot or General Stellar Parametrizer from Spectroscopy (GSP-Spec) combined with the Gaia photometry and parallaxes. We retrieved the radii of the sample stars from FLAME and GSPPhot from Gaia DR3, which we retrieved from the Gaia Archive⁴, and calculated new radii using planet-to-star ratios that were queried from The Q1-Q17 Data Release 25 (DR 25) Kepler Objects-of-Interest (KOI) (Hoffman & Rowe 2017) if available, otherwise, the most recent values were queried from the NASA Exoplanet Archive. We used this new radii sample to check if the results of our analyses would differ from those performed on our initial sample. To address inhomogeneity in mass in our sample, we identify the method of mass determination for all planets in the sample with known mass, using both The NASA Exoplanet Archive and Extrasolar Planets Encyclopaedia. We checked that the distributions of planets with TTVs and RV-determined masses are similar in mass-radius space, see Fig. A.1, right. In this figure, all mass bins from smaller than 1 M_⊕ up to 10³ M_⊕ have planets whose masses were detected either with one or another method. There are no visible clusters in the distribution of planets’ mass determined by these two methods, and therefore, we assume any method-related mass bias to be small.

As a result, we reported our findings based on the best and most recent data to our knowledge. The data of the planets included in our main sample catalogue is accessible in ASCII format at the CDS. The code prepared here for assembling and analysing the sample can be found on GitHub⁵.

2.2 Statistical analyses

The Pearson R test measuring the strength of correlations between two variables (Pearson 1905) is a powerful method to investigate the correlation in two distributions. R values vary between −1 and +1 with zero implying no correlation and ±1 implying an exact linear relationship. The probability that an uncorrelated sample of the same size will show the same R is denoted by the p-value and indicates the significance of the results. The test assumes a normal distribution and a linear relationship between two variables, and therefore suffers a setback for small samples (N < 30; Jiang et al. 2020). Non-normality caused the correlation coefficient inflation, particularly in the case of extremely skewed distributions, but for large samples, this effect is less prominent (Bishara & Hittner 2015). The Pearson test can also fail if the dataset consists of data points gathered around the unity line at the low and high tail ends. If the pairs of very similar values belonging to the two populations, comparing in the test, occupy the high and low tails of distributions, it would lead to apparent correlation, and in the opposite, highly dissimilar values on the one tail of distribution would ‘smear’ out the possible correlation. Additionally, the Pearson correlation coefficient as a method has its limitations for searching for correlations in sufficiently large samples (typically N ≳ 500) or insufficiently small samples (N ≲ 30), such that the significance levels of Pearson’s R may be biased due to the sensitivity of p-value to N. Inherently, the p-value calculated along with Pearson R as a function solely depends on sample size and value of R, and can therefore be very small with large R and N. The confidence intervals, however, can be a reality check, showing the interval where the actual value of R might be at a certain degree of confidence. That is why we calculated them for R using Fisher’s z′ transformation (Fisher 1915), which is used due to the non-normality of the distributions.

Nevertheless, this method was suggested by multiple studies for finding a correlation in fundamental parameters of two adjacent planets within observed systems (Weiss et al. 2018; Jiang et al. 2020; Millholland & Winn 2021; Otegi et al. 2022).

Therefore, we also used the Pearson correlation test in this study and further discussed its caveats in Sect. 4.

First, we performed this test in the initial sample, comparing planet parameters such as mass, radius, density in planet pairs and period ratio in triples. In forming our samples, we chose to maximise the number of planets included, therefore imposing an uncertainty threshold that is larger than previous authors for similar studies (for example, Otegi et al. (2022) suggested as a cut-off σ_msum ≤ 1m and σ_reum ≤ 0.32r for mass and radius uncertainties, respectively). The number of planets in our sample with uncertainties larger than in previous studies is relatively small (see Fig. A.1). To account for the impact of large uncertainties, we implemented a routine that ‘weights’ the data points according to the measurement precision by performing multiple (10⁵) Monte Carlo experiments when the results are computed based on a repeated random uniform sampling of data points (x, y) from the observational error intervals (x ± dx, y ± dy). At every step during the simulations, we also performed the Pearson correlation test. Then we performed subsequent statistical analysis to determine the median values of the distribution. These median values can be considered the true values of the correlation coefficients. We note that the computation of the Pearson coefficients initially does not include the uncertainties, but our weighting or ‘error-accommodation’ technique, allows us to overcome this without complicated calculations.

To establish the null hypothesis as a comparison, we constructed test populations from the non-parametric model in Ning et al. (2018) which relates the planet’s radius to its mass. Using observed radii in our sample, we randomly draw test populations of mass and density in such a way, that any adjacent pair of planets do not violate stability in the system. For that, we assumed that mock planets occupy the same orbits as in the real systems, and calculated stability criteria as formulated in Deck et al. (2013): in the case of initially circular orbits two bodies in a planetary system should satisfy the following equation for semi-major axes: $$\frac{a_2-a_1}{a_1}\gtrsim 2.4{ϵ}_{\text{p}}^{1/3},$$(1)

where ϵ_p = (M_p(i)+M_p(i+1)))/M_★ is the total mass of the planets relative to the mass of the host star. This equation was used to determine if we produced any unstable planet pairs, and if so, we eliminated them from the calculation. This way of constructing the mock populations allowed us to preserve observed structure in the most precisely up-to-date determined fundamental parameters of exoplanet: in radius and period. After acquiring 10⁵ test populations, we performed statistical analyses and calculated mean or median values (in case of heavily skewed distributions). We used the test population to compare with the Pearson test results and consequently in further analyses of the sample.

Another way to address the intrasystem similarity as a trend is an evaluation of the mean normalised parameter dispersion (first formulated by Goldberg & Batygin (2022) as adopted a modified version of the mass uniformity metrics from Wang (2017), and subsequently used by Lammers et al. 2023): $$D=\frac{1}{N_{\text{sys}}}{\displaystyle \sum _{i=1}^{N_{\text{sys}}}\frac{{\sigma }_{p_i}}{\overline{p_i}}}=\frac{1}{N_{\text{sys}}}{\displaystyle \sum _{i=1}^{N_{\text{sys}}}\sqrt{\frac{{\displaystyle {\sum }_{j=1}^{N_{\text{p}{l}_i}}{\left(p_j-\overline{p_i}\right)}^2}}{{\overline{p_i}}^2\left(N_{\text{pl}{l}_i}-1\right)}}},$$(2)

where $N_{\text{p}{l}_i},\overline{p_i}$, and σ_pi are the number of planets, mean parameter value, and parameter’s standard deviation in the ith system, respectively. If the value of D is zero, the parameter’s dispersion is zero, and with rising values of D similarity in the system dissipates. Essentially, D represents the sum of standard deviations divided by means or coefficients of variation of parameters in all systems averaged by a number of systems. For example, a D = 0.5 signifies that the system-wide spread in the parameter of interest is on average half the mean value of that parameter in the given system. We calculate intrasystem parameter dispersion $\frac{{\sigma }_{p_i}}{\overline{p_i}}$ for masses, radii, densities and period ratios for each system in our sample, and D for these parameters for the entire sample and several subsamples.

The fundamental parameters of the protostar can lead to formation of planets with similar properties. For example, Millholland & Winn (2021) suggested that low-metallicity stars are harbouring planets more equal in size. We probed the hypothesis, that subsets of data points are correlated stronger than the entire sample, and we tested it in two ways. First, from the ratio of adjacent planet parameters (P_i+1/P_i) distributions we selected those planet pairs whose parameter ratio is closest to unity. By removing values farthest from unity, and performing the Pearson test at every step, we obtain the high-correlated (HC) subsam-ples of masses, radii, and densities in pairs and period ratios in triples with Pearson R ~ 0.95 (HCR, HCM, HCD, and HCP subsamples). We used these new subsamples for cross-checking analyses for mass and radius, and interior composition similarity for planets belonging to the same subsample.

Second, we performed a ‘moving window’ test, similar to Otegi et al. (2022): (i) we worked with separated subsamples for masses, radii, densities and period ratios of adjacent planets; (ii) we sorted all planet pairs in each subsample ascending for a given stellar parameter, then selected the first 40 data points and perform the Pearson test; (iii) we repeated the same procedure with a continuously moving 40-pair subsample window until the entire subsample is covered. This technique allows us to identify the ranges in stellar parameters, where subsamples show an unusual correlation. The significance of this tendency was assessed by a bootstrap test: the correlation coefficients of the partial sample were compared with a median value of those of distributions from 10⁵ simulations, in which we calculated the Pearson coefficients for a random subset of the same size drawn from the entire sample. The bootstrap technique is typically used to assess non-normal data distributions, as in this study, and we used this bootstrap experiment as an additional test for the previously obtained Pearson correlation coefficients for mass, radius, density and period ratio subsamples.

3 Results

3.1 Mass, radius, density, and period ratio analyses

We analysed the peas in a pod trend in our sample. First, we tested whether the pairs of adjacent planets in our sample have similar mass, radius and density, or whether the planet triples show similarity in period ratios, by using the Pearson correlation test. Second, we assessed intrasystem similarity, calculating the parameter dispersion for each individual system and the mean value for the entire sample. As we aimed to assess similarity patterns in the architecture of planetary systems in the general interpretation formulated by Mishra et al. (2021), we needed to first evaluate the sample of all multi-planet systems with known fundamental parameters for at least two planets, and then test our findings by posing multiple restrictions on the sample, including the error truncation. Although it seems more practical only to study planets with well-determined (>3σ) masses, imposing such a significance criterion biases the sample toward more massive planets at a given radius. This is why we suggested the error-accommodation routine that “weights” the data points according to the measurement precision to fully accommodate the error of the different parameters. The routine allowed us to consolidate a sample large enough to represent multi-planet systems statistics. Our sample consists of data for 306 planets in 122 systems with determined masses, 414 planets in 145 systems with determined radii, and 353 planets in 93 systems with determined periods.

3.1.1 Correlation in the parameters of adjacent planets

In Fig. 1, the scatter graphs represent the inner planet parameter plotted against the outer planet parameter: mass, radius and density, as well as period ratio between the middle and inner planet plotted against period ratio between the outer and middle planet in triples. The outer planets with larger radius comprise 62% of the sample, and this “ordered in size” trend was previously suggested by Weiss et al. (2018). 52.7% of the mass pairs are above the unity line, indicating that there is no obvious trend. Only 35.3% of the density pairs are above the 1:1 line, which means that more distant planets tend to be less dense. For period ratios, spacing in triples lies almost equally above (50.3%) and below (49.7%) the unity line.

We find that there is a moderate correlation for radii of adjacent planets with the Pearson R = 0.516 and p-value ~10⁻¹⁹, trumpeting R as significant, and a weak correlation for densities R = 0.265 with p = 0.0005, which is less than adopted α = 0.05. There is no correlation for masses R = −0.007 or period ratios R = −0.017, and the 95% confidence intervals for these parameters do not give grounds to assume even a weak correlation for the planets in the main sample. The results for the sample and after running 10⁵ simulations using the error-accommodating technique (see Sect. 2) do not differ much and are represented in Table 1 and in Fig. A.2. We repeated the analysis with a restricted subsample with higher parameter precision of σ_Psum ≤ 0.5P, and it showed similar results, except for the density subsample, which revealed a moderate correlation with R = 0.654 in the small subsample (n = 44). We created 10⁵ bootstrapped samples with 40 random pairs from each parameter subsample and subsequently performed correlation tests. The bootstrap test results confirmed our initial outcome, and along with the error precision test results are also represented in Table 1. and additionally in Fig. A.6. As the last test, we constructed 10⁵ mass and density test populations using the non-parametric model (Ning et al. 2018) for the observed radii in our sample and repeated the calculation of the Pearson coefficients. For mass and density, the expected values in the main sample are within 1σ, and for the error precision of σ_Psum ≤ 0.5P we found that the expected value for density correlation is larger than for the main sample (see Table 1).

To examine the robustness of our results concerning the similarity of adjacent planets in the sample, we conducted several tests along with those already mentioned. Our original sample consists of data for 246 planet pairs from the Kepler survey and 47 pairs of planets discovered by TESS and the rest discovered by others, and therefore the sample could suffer from instrument-related biases. To eliminate this effect, we analysed the data from Kepler separately and listed results in Table 1). The Kepler sub-sample shows similar results as the main sample in the Pearson test. Then, we tested the sensitivity of the results on the sample size and parameter error precision by restricting the entire sample to σ_Psum/P ≤ 1/2. This test clarifies that the results do not deviate much with increasing data accuracy, except for densities, which show a moderate correlation in the restricted sample. Moreover, the test population produced from the non-parametric MRR shows analogous results in mass and density for the Kepler sample.

Finally, we tested the sample restrictions often found in the literature: M_p < 100 M_⊕ and R_p < 10 R_⊕ (see e.g. Wang 2017; Otegi et al. 2022). This restricted subsample consists of 73 pairs of planets with error precision adopted in this study and the Pearson coefficients show a medium to strong correlation in radius (R = 0.69), a weak to medium correlation in mass (R = 0.35) and very weak correlation in density (R = 0.16).

Fig. 1 Masses, radii, densities, and period ratios of adjacent planets. The scatter graph for mass (top left), radius (top right), and density (bottom left) of a planet against the same parameter of the next planet farther from the star is plotted with observational errors indicated by the grey lines. The bottom right panel shows the period ratio between the middle and inner planet against that between the outer and middle planet in triples. The blue and red dot-shaped markers correspond to the sample and Solar System planet pairs, respectively. The black dashed line is the 1:1 line. The observational errors could be smaller than the marker size. In the bottom right corner of each plot, the cumulative distribution function is plotted for the decimal logarithm of each planet parameter ratio (x-axis) versus counts (y-axis). In these subplots, the black vertical dashed line is unity and the red horizontal dashed line is where the cumulative distribution function reaches the 50th percentile.

3.1.2 Intrasystem dispersion

For the entire sample of systems with at least two planets, we found that the average intrasystem dispersion D in radius is 0.317, and it is larger for mass (D = 0.54), for density (D = 0.584), which indicates that radii of planets belong to one system are more uniform than masses and densities. We run 10⁵ simulations using the error-accommodating test for these parameters, and the results show a slightly larger dissimilarity in mass and density. Using the MRR from the non-parametric model (Ning et al. 2018) we simulated 10⁵ mass and density test populations for the observed radii in our sample and repeated the calculation of D. These results show that the expected values of D point towards an even larger dispersion than the error accommodation test, but it is still less than the simulated mean for mass and density (see Fig A.4, two left plots. The repeated analysis with a restricted subsample with the parameter precision of σ_Psum ≤ 0.5P showed similar results.

From the main sample, we selected systems with at least 3 planets and repeated the routine, adding calculations for the period ratio between two consecutive pairs of planets in the system. Again, the radius is the only parameter that is consistently more uniform than the other parameters we investigated. Remarkably, the systems with at least three planets show less uniformity in mass and density, than the entire sample (D≈0.54–0.58 for N_pps>1 and D≈0.71–0.74 for N_pps>2, where N_pps is the number of planets per system). It could be the indication, that 2-planet systems, which is a majority in our sample, are more uniform than systems with higher multiplicity, but we noted that such systems can be incomplete, and we discuss the completeness of systems in our sample in Sect. 4. It also confirms our conclusion based on the Pearson correlation test and shows that the radii tend to be similar both in pairs of adjacent planets and within the system they belong to. However other parameters show more dissimilarity in both cases. The calculation of D for the main sample and subsamples is represented in Table 2.

In Fig. 2, the scatter graphs represent all systems in our sample, for which intrasystem dispersion in mass, radius and density was calculated. The number of systems which have low dispersion in radius is larger than that in mass. Low dispersion in density appeared in systems with large dispersion in mass. Among systems with both small dispersion in mass and radius, there are systems with large dispersion in density, and in systems with dissimilar mass or radius, density can be similar. More systems are below the unity line, which not only indicates that more systems are similar in radius than in mass, but the systems with similar planet masses differ from those hosting planets with similar radii.

Usually, the significance of the intrasystem uniformity D is assessed by comparing its value to a Monte Carlo distribution of values generated by randomly shuffling planets among the set of systems or drawing one planet parameter randomly from the observed population (e.g., Millholland et al. 2017; Wang 2017; Otegi et al. 2022; Lammers et al. 2023). The former approach can lead to instability in the modelled systems due to the adjacent pair of mock planets ending up too close to each other. The latter approach could produce non-physical systems of planets with an unnatural set of fundamental parameters unlike any real system, and the drawback of this was described by Millholland & Winn (2021) as such: the erasing of the structure in period and radius distributions would lead to an overestimate of the strength of intrasystem uniformity. We therefore chose another method of producing test populations which is necessary for the null hypothesis testing. Instead, for mass and density, we produced multiple (10⁵) test populations using the non-parametric MRR (Ning et al. 2018) and observed radii in our sample and subsamples, retaining only planets that satisfied dynamical stability requirements (see Sect. 2). We tested intrasystem similarity by measuring the dispersion in fundamental parameters in our sample, including only the systems with planets M_pl < 100 M_⊕. In addition to the main and mass-restricted sample, we also computed intrasystem dispersions in the TTV and RV subsamples separately. For masses, densities and period ratios the values of D in all cases were consistently larger than 0.5 with the exception of D in mass in restricted M_pl < 100 M_⊕ subsample with D = 0.476±0.013. In this study, for most subsamples and the entire sample, D in mass is typically larger than 0.5, whereas our calculations show that dispersion in density is typically larger than 0.6 and exceeds the D value for mass and radius.

Fig. 2 Intrasystem dispersion in mass, radius, and density. The x-axis shows dispersion in mass for systems in the main sample, the y-axis represents dispersion in radii, and the colour of the dots corresponds to dispersion in density. The red colours indicate similarity in density, and the blue colours indicate dissimilarity. The black dashed line is 1:1.

Fig. 3 Sample on the mass–density diagram. The black spheres represent the sample’s planets and the red spheres are the Solar System planets, plotted for comparison. Dashed black lines connect a pair of adjacent planets in the “high correlated in radius” subsample. The coloured interior curves (from red to light blue) are taken from Zeng et al. (2016) and correspond to different calculated MRR for several potential 2-layer structures of exoplanets, representing different combinations of constituents of the interior: Fe, MgSiO3, and H2O. The error bars are plotted for planets comprising the “high correlated in radius” subsample and omitted for others for the sake of clarity.

3.2 Planets with similar sizes can have dissimilar masses, densities, and interior structure

From the distributions of parameter ratios (P_i+1/P_i) of adjacent pairs, by sequentially removing values farthest from unity and calculating the Pearson coefficients at every step, we obtained the highly correlated subsamples of masses (HCM), radii (HCR), and densities (HCD) in pairs and period ratios (HCP) in triples with Pearson R ~ 0.95 (see Table 1 under “High correlated subsamples” section). The HCM subsample consists of only 67 pairs out of 184, which indicates that a small number of planets are similar in mass with their subsequently located neighbour. The HCR subsample, showing a high correlation in radii, consists of 149 pairs out of 269, so more than half of the sample has a small dispersion around the unity line for radius ratio. The HCD subsample consists of 69 pairs out of 167, and the HCP subsample is even smaller (43 pairs out of 167).

We proceed to study the HCR subsample further. The HCR subsample correlation in radii of adjacent planets is represented in Fig. A.7, top left, along with one random realisation of bootstrapping pairs of planets from the main sample. We further investigated intrasystem dispersion in mass and density in the HCR subsample by performing the error-accommodation routine and comparing obtained D_err with the distribution of D_sim from 10⁵ random simulated mass and density from non-parametric MRR, and plotted the results in Fig. A.4. The observed dispersions in density and mass (D_err=0.437 and D_err=0.388, respectively) follow the same trend as the entire sample showing smaller values than the expected values D_sim, however, at the same time being larger than dispersion in radius in this subsample, which is the lowest among all samples (D_err=0.116).

In Fig. 3, we plotted the entire sample in the mass–density diagram and connected a pair of planets from the HCR subsample with dashed lines. To guide the eye, the coloured curves in this figure correspond to several representative theoretical MRRs for 2-layer exoplanet structural models including rocky and volatile components (Zeng et al. 2016). Even when considering the radius and mass measurement errors, adjacent planets connected by dotted lines from the HCR subsample appear to be to a large extent well separated in their mass and bulk densities, and often straddle across distinct bulk composition curves.

To assess how strongly each parameter differs among a pair of planets, we introduce a parameter gap, ɡ, to characterise the relative difference between a pair of values of a certain parameter. The term relative difference is used in quantitative science to compare two quantities while taking into account the parameter “sizes” being compared, that is dividing by a standard, reference or starting value, or as formulated in Törnqvist et al. (1985), an indicator of relative difference, that is a function dependent only on the ratio of parameter values, usually the arithmetic, geometric, harmonic mean, etc. The result is expressed as a ratio and is unit-less and independent of values in the sample. Therefore, the parameter gap ɡ is an absolute value of the difference between parameter values divided by the indicator of relative difference. In this case, we choose a moment mean change of second order f(x, y) = (1/2(x² + (y²))^1/2 as an indicator of change: $$g=\frac{|x-y|}{{\left(x^2+y^2\right)}^{1/2}},$$(3)

where x and y denote the values being compared. The metric ɡ varies between 0 and 1 and represents a scale-independent, normalised distance in the parameters. At zero, the parameter gap indicates maximal similarity, and the closer this parameter is to unity, the higher is dissimilarity in the planets’ pair.

To quantify the degree of mass and density correlation in the HCR subsample we first calculated the parameter gap ɡ using Eq. (3), where x and y are parameters of the first and the second planet, respectively.

In Fig. 4, we plotted empirical cumulative distribution functions (ECDF) of g in radius, mass and density of adjacent planets in the HCR subsample. ECDF represents the proportion of observations falling below each unique value in a dataset and has the advantage that each parameter ɡ is visualised directly and aids in direct comparisons between these distributions. From this figure, it is apparent that the parameter gap in radius is small for the majority of the sample. However, mass and density gaps vary in the subsample significantly, and only a third of the subsample pairs show similar gaps in all three parameters. We further probed the mass-density relationship by performing the Pearson correlation test for the HCR subsample for masses and densities of planet pairs (91 and 84, respectively) and found a weak-to-moderate correlation for both parameters instead of an expected high correlation, as for radii (see Table 1 under “HCR subsample correlation” section). That is, in most cases, in a pair of planets from these samples, masses and bulk densities appear different, and only a very small population of adjacent planets, relative to the whole sample, will show uniformity in all three parameters.

If we place the planet pairs from the HCR sample on a fine mesh of interior structure curves in mass-density space, as shown in Fig. 3, it becomes visually apparent that most pairs of planets with very similar radii straddle widely separated compositional curves, implying distinct interior structures and bulk compositions. Cases like Earth-Venus and the TRAPPIST planets, which are strung along roughly the same curves, are the exceptions.

To quantify the separation between planets in pairs in the MRR grid, we calculate the interior distance between the MRR curves that best describe each planet. For planets in our sample, that belong to the area in the mass–density diagram described by the aforementioned MRR curves, we first calculated the closest curve to each planet. To illustrate the concept, in Fig. 5 we plotted a pair of hypothetical planets with similar sizes⁶. We locate the curve nearest each planet, then designate the closer of the two as the reference curve for the pair. The density values corresponding to the planet’s masses along the reference MRR curve are taken as reference densities, from which we calculate the signed difference to the actual planet densities using the gap parameter as defined in Eq. (3) (without the absolute value sign in the numerator). The mean of the two gap distances is the interior distance (I) for the planet pair: $$I=\frac{\left|I_{\text{p}1}+I_{\text{p}2}\right|}{2},$$(4)

where I_pi is the density gap between the planet and the reference MRR curve.

We compute I for all planet pairs with at least one planet within the area covered by the grid of two-layer interiors models from Zeng et al. (2016), which are valid for planets with M_pl ≤ 32 M_⊕. As a result, we excluded giant planets, as well as pairs where both planets have I_pi ≥ 0.1. This value was chosen because if both planets are farther than 0.1 from any MRR curve in our set, there is no particular curve which can represent the interior composition of any planet in the pair. In this case, the resulting average I would be large and would not indicate a difference between planets, but rather their mutual separation from the nearest interior curve.

The results of the I calculation are presented in Fig. 6. We noted that only TRAPPIST-1 shows intrasystem similarity in the interiors of four planets (with I ~ 0.01–0.017), but several pairs of more massive planets also show similarity. We used the mesh of MRR curves that has a 10% increase in each given chemical compound content from 100% MgSiO₃ to 100% Fe in the direction of increasing density and to 100%H₂O in the opposite direction. In this case, the values of I ≥0.1 effectively put a pair of planets on curves that will be at least three curves apart and therefore differ in composition by ≥30%. There are 15 pairs of planets with I < 0.1, indicating they could have relatively similar interiors. There are 26 pairs that have I ≥ 0.1, and therefore probably differ substantially in the interior structure. This suggests that neighbouring planets with similar radii often have vastly different compositions, interior structures, and volatile contents.

We proceed with the error-accommodation routine, and this test showed that the mean value I_err and the corresponding standard deviation for some of the pairs similar in the interior could have I>0.1, but for pairs with dissimilar interior I_err can be even larger than the observed I. A couple of planet pairs exhibiting the greatest compositional similarity, Kepler-3 b and c and Kepler-176 c and d, have their interior distance I_err inflated to >0.1 due to the large observational errors. Otherwise, the error-accommodation routine shows similar results. We tested the significance of our results by comparing them to the 10⁵ test populations in the same way as previously and obtained I_sim. In Fig. 7, we plotted observed I against expected I_sim. We examined trends with host star T_eff and noted that cold stars show more or equal similarity in their planets’ interior than expected, and hotter stars hosted planets with less similarity in the interior than expected. All results of the I calculation are presented in Table 3, where along with the mean value of I_sim, we listed also the median values due to most distributions in the simulations being non-symmetric.

Fig. 4 ECDF of the parameter gap ɡ in mass, radius, and density in the HCR subsample. The x-axis shows ɡ in a certain parameter, the y-axis represents the proportion of the gaps in pairs that reached a certain value.

Fig. 5 Interior difference between a pair of adjacent planets. In this illustrative example, planet P1 lies closest to the blue curve (distance L0) whereas planet P2 is closer to the red curve (distance L1). The blue curve, corresponding to 30% Fe and 70% MgSiO3, is chosen as the reference curve for the pair. This means if P2 had the same bulk composition as P1, P2’s bulk density should be given by the blue curve corresponding to its mass (red circle). Ip1 and Ip2 (shown in the thin red line) are measured as the vertical gap distance between planet P1 and planet P2 to the reference curve at their respective masses. This figure shows LTT 1445A b and c planets, and they are separated by an interior distance I = 0.053.

Fig. 6 HCR sample on the mass-density diagram. The dots, connected by black dashed lines, are planet pairs that belong to this sample. The coloured dots represent those planets, for which we calculated interior distance I. The red colours indicate a pair of planets at a short interior distance, which means they are probably similar in the interior structure, and blue colours indicate planets that show larger separation in the interior from their neighbour. Light-grey curves represent the MRR interior curves calculated by Zeng et al. (2016). From the baseline of a silicate planet with 100% MgSiO3, the red points in the direction of increasing Fe content, and each curve represents a 10% increase in Fe content up to 100% Fe. Similarly, the blue arrow points in the direction of increasing water content at 10% increments per curve. Solar System planets are added for comparison and plotted as purple stars.

Fig. 7 Interior distance in HCR sample systems. Observed I in the x-axis is plotted against expected I from non-parametric MRR (Ning et al. 2018) at the y-axis. The colours of the dots indicate stellar Teff of the system’s host star. The black dashed line is unity.

3.3 Similarity in the systems with similar stellar properties

We explored the possibility that the strength of the peas in a pod trend for adjacent planet pairs masses, radii, bulk densities and period ratios in triples depends on the properties of the host star. Such a correlation could arise from detection biases or from astrophysics. For example, Zhu (2020) mentioned that smaller planets are more easily detected around more quiet stars. On the other hand, the host star properties, such as stellar mass or metal-licity should, in principle, affect planet formation and evolution leading to certain planet interior structures as the end product. We used the moving window test and found that certain parts of the sample show signs of possible enhanced underlying correlation for masses and densities of adjacent planets. For planetary radii, the R and p-value somewhat fluctuate around the previously obtained values for the main radius subsample, giving R values close to or inside of 95% CI, and therefore the ‘40 radius pairs’ subsamples do not correlate any stronger than the entire dataset (see Fig. A.5). The same test revealed no dependency on stellar parameters for the period ratios in triples, with the maximum value of R < 0.5.

Masses of planet pairs around cool host stars (T_eff ∊ [2566; 4398] K) are correlated especially well (see Table 1), otherwise, the high R values also are typical for planet pairs around small, low mass, slightly metal-poor stars. Densities of planet pairs orbiting host stars with metallicity lower than solar (MT ∊ [−0.210; −0.010] [Fe/H]) and of old age (A ∊ [6.03; 12.3] Gyr) show a high correlation. Figure 8 shows p-value dynamics for these subsamples in coloured lines for density depending on host star age or metallicity and for masses depending on host star T_eff. Figure A.7 demonstrates these subsamples along with one random run of the bootstrap from the whole sub-sample. Additionally, we compared the resulting R values with the bootstrap test distributions in Fig. A.6. These R are far removed from the bootstrapped distributions and are therefore highly significant.

Fig. 8 p-value dynamics for subsets of data points. Coloured lines indicate p-value changes during the moving window test for density depending on host star age (green) or metallicity (grey) and for masses depending on host star Teff (blue). The light-coloured shades under the curves roughly indicate these subsamples in corresponding colours.

3.4 Assessment of data homogeneousness

Here we aimed at using a large and diverse sample, which is why we used the values collected by the databases from diverse sources in the literature. While this maximises the size of our sample, it also means the measurements are inhomogeneous. Homogenous stellar data could greatly help to minimise selection biases from the observations. The planetary radii uncertainties for transiting planets are dominated by the uncertainties in the stellar radii due to the transit depth δ = (R_p/R_✶)², where R_p and R_✶ are radii of planet and star respectively, and stellar radius determination can greatly affect the known value of the individual planet radius (Muirhead et al. 2012). To investigate the impact of using a set of homogeneously determined system parameters on our statistical results, we calculated new planetary radii from Gaia DR3 stellar data. Namely, the radii of the sample stars from FLAME and GSPPhot and associated uncertainties were retrieved. A planet-to-star ratio for every data point of the constructed radius subsample was queried from the NASA Exoplanet Archive, and the most recently reported value was used. The majority of the Kepler planet planet-to-star radius ratios came from the DR 25 KOI catalogue (Hoffman & Rowe 2017).

From 256 initial data points for the radii correlation analysis, we recalculated radii for 227 pairs for all systems that have radius data of the star from FLAME and GSPPhot. The results for this 227 pairs sample are represented in Table 1 under the section “Radii from Gaia DR3 correlation”. 62% of the radius pairs are above the 1:1 line, as for previous radii values, the retrieved Pearson R is lower than our previous result for the initial sample, and the p-values indicate slightly lower confidence in the results. The newly obtained radii of adjacent planets still appear to be moderately correlated, but with R ≤ 0.454 which is less than the previously obtained value. We proceed to calculate dispersion in radius D, using these newly obtained radii from Gaia DR3. We found that D shows similarly low values for both FLAME and GSPPhot inferred radii, and the error-accommodation test shows the results very close to the obtained values. The results can be found in Table 2 under the section “D for radius from Gaia DR3”. Therefore, based on these results from the Gaia homogenous sample, we consider our initial sample to be relatively homogenous in radius as well. Moreover, as we found no evidence that masses of planets in our sample determined by two different methods: RV and TTVs, are distributed differently (see Fig. A.1) we consider, albeit with less confidence, our sample relatively homogenous in mass.

4 Discussion

4.1 Radii are more uniform in pairs and systemwide than any other parameter

We compared key properties of adjacent planets, as well as their intrasystem dispersion, in a large sample of multi-planet systems. We found that mass, density, and period ratios are generally less uniform than radius. Using the Pearson correlation test, we found moderate correlation for radii, weak correlation for densities, and no correlation for masses of adjacent planets and period ratios in triples in the entire sample. On the other hand, the dispersion metric indicates increased system-wide uniformity in mass and density compared to the null hypothesis (significant to ~3σ). Moreover, among systems with high levels of radius correlation, masses and densities are not similar to nearly the same degree. The results of this study could further hint at why the MRR diagram presents a large scatter. Below we compare our results with existing observational and theoretical studies on peas in a pod patterns from the literature.

4.1.1 Radius trend

The radii of planet pairs in our main sample as measured by the Pearson coefficient have consistently shown medium to strong correlation (R = 0.516). They are fully consistent with those of Weiss et al. (2018, R = 0.62 for 355 multi-planet systems) and Millholland & Winn (2021, R = 0.444 for 520 systems). Furthermore, when using planetary radii based on homogeneously derived stellar radii from Gαiα DR3, we show that the radii correlation is slightly diminished (R ≥ 0.45), but still highly significant.

4.1.2 Period ratio trend

Weiss et al. (2018) found consecutive planets tend to have regular period spacings, as indicated by the Pearson coefficient (R = 0.46 in adjacent period ratios for a sample of 104 systems with more than 3 planets). This is apparently in tension with the lack of correlation in period ratios in our main sample. However, it is worth noting that the systems studied by Weiss et al. (2018) excluded planet pairs with period ratios >4. As demonstrated in Jiang et al. (2020) and Otegi et al. (2022), imposing such a restriction tends to inflate the significance of the correlations.

When focusing on compact systems with small planets where all period ratios pr ≤ 4 and radii R_p ≤ 6 R_⊕, Jiang et al. (2020) reported significant spacing correlation. We repeat our analysis on subsamples formed by enforcing the same upper bounds in period ratio and radii. Table 4 presents the corresponding R and p-values, which show no evidence of significant correlation. Such discrepancy could be in part due to differences in the underlying samples: our sample is constructed from systems containing planets where at least two of them have known radii and masses, whereas other authors have generally not demanded that masses are known. On the other hand, we do observe a reduction in the system-wide dispersions for the subsamples with pr ≤ 4. The TTV subsample also presents lower D in pr, which could be explained by the systems being preferentially near mean motion resonances (Jiang et al. 2020).

In planet synthesis models of Mishra et al. (2021), the underlying population of simulated planets with pr <4 exhibits a high degree of spacing correlation (R = 0.55). However, Keplerlike transit detection biases substantially reduce the apparent correlation to R = 0.25. Therefore, it is possible that the true spacing regularities are greater than they seem in the observed sample.

4.1.3 Mass trend

Since precise masses are difficult to measure, the peas in a pod trend in mass has been more challenging to study. The absence of effective ways to deal with measurement errors in most similarity metrics means samples are often formed with strong precision constraints. This leads to limited sample sizes, and potentially biases samples towards more massive planets at a given radius. Furthermore, the use of TTV versus RV masses is a topic of ongoing debate. The RV signal induced by small planets is often less than 1 m s⁻¹, which is below the detection threshold of current state-of-the-art instruments and below typical noise floors from stellar activity (~2 ms⁻¹) and photon noise (~2 m s⁻¹, Weiss & Marcy 2014). The TTV method is also problematic in that it is only useful for planet pairs in near-resonance, and can overestimate masses for planets with high orbital eccentricities (Hadden & Lithwick 2014; Yoffe et al. 2021).

Many authors have explored intrasystem trends in planetary mass and found evidence of relatively significant similarity (from ~2.5 to 8σ). Millholland et al. (2017) analysed 37 systems with planet masses derived from TTV (all below 50 M_⊕) and found a smaller dispersion within systems relative to random shuffling to high significance (~8σ). Citing concerns about selection biases towards resonant orbital configurations in TTV systems, Wang (2017) used 29 systems with RV planets and showed that a mass threshold (at 30 M_⊕) is required to reproduce the mass uniformity to 3σ. When systems with gas giants ≳100 M_≳ are included (representing half of the sample), the system-wide mass dispersion is indistinguishable from random. Goyal & Wang (2022) studied 17 systems with masses derived exclusively from RV measurements without imposing mass limits. They used the Gini index⁷ as their similarity metric and found a hint of intrasystem mass uniformity, albeit only to modest significance (~2.5σ). Recently, Lammers et al. (2023) computed mass dispersion for 100 systems harbouring only planets with M_p < 30 M_⊕ and found uniformity significant to ~6σ.

For this study, we assembled a relatively large mass sample using a combination of TTV and RV measurements and lenient error precision criteria. In computing the similarity metrics, the errors are accounted for using the error-accommodation routine. Our main sample does not appear to present mass uniformity to the same extent as previous works. However, removing the systems harbouring massive planets (> 100 M_⊕) from our analysis does strengthen correlations in mass. We tested this by selecting 73 pairs of planets with M_p < 100 M_⊕ and R_p < 10 R_⊕ and error precision thresholds as used for the main sample. In this subsample, the Pearson coefficients show a medium to strong correlation in radius (R = 0.69) and a weak to medium correlation in mass (R = 0.35). Correspondingly, the system-wide dispersion in mass also contracted compared to the main sample (see Table 2). Nevertheless, they remain far below the degree of similarity observed in radius. Intriguingly, the density correlation is weakened (R = 0.16) in the mass-restricted sample compared to the main sample, and the dispersion is not notably changed.

On the other hand, the physical motivation for excluding a certain type of planet from characterising general trends in planetary system architecture is less clear. Gas giants may play a fundamental role in sculpting planetary systems (e.g., Bitsch & Izidoro 2023), and large planets are commonplace, especially around FGK stars (e.g., Mayor et al. 2011; Hsu et al. 2019; Fulton et al. 2021) and including in our own system. In Wang (2017)’s entire sample, half the systems contained planets ≳ 100 M_⊕. Therefore, the results from our main sample are also instructive. Moreover, the population synthesis simulations of Mishra et al. (2021) did not require a mass cut-off to produce a mass correlation in adjacent planet pairs comparable to that in radius (both with R ~ 0.6). By comparing this prediction with the observed results from our main sample, we can test the planet formation models. The fact that they contradict our observed results seems to undermine Mishra et al. (2021)’s conclusion that radii similarity is rooted in mass similarity. Nevertheless, it is also possible that the observed and theoretical samples are not directly comparable because detection biases associated with planetary mass measurements were not modelled in Mishra et al. (2021).

4.1.4 Taken together: Comparison with Otegi et al. (2022)

The most comparable study to ours is Otegi et al. (2022). In a sample of 144 planet pairs from 48 multi-planet systems with well-determined masses and radii, Otegi et al. (2022) investigated similarities in radii, masses, as well as bulk densities, finding significant pair-wise correlations in all three parameters. However, the correlation is much weaker in mass than it is in radii and density, which is in qualitative agreement with our conclusions.

A key difference in our samples could contribute to the discrepancy in the mass correlations. Our sample is constructed using newer data and a more relaxed selection criteria in parameter errors (which we marginalise over in our analysis) and is therefore also larger (>180 systems). To test the effect of using data with the same precision as Otegi et al. (2022), we repeated our analysis on a subsample selected with Otegi et al. (2022)’s error criteria. This resulted in a greatly strengthened density correlation, but no meaningful change to that in mass (see Table 1 under “Error precision tests”)⁸. We suggest that recent updates to the exoplanet inventory or differences in the definitions of the similarity metrics could be responsible for the remaining discrepancy.

Otegi et al. (2022) found that most systems display higher dispersion in mass compared with radius, which they posited to be a consequence of densities tending to be similar in a given planetary system. We directly investigated density dispersions in all the systems (see Fig. 2). We found that systems with low dispersions in mass and radius tend to have low dispersions in density, but systems with higher dispersions in mass and/or radius show a variety of density dispersions. Moreover, our findings are compatible with Otegi et al. (2022)’s realisation that systems with the greatest intrasystem uniformity in radius are often not the most uniform in mass, and vice versa.

Similar to Wang (2017), Otegi et al. (2022) found a transition in the peas in a pod behaviour at around Saturn-mass (~ 100 M_⊕, ~10 R_⊕), such that excluding the giant planets improved uniformity in radius and mass. This is also observed in our sample, as discussed in Sect. 4.1.3. However, we find no evidence for an analogous behaviour in density, and it is not mentioned by Otegi et al. (2022).

4.1.5 Dependence on stellar parameters

Millholland et al. (2017) stated that previous works have shown that intrasystem radius uniformity depends only weakly on stellar properties. On the contrary, Otegi et al. (2022) found that planetary systems around less massive stars tend to be more uniform in mass, radius, and density, since the planets tend to be less massive. In our study, we found that the correlation in radius fluctuates around R ~ 0.5, but some of the subsamples show significantly enhanced correlations in mass (R > 0.85) for planets orbiting cool stars (T_eff < 4400 K) and in density for slightly metal-poor (−0.2 < [Fe/H] < 0) or older (age > 6 Gyr) host stars (see Fig. A.6), suggesting that planet formation mechanisms around cool, old dwarfs conspire to produce more uniform planets. It would be interesting to investigate whether this hypothesis holds in a larger sample of systems with well-determined stellar parameters.

4.2 Light and heavy ‘peas’ in the same ‘pod’

In agreement with previous studies, we find significant intrasys-tem uniformity in planet radii in our sample. This trend could be due to underlying astrophysical reasons, such that similarly sized planets are common outcomes from planet formation processes in a given system (see e.g. Weiss et al. (2022)). However, these trends must be interpreted in the context of selection and observational biases associated with the transit method, which is required to measure planet radii. In addition to the effect of detection bias related to the SNR floor, as pointed out by Zhu (2020), we must remember that our knowledge of the planet systems is inherently incomplete. The transit technique is only sensitive to close-in planets (i.e. ‘inner system architecture’) that transit, which means a large part of the system architecture could be amiss (Zhu & Dong 2021). These missing planets could distort the similarity and spacing trends. An indication of this is the fact that, for mass and density, we find the intrasystem dispersion in systems with more than 3 planets to be greater than that in systems with more than 2 planets.

In contrast with radii, the evidence of system-wide correlation in planet masses is far more tenuous. As discussed in Sect. 4.1.3, there have been few studies on peas in a pod trends in planet mass. While they have found correlations, the results are based on small and possibly biased samples. This is because mass measurements tend to be more difficult than radii, with larger typical errors and more susceptibility to bias. In this work, we used a large sample of the latest available planet mass measurements in multi-planet systems, from both RV and TTV. We find that intrasystem correlation in mass is less prominent than in radius, otherwise proven uniformity in both parameters will be an important constraint for planet formation theories that model planet masses.

One reason for the difference in correlation patterns between planet radii and mass could be purely mathematical. As we imagine the planets as spherical, for a given bulk density, the radius variations are less than those in mass, and radii fluctuate in a much smaller range than masses. For example, in this study’s sample, the radii of planets occupy the range of 0.28–13.2 R_⊕ and masses those of 9.85×10⁻² to 2.16×10³ M_⊕, which leads to great variety in bulk densities. So one can deduce that radii will fluctuate in a much smaller range than masses in the population. As explained in Sect. 2, extreme outliers can affect the behaviour of the Pearson R. Since masses scale as radius cubed, correlations in mass are more prone to being skewed by a few points on the extremity. For example, if we added the massive (and relatively highly correlated) direct imaging planet systems HR 8799, PDS 70, TYC 8988-760-1 to the sample, the Pearson coefficient for masses grows to R = 0.91, whereas it is only R = 0.65 for radius. On the other hand, the intrasystem dispersion and parameter/interior gap calculations, two other methods we used in this work, are not affected by the range of variables due to normalisation. Using all these methods combined, we showed that parameters of the planets in the exoplanetary systems other than radius vary greatly.

A lack of mass correlation despite radius correlation can also imply non-uniformity in bulk density. Indeed, we find only weak evidence for intrasystem correlation in bulk density. Our sample analyses showed that at a certain radius, the planet could have a variety of masses due to widely varying bulk densities. Only a small subsample of planet pairs is highly correlated in both mass and radius (the number of pairs N_P=49, the number of systems N_s=47, the number of planets N_p=97) and even fewer are also similar in density (N_P=35, N_s=34, N_p=70). We acquired these two sub-samples by cross-checking planets belonging to HCR and HCM subsamples in the former, and all three HCR, HCM and HCD subsamples in the latter case.

Within a given system, where formation conditions are expected to be similar (especially in the inner regions probed by the compact multi-planet systems studies), one might expect planets with similar radii to also be similar in composition and interior structure, like Venus and Earth. Therefore, the fact that the exoplanets in the HCR sample tend to have dissimilar bulk densities is surprising. While bulk densities are a proxy for interior structure, their correspondence is not 1:1. To understand the degree of fundamental physical similarities between adjacent planet pairs with superficially similar radii, we locate the HCR sample on planetary interior models of Zeng et al. (2016) in bulk density-mass space (Fig. 3). Planet pairs that are alike in composition and structure should cluster along one specific interior structure contour like “strung beads”. This should be especially prominent for the denser rocky planets, which do not have significant atmospheres to modify their radii. We find that, with few exceptions, the planet pairs straddle across composition contours. Indeed, apparent similarity in planet size does not mean similarity in composition or interior structure.

In their work, Zeng et al. (2019) hypothesised that terrestrial planets, including those that are less dense than the Earth, can be divided by interior composition for two types: predominantly MgSiO₃ with Fe cores (rocky worlds) or predominantly H₂O with MgSiO₃ layers (water-worlds). We include these prescribed interior curves in our analyses (see Fig. 6). We noted, that two adjacent planets can differ in density and mass in such a way, that one planet will belong to rocky worlds and its neighbour will belong to the water worlds. These results contradict Millholland et al. (2017), as they found that astrophysical scatter in general the MRR is dominated by or at least largely linked to different systems with different available planet-forming constituents or, in other words, inter-system differences, rather than intrasystem differences in interiors of planets. On the other hand, there are several notable examples of planets exhibiting a high degree of intrasystem compositional uniformity. One is the TRAPPIST-1 planets, four of which (b, c, f, g) are in our HCR sample. Our results are consistent with Agol et al. (2021), who studied the TRAPPIST-1 system in detail and found that all the planets are consistent with one interior model contour: either one with depleted iron relative to Earth or an Earth-like composition with an enhancement in light elements, such as a surface water layer.

We introduced interior distance as a metric to quantify the compositional difference between similar-sized planets in the same system and suggested I = 0.1 as a threshold to distinguish between planet pairs with similar and dissimilar interior structure. This interpretation compares well with Rodríguez Martínez et al. (2023), who examined core/water mass fractions of planets with M_p < 10 M_⊕ in multi-planet systems around M dwarfs. For example, for the HD 260655 b–c pair, our calculated I = 0.117 would indicate dissimilar structure. This is consistent with Rodríguez Martínez et al. (2023), who determined that the system consists of a rocky and water-rich world. On the other hand, the LTT 1445A b–c pair has I = 0.053 to indicate similar structures, which is consistent with the assessment of Rodríguez Martínez et al. (2023) (rocky planets with comparable core-mass fractions).

The volatile-rich population in a range of Neptune and Uranus-like planets in means of radius and mass is widely represented in our HCR subsample, but their constituents are currently only theorised (Dorn et al. 2015, 2017; Dorn & Lichtenberg 2021) due to our ability to constrain interior structure is directly limited by data uncertainties and lack of such short-period planets in the Solar System. Further constructing a mesh of potential MRR based on interior structures not only for terrestrial planets but for the volatile-rich population by considering their irradiation environment and possible gaseous envelopes, following the work of Zeng et al. (2019), could aid in exploring the interior composition and mass range of a potential atmosphere.

4.3 The extent of our knowledge about exoplanets

Unlike in the Solar System where largely no radius correlation between pairs exists and high intrasystem dispersion in this parameter is present, the radius correlation in exoplanetary systems and intrasystem dispersion persists as an observational pattern. Although many authors report similar findings, the true nature of this pattern was questioned before by Zhu (2020) who attributed it to Kepler observational biases, but also by Murchikova & Tremaine (2020), who used the same sample as Weiss et al. (2018) and showed that similar patterns could arise in a mock universe with unnatural assumptions about the MRR of the planets. Since the radius is the only planetary parameter that shows a consistent correlation pattern in our entire sample and significantly lower intrasystem dispersion compared to other fundamental parameters, it is reasonable to assume that our limited knowledge of the exoplanet occurrence rates, impacted by geometric and detection biases, affects this. Currently, we can determine the radii of transiting planets only, which means we are only sensitive to planets above a certain size threshold, on relatively short orbits and are coplanar with the others in a given system. Further, the Kepler mission had its observational limitations, which we discussed previously. It has been shown that small, Earth-sized planets are very common and greatly outnumber larger planets (e.g., Petigura et al. 2013; Burke et al. 2015; Fulton et al. 2017).

Using Kepler data and simulations, Mulders et al. (2018) predicted that at least half of known exoplanets within the habitable zone (HZ) are accompanied by non-transiting planets at shorter orbital periods. Inversely, it means that at least some of the known systems with close-in planets also are harbouring smaller planets in the HZ and at further separation, transiting or not. The presence of small planets will break the similarity pattern, as was discussed by Zhu (2020). Zhu & Wu (2018) stated that outside the Kepler period limit (~1 yr), cold giant planets preferentially coexist with inner small planets. Recently, Millholland et al. (2022) showed that if the compact multi-planet systems with the prominent peas in a pod patterns in radius and period ratios would have additional transiting planets, consistent with these patterns and located beyond the known planets, the existing observation facilities would already be able to detect them. Thus, if these known compact systems are in fact harbouring additional planets at larger separations, these planets would need to have radii smaller by a factor of ≲0.5 that of super-Earth, which would be Earth-size planets. In recent work, Hatalova et al. (2023) show that typically simulated systems, formed for example around K-star, would harbour several Earth-size planets at larger separations from the star in addition to several close-in Super-Earths. The incompleteness of observed systems could be connected to our findings that systems with three or more planets show large dispersion in mass and density. Therefore, the incompleteness poses a question mark on the apparent similarity in parameters within multi-planet systems.

Finally, the measured parameters of exoplanets depend highly on the precise determination of stellar masses and radii. Gaia DR3 release provides homogeneous data for hundreds of millions of stars but has its shortcomings in characterisation of dim sources (Smart et al. 2021; Riello et al. 2021), i.e. cool M- and K-dwarfs, the stars most interesting in means of finding terrestrial planets. The new Euclid ESA mission (Euclid Collaboration 2022), launched in 2023, will produce a star catalogue based on Gaia data, providing more accurate stellar parameters. Due to the limitation of the methods of exoplanet mass and radius determination, our knowledge of the entire population of exoplanets is incomplete, and it is possible that the apparent trends that we observe in known systems appear only in parts of the intrinsic distributions of fundamental parameters. Future missions such as PLATO (Rauer et al. 2014) along with ongoing missions and growing knowledge about exoplanets’ properties will contribute to overcoming these obstacles.

5 Conclusions

In this paper, we study the ‘peas in a pod’ pattern of radii, masses, densities of adjacent planet pairs, and orbital period ratios of triples belonging to the same system. Our aim is to investigate whether or not the apparent similarity in radius implies similarity in mass, density, and interior structure. We were motivated to revisit these patterns by the fact that the masses of many planets have recently been determined. As in previous works, we find that planet radii within the same system consistently and robustly show similarity as opposed to a random arrangement. We additionally tested radius similarity using the homogenous dataset of stellar radii from Gaia DR3 and come to the same conclusion. This is probably due to one or more of the following reasons: (1) the physical processes creating the planets in the same system with similar radii are probably the same; (2) there are physical constraints that mean that radius will vary less than other planet parameters, (e.g. uniformity of the protoplanetary disc could lead to similar density and interior structure of produced planets); or (3) the current sample is still not reflecting the entire exoplanets population.

Planetary masses within a system are not as uniform as radii and reflect the overall scatter in the MRR for the whole variety of exoplanets, covering a wide range of bulk densities. Using a fine grid of theoretical MRR compositional curves, we investigated the extent to which the peas in a pod pattern propagates into the chemical compositions of planets in compact multiple systems and find that apparently similar-sized planets are often not alike in composition or interior structure. Two adjacent planets similar in radius can belong to the rocky- and water-world categories, respectively. In other words, even within the same system, it cannot be assumed that similarity in radius necessarily implies similarity in mass, density, composition, or internal structure. The principle of separating the data by prescribed interior structures is preferable to separating by the density of a uniform sphere because the constant density would not account for the self-compression of an object. As we gather more information on the detailed nature of planets in these systems, characterisation of this population beyond the superficial parameters of size and mass will help to shed light on the physical origin of these intriguing intrasystem trends, or the lack thereof.

Interestingly, the degree of likeness in adjacent planet masses seems to depend on the stellar properties of the host. For example, small, less metal-rich stars tend to host planets with similar mass and density, which probably suggests there is comparatively less building material around such stars with which to form large planets, meaning that all the formed and detected planets tend to be small and more similarly sized. The similarities in masses and densities within one system could originate from the epoch of planet formation, potentially prescribing a certain amount of mass, and the constituents of a protoplanetary disc, as well as its structure, shape the architecture of the system. Further study of the similarities in mass and density among systems with similar of a host-star parameters using a much larger sample could further our understanding of the planets’ internal structure and the processes of planet formation.

Acknowledgements

We acknowledge financial support from the Research Council of Norway (RCN), through its Centres of Excellence funding scheme, project 223272 (CEED) and project number 332523 (PHAB). We thank Pavel Dubrovine (PHAB), Anders Erikson (German Aerospace Centre, Berlin-Adlershof, Germany), Mathew Domeier (PHAB), and Alexey Pankine (Space Science Institute, USA) for the fruitful discussion. We thank the anonymous reviewer for the thoughtful feedback that helped to improve this manuscript.

Appendix A Supplementary figures

Fig. A.1 Sum of errors normalised to the parameter value in the main sample (left). Black dots represent the normalised sum of errors in radius, white dots are the same for mass. Red dashed and blue dashed lines are for the sum of errors corresponding to the construction of the sample from Otegi et al. (2022) for mass and radius, respectively. Planets with mass determined from RV (white upwards-pointing triangles) and TTVs (black left-pointing triangles) methods in the main sample.

Fig. A.2 Pearson coefficients distributions from the error-accommodation method. Each graph shows R distribution for parameters of adjacent planets, which were tested for linear correlation by running 105 random uniform simulations. From left to right the correlation is calculated for radii, masses, densities of adjacent planets in pairs and period ratios in triples. The blue dash-dotted and black dashed lines indicate mean and median values, respectively. The solid blue line indicates the best-fitted probability density function (PDF) of the normal distribution. The y-axis shows the probability density (the probability per unit on the x-axis).

Fig. A.3 Dispersion in parameters from the error-accommodation method. Each graph shows intrasystem dispersion distribution in the main sample acquired by running 105 random uniform simulations within the error intervals. From left to right the distributions are obtained for radii, masses, densities of adjacent planets in systems with at least two planets and period ratios in systems with at least three planets. The blue dash-dotted and black dashed lines indicate mean and median values, respectively. The solid blue line indicates the best-fitted probability density function (PDF) of the normal distribution. The y-axis shows the probability density (the probability per unit on the x-axis).

Fig. A.4 Expected dispersion from the non-parametric MRR (Ning et al. 2018) compared to dispersion in mass and density from the error-accommodation method. Red solid lines show intrasystem dispersion Derr and the corresponding standard deviation (light-red shadowed regions). From left to right the 105 random distributions from MRR are obtained for masses, and densities of adjacent planets in the main sample, and for masses, and densities of adjacent planets in the HCR sample. The blue dash-dotted and black dashed lines indicate mean and median values, respectively. The solid blue line indicates the best-fitted probability density function (PDF) of the normal distribution. The y-axis shows the probability density (the probability per unit on the x-axis).

Fig. A.5 Moving window test of the “40 data points” moving subsamples for masses (green line), radii (blue line), and densities (grey dashed line) of adjacent planets pairs. R and p-values from the repeated Pearson correlation analyses are plotted against ranges of stellar parameters. The values of the certain stellar parameter on the x-axis are calculated as the mean values for all planets in each subsample.

Fig. A.6 Bootstrap test statistics. Each graph shows the Pearson coefficient distribution for parameter pairs for bootstrapping 40 planet pairs in 105 simulations. From left to right the correlation is calculated for radii, masses, densities of adjacent planets in pairs and period ratios in triples. In the “mass” graph the blue dashed line indicates the “high-correlated” subsample for planets around cool stars (2566-4398 K). In the “density” graph the blue dash-dotted and black dashed lines indicate the “high-correlated” subsamples for planets around stars with metallicity ([Fe/H] ∈[−0.210; −0.010]) and old stars (6.03-12.3 Gyr) respectively.

Fig. A.7 HCR (149 pairs), and stellar 40-pair HCM and HCD subsamples. The HCR sample and comparison to a random bootstrap of 40 pairs (top left).The planet pairs (blue dots) around host stars with similar properties: around cool stars, Teff ∈[2566; 4398] K (masses, bottom left), around those with [Fe/H] ∈[−0.210; −0.010] and old stars, age ∈[6.03; 12.3] Gyr (top and bottom right, densities). The grey dots represent one realisation of the bootstrapped 40-pair samples out of 105 simulations. The red dots correspond to the Solar System planets. The black dashed line is the 1:1 line.

References

All Tables

All Figures

Fig. 1 Masses, radii, densities, and period ratios of adjacent planets. The scatter graph for mass (top left), radius (top right), and density (bottom left) of a planet against the same parameter of the next planet farther from the star is plotted with observational errors indicated by the grey lines. The bottom right panel shows the period ratio between the middle and inner planet against that between the outer and middle planet in triples. The blue and red dot-shaped markers correspond to the sample and Solar System planet pairs, respectively. The black dashed line is the 1:1 line. The observational errors could be smaller than the marker size. In the bottom right corner of each plot, the cumulative distribution function is plotted for the decimal logarithm of each planet parameter ratio (x-axis) versus counts (y-axis). In these subplots, the black vertical dashed line is unity and the red horizontal dashed line is where the cumulative distribution function reaches the 50th percentile.

In the text

	Fig. 2 Intrasystem dispersion in mass, radius, and density. The x-axis shows dispersion in mass for systems in the main sample, the y-axis represents dispersion in radii, and the colour of the dots corresponds to dispersion in density. The red colours indicate similarity in density, and the blue colours indicate dissimilarity. The black dashed line is 1:1.
In the text

Fig. 3 Sample on the mass–density diagram. The black spheres represent the sample’s planets and the red spheres are the Solar System planets, plotted for comparison. Dashed black lines connect a pair of adjacent planets in the “high correlated in radius” subsample. The coloured interior curves (from red to light blue) are taken from Zeng et al. (2016) and correspond to different calculated MRR for several potential 2-layer structures of exoplanets, representing different combinations of constituents of the interior: Fe, MgSiO3, and H2O. The error bars are plotted for planets comprising the “high correlated in radius” subsample and omitted for others for the sake of clarity.

In the text

	Fig. 4 ECDF of the parameter gap ɡ in mass, radius, and density in the HCR subsample. The x-axis shows ɡ in a certain parameter, the y-axis represents the proportion of the gaps in pairs that reached a certain value.
In the text

Fig. 5 Interior difference between a pair of adjacent planets. In this illustrative example, planet P1 lies closest to the blue curve (distance L0) whereas planet P2 is closer to the red curve (distance L1). The blue curve, corresponding to 30% Fe and 70% MgSiO3, is chosen as the reference curve for the pair. This means if P2 had the same bulk composition as P1, P2’s bulk density should be given by the blue curve corresponding to its mass (red circle). Ip1 and Ip2 (shown in the thin red line) are measured as the vertical gap distance between planet P1 and planet P2 to the reference curve at their respective masses. This figure shows LTT 1445A b and c planets, and they are separated by an interior distance I = 0.053.

In the text

Fig. 6 HCR sample on the mass-density diagram. The dots, connected by black dashed lines, are planet pairs that belong to this sample. The coloured dots represent those planets, for which we calculated interior distance I. The red colours indicate a pair of planets at a short interior distance, which means they are probably similar in the interior structure, and blue colours indicate planets that show larger separation in the interior from their neighbour. Light-grey curves represent the MRR interior curves calculated by Zeng et al. (2016). From the baseline of a silicate planet with 100% MgSiO3, the red points in the direction of increasing Fe content, and each curve represents a 10% increase in Fe content up to 100% Fe. Similarly, the blue arrow points in the direction of increasing water content at 10% increments per curve. Solar System planets are added for comparison and plotted as purple stars.

In the text

	Fig. 7 Interior distance in HCR sample systems. Observed I in the x-axis is plotted against expected I from non-parametric MRR (Ning et al. 2018) at the y-axis. The colours of the dots indicate stellar Teff of the system’s host star. The black dashed line is unity.
In the text

	Fig. 8 p-value dynamics for subsets of data points. Coloured lines indicate p-value changes during the moving window test for density depending on host star age (green) or metallicity (grey) and for masses depending on host star Teff (blue). The light-coloured shades under the curves roughly indicate these subsamples in corresponding colours.
In the text

Fig. A.1 Sum of errors normalised to the parameter value in the main sample (left). Black dots represent the normalised sum of errors in radius, white dots are the same for mass. Red dashed and blue dashed lines are for the sum of errors corresponding to the construction of the sample from Otegi et al. (2022) for mass and radius, respectively. Planets with mass determined from RV (white upwards-pointing triangles) and TTVs (black left-pointing triangles) methods in the main sample.

In the text

Fig. A.2 Pearson coefficients distributions from the error-accommodation method. Each graph shows R distribution for parameters of adjacent planets, which were tested for linear correlation by running 105 random uniform simulations. From left to right the correlation is calculated for radii, masses, densities of adjacent planets in pairs and period ratios in triples. The blue dash-dotted and black dashed lines indicate mean and median values, respectively. The solid blue line indicates the best-fitted probability density function (PDF) of the normal distribution. The y-axis shows the probability density (the probability per unit on the x-axis).

In the text

Fig. A.3 Dispersion in parameters from the error-accommodation method. Each graph shows intrasystem dispersion distribution in the main sample acquired by running 105 random uniform simulations within the error intervals. From left to right the distributions are obtained for radii, masses, densities of adjacent planets in systems with at least two planets and period ratios in systems with at least three planets. The blue dash-dotted and black dashed lines indicate mean and median values, respectively. The solid blue line indicates the best-fitted probability density function (PDF) of the normal distribution. The y-axis shows the probability density (the probability per unit on the x-axis).

In the text

Fig. A.4 Expected dispersion from the non-parametric MRR (Ning et al. 2018) compared to dispersion in mass and density from the error-accommodation method. Red solid lines show intrasystem dispersion Derr and the corresponding standard deviation (light-red shadowed regions). From left to right the 105 random distributions from MRR are obtained for masses, and densities of adjacent planets in the main sample, and for masses, and densities of adjacent planets in the HCR sample. The blue dash-dotted and black dashed lines indicate mean and median values, respectively. The solid blue line indicates the best-fitted probability density function (PDF) of the normal distribution. The y-axis shows the probability density (the probability per unit on the x-axis).

In the text

Fig. A.5 Moving window test of the “40 data points” moving subsamples for masses (green line), radii (blue line), and densities (grey dashed line) of adjacent planets pairs. R and p-values from the repeated Pearson correlation analyses are plotted against ranges of stellar parameters. The values of the certain stellar parameter on the x-axis are calculated as the mean values for all planets in each subsample.

In the text

Fig. A.6 Bootstrap test statistics. Each graph shows the Pearson coefficient distribution for parameter pairs for bootstrapping 40 planet pairs in 105 simulations. From left to right the correlation is calculated for radii, masses, densities of adjacent planets in pairs and period ratios in triples. In the “mass” graph the blue dashed line indicates the “high-correlated” subsample for planets around cool stars (2566-4398 K). In the “density” graph the blue dash-dotted and black dashed lines indicate the “high-correlated” subsamples for planets around stars with metallicity ([Fe/H] ∈[−0.210; −0.010]) and old stars (6.03-12.3 Gyr) respectively.

In the text

Fig. A.7 HCR (149 pairs), and stellar 40-pair HCM and HCD subsamples. The HCR sample and comparison to a random bootstrap of 40 pairs (top left).The planet pairs (blue dots) around host stars with similar properties: around cool stars, Teff ∈[2566; 4398] K (masses, bottom left), around those with [Fe/H] ∈[−0.210; −0.010] and old stars, age ∈[6.03; 12.3] Gyr (top and bottom right, densities). The grey dots represent one realisation of the bootstrapped 40-pair samples out of 105 simulations. The red dots correspond to the Solar System planets. The black dashed line is the 1:1 line.

In the text