Sandwiched planet formation: restricting the mass of a middle planet Open Access

ABSTRACT

We conduct gas and dust hydrodynamical simulations of protoplanetary discs with one and two embedded planets to determine the impact that a second planet located further out in the disc has on the potential for subsequent planet formation in the region locally exterior to the inner planet. We show how the presence of a second planet has a strong influence on the collection of solid material near the inner planet, particularly when the outer planet is massive enough to generate a maximum in the disc’s pressure profile. This effect in general acts to reduce the amount of material that can collect in a pressure bump generated by the inner planet. When viewing the inner pressure bump as a location for potential subsequent planet formation of a third planet, we therefore expect that the mass of such a planet will be smaller than it would be in the case without the outer planet, resulting in a small planet being sandwiched between its neighbours – this is in contrast to the expected trend of increasing planet mass with radial distance from the host star. We show that several planetary systems have been observed that do not show this trend but instead have a smaller planet sandwiched in between two more massive planets. We present the idea that such an architecture could be the result of the subsequent formation of a middle planet after its two neighbours formed at some earlier stage.

1 INTRODUCTION

The discovery of the Kepler planets has hugely increased the number of observed exoplanets, and with that our ability to statistically understand the properties of planets and the architectures of planetary systems has grown. One property that is particularly intriguing from a planet formation perspective is the mass ordering of neighbouring planets. Given that the majority of planets known to date have been discovered through the transit technique, most of the data on neighbouring planets in a system appears to be on their radii. Weiss et al. (2018) analysed a set of 909 Kepler planet candidates from the California Kepler Survey from 355 multiplanet systems and found that in general, neighbouring planets seem to be reasonably similar in size with a trend towards larger planets existing at larger orbital separations in |$65{{\ \rm per\ cent}}$| of the cases (also see Lissauer et al. 2011, Ciardi et al. 2013; Chevance, Kruijssen & Longmore 2021). The former trend has been termed Peas in a Pod (see Weiss et al. 2023, for a recent review) and it has been suggested that this phenomenon likely occurs from birth (Chevance, Kruijssen & Longmore 2021). However the Peas in a Pod trend is still under debate given that it is thought that giant planets further out than what can be observed by Kepler are thought to exist with smaller inner planets (Zhu & Wu 2018; Bryan et al. 2019; Herman, Zhu & Wu 2019) as well as the claim that the Peas in a Pod trend may suffer from detection biases (Zhu 2020; Murchikova & Tremaine 2020, though this has also been contested by Weiss & Petigura 2020). This trend may also be affected if small, undetected planets exist in these systems (Zhu 2020).

Understanding the mass trends are more challenging. Millholland, Wang & Laughlin (2017) analyse 55 multiplanet systems consisting of a total of 145 Kepler planets whose masses were determined using Transit Timing Variations (TTVs; Hadden & Lithwick 2017) and showed a similar Peas in the Pod trend for planet mass, whereby neighbouring planets are similar in mass (compared to a randomized distribution of planets) and that in a single system they are ordered in mass so that the outer planets are more massive (also see Bashi & Zucker 2021; Chevance, Kruijssen & Longmore 2021; Goyal & Wang 2022 for non-TTV planets). However, Zhu (2020) also argues that these results may be due to observational biases. More recently, a theoretical study by Mishra et al. (2021) considered whether the Peas in the Pod trends exist in the results of population synthesis models of planet formation and evolution, that are intended to be mock Kepler observations. They find that neighbouring planets do have similar masses and that 55 per cent of neighbouring planets are ordered in mass. Nevertheless this means that a significant proportion of neighbouring planets do not adhere to this mass ordering trend (also see Mishra et al. 2023).

From a planet formation perspective, an increase in planet mass with orbital distance is indeed expected: the isolation mass – which is the amount of solid material that can be fed to the planet before the dust supply is reduced due to the planets’ changing gas surface density (and hence pressure profile) – in a protoplanetary disc increases with radial distance from the star (Lambrechts, Johansen & Morbidelli 2014; Bitsch et al. 2018; Ataiee et al. 2018). In addition, the Hill radius also increases with increasing orbital separation, which means that even for similar mass planets, the amount of material available in an outer planet’s feeding zone is larger due to an increase in the local disc mass, M_{d, local} ∼ ΣR², at larger radii. Moreover, the inside–out planet formation model also leads to an increase in planet mass with distance (Chatterjee & Tan 2014). In this model, dust collects into a ring and grows into a planet at the pressure maximum at the inner edge of the dead zone, after which the planet opens up a gap causing the dead zone to retreat and for the cycle to continue at larger disc radii.

Despite this, there are some ways in which planetary architectures may not follow this typical trend. From an observational perspective, the NASA exoplanet archive¹ shows a substantial fraction of cases where a planet in between two other planets does not follow the trend of increasing mass with increasing orbital distance. We discuss these in more detail in Section 4.2.

Generally sequential planet formation is thought to occur from inside–out due to the time-scales involved in the core accretion process, or naturally through the inside–out planet formation process (Chatterjee & Tan 2014) described above. However in this paper, we explore what happens to the dust flow inwards in a protoplanetary disc when two planets have already formed, and the impact that the outer planet can have on the collection of dust in the ring just exterior to the inner planet. We present the idea of a sandwiched planet formation process – whereby two planets form initially in a disc and then we hypothesize if a third (smaller) planet has some chance of forming in between the two, with a mass that does not follow the expected trend of increasing mass with increasing orbital distance.

How might such a process present itself in the formation process whereby two planets have already formed before the intermediate third one? The inner planet may form as expected from the aforementioned traditional view of sequential planet formation. We may then hypothesize that a second planet may form either preferentially at a special location such as an ice line (Ros & Johansen 2013; Cridland, Pudritz & Birnstiel 2017), or that a second planet may form exterior to the first and migrate outwards (e.g. to a heat transition or ice line, Hasegawa & Pudritz 2011; Speedie et al. 2022). In these cases, the dust flow into and out of the region in between the two planets will be determined by the planetary masses. Alternatively if planets were to form quickly by gravitational instability (which may produce Neptune-mass planets if the discs are magnetized, Deng, Mayer & Helled 2021), we may find ourselves with a similar scenario of two planets having formed with space in between for dust evolution to occur. Finally, whilst not thoroughly explored, the idea of core accretion occurring in a gravitationally unstable disc may also lead to the formation of multiple planets, which could potentially be low in mass depending on the efficiency of this process.

In this work, we make no assumptions about the way in which the two initial planets form but simply assume that they have formed, and observe the resulting dust flow that occurs due to the second planet being present. In Section 2, we present our method. In Sections 3 and 4, we present the results of our simulations and discuss them, respectively. Finally, we present our conclusions in Section 5.

2 METHODS

We perform 2D hydrodynamical simulations of gas and dust in protoplanetary discs with embedded planets, using the fargo3d code (Benítez-Llambay & Masset 2016). Due to the isothermal equation of state, this is a scale-free code whose variables need to be rescaled to give physical quantities. We treat dust species as pressureless fluids, an approximation that is valid for Stokes numbers |$\rm St\lesssim 1$|⁠; dust coagulation and fragmentation as well as backreaction onto the gas are neglected, thus all dust species are independent. Dust diffusion is included with a Schmidt number |$\rm Sc = 1$|⁠. This is largely the same code as that used by Rosotti et al. (2016) and Meru et al. (2019), with the exception of slightly different boundary conditions that we describe below.

The disc is simulated on a N_ϕ × N_R cylindrical grid with N_ϕ = 1024 and N_R = 638, with logarithmically spaced grid cells in the radial direction to ensure square cells across the computational domain. This resolution resolves the Hill radius of the inner and outer planets by 7–8 and 8–10 grid cells, respectively, in both the radial and azimuthal directions. In code units, the inner and outer disc boundaries are located at R = 0.2 and 10, respectively, and we ensure that the choice of the inner boundary does not cause significant artefacts in the inner disc (see Appendix  A for a comparison with an inner boundary of R = 0.5). Closed (antisymmetric) boundary conditions are applied to the gas radial velocity, with gas not being allowed to flow through the inner or outer boundaries. The gas density is extrapolated according to a Keplerian power law at both the inner and outer boundaries. We apply Stockholm wave-killing zones to the gas at both boundaries (de Val-Borro et al. 2006), but with the locations of the wave-killing boundaries as described by McNally et al. (2019). We also ensure that the 2:1 resonant locations for each planet are contained within the undamped region of the disc (Benítez-Llambay et al. 2016). Note that the damping is only applied to the radial velocity of the gas in order to conserve mass and angular momentum as in Dempsey, Muñoz & Lithwick (2020). Inner and outer boundaries are open for the dust (using a symmetric boundary condition for the dust radial velocity), simulating dust flow from the outer disc through the outer boundary and from the simulated region to the inner disc through the inner boundary. The dust density is held fixed at the outer boundary according to the initial dust distribution and extrapolated at the inner boundary according to a Keplerian power law.

In all the simulations, the planets do not feel the gravitational potential of the other planet; and for all simulated systems, there is no overlap of planetary Hill spheres and hence planet–planet interactions are expected to have a limited effect. The planets are introduced into the disc simultaneously and do not accrete material. Furthermore, planet migration is not included, and the planets are held on fixed circular orbits.

2.1 Disc parameters

We use the same disc surface mass density and temperature parameters as that used by Rosotti et al. (2016) and Meru et al. (2019). The disc’s gas distribution is modelled as:

where Σ is the vertically integrated gas density with Σ₀ = 1 × 10⁻³ and R₀ = 1 in code units is a normalization constant. These simulations are scale-free but as an example, for a 1 M_⊙ star and R₀ = 10 au, Σ₀ would be |$89\,\rm g\,cm^{-2}$| and the disc mass simulated would be 0.06 M_⊙. The disc’s aspect ratio varies according to

We choose a fiducial viscosity parameter (Shakura & Sunyaev 1973) of α = 10⁻³ for our main body of simulations. To evaluate the influence of viscosity on our conclusions, additional simulations are performed with α = 10⁻⁴ (see Appendix  B).

We simulate three dust species with logarithmically spaced Stokes numbers between 0.2 to 2 × 10⁻³. The initial dust distribution has the same power-law profile as that of the gas. Given that we simply model how the dust moves in response to the gas without any backreaction, the normalization of the dust density is arbitrary; in what follows we assume a canonical value of the dust-to-gas ratio of 0.01.

2.2 Suite of simulations

We focus on four cases corresponding to inner–outer planet mass combinations where each planet has a weak or relatively strong influence on the the disc’s gas pressure profile; that is, either generates a point of inflection or a pressure maximum. These four cases are illustrated schematically in Fig. 1. The planets are placed at R = 1 and 2 (code units). Assuming a stellar mass of 1 M_⊙, the chosen inner planet masses correspond to 12 and 20 M_⊕ (equivalent to planet-to-star mass ratios of q = 3.6 × 10⁻⁵ and 6 × 10⁻⁵, respectively) resulting in a point of inflection and a maximum in the disc’s pressure profile, respectively. For the outer planet, we simulate masses of 20 and 35 M_⊕ (equivalent to q = 6 × 10⁻⁵ and 1 × 10⁻⁴, respectively). Each of our main simulations were run for 3000 orbits and the low viscosity simulations in Appendix  B were run for 9000 orbits (at the inner planet’s location). For comparison purposes, we also conduct single-planet simulations for each system, that is, simulations of just the inner planet alone. We note that the single 12 and 20 M_⊕ planet cases reproduce the simulations of Rosotti et al. (2016) and give similar results. The two-planet simulations are summarized in Table 1.

In addition to the further simulations with α = 10⁻⁴ mentioned above in Section 2.1, we also conduct additional test simulations with varied resolution and outer boundary locations, and find that this has no effect on the conclusions presented. Changing the inner boundary position did have an effect on the disc profiles, however this was only in the inner disc. For more details, see Appendix  A.

3 RESULTS

Fig. 2 summarizes the results of our simulations by comparing the disc profiles between the single- and two-planet cases. Each row corresponds to the azimuthally averaged profiles for the gas surface mass density, the dust surface mass densities for each dust species, and the quantity dlnP/dlnR, for each of the four cases described above. The latter quantity describes the pressure gradient such that the presence of positive values indicates that a pressure maximum exists, while changes in negative values results in a point of inflection. For these low-mass planets, the effect on the disc profiles are restricted to the regions locally interior or exterior to the planet, with the profiles being close to the initial unperturbed profiles in the outer parts of the disc.

We present the results associated with each of the scenarios referred to in Fig. 1, though Section 3.4 presents the main results associated with the sandwiched planet formation that we present in this paper.

3.1 Small–Small (two points of inflection)

In the single-planet case with only the inner planet, the traffic jam caused by the small influence on the disc pressure profile results in only a modest build up of dust, with the local dust density increasing by a factor of a few compared to the initial value for the largest dust. The outer planet’s small influence on the gas pressure profile also only results in a modest dust enhancement exterior to it. Crucially, the outer planet’s small pressure perturbation has a very small impact on the flow of dust into the inner disc. As a result, the two-planet disc profiles (blue lines in the left column of Fig. 2) in the interplanetary region and close to the inner planet show only small differences from the single inner planet case (red lines).

We note that the reduction in the peak for |$\rm St = 0.02$| is greater than for |$\rm St = 0.2$|⁠. At the time plotted, the amount of dust in the peak has not reached a steady state and continues to increase slowly over time. However, the conclusion that the low-mass outer planet only has a small impact on the dust flow remains unchanged. As seen in Section 3.4, the main focus of this paper is the impact that two massive planets have on the architecture of the planetary system. Therefore, we do not explore this reduction in detail, but note that a future study that quantifies the impact of a small outer planet would be needed to consider this in more detail.

3.2 Big–Small (pressure maximum and point of inflection)

In this case, the inner planet now has a significant impact on the disc profile. Generating a pressure maximum, the inner planet causes the formation of a dust trap exterior to its orbit. With continued inflow from the outer disc, significant dust collection results, with peak increases of two orders of magnitude for the largest dust (see the second plot in the fourth row of Fig. 2). As in Section 3.1, the outer planet’s small influence on the gas pressure profile has little effect on the amount of dust that can drift past it into the interplanetary region, hence the flow of dust into the dust trap is largely unchanged, resulting in a similar dust disc profile in the region exterior to the inner planet to that of the single inner planet case. This is seen by comparing the red lines (inner planet only) to the blue lines (two-planet case) in the second column of Fig. 2.

3.3 Small–Big (point of inflection and pressure maximum)

In this case, it is the outer planet that generates a pressure maximum and the inner planet that creates a traffic jam effect. The resulting dust trap in the outer disc greatly reduces the amount of dust that can drift past the planet and into the interplanetary region, resulting in a large dust depletion compared to the initial profile. For the largest two Stokes numbers, this effect greatly overshadows the modest dust enhancements the inner planet causes in the single-planet case. To see this, compare the red line (inner planet only) with the blue line (both planets) in the third and fourth plots in the third column of Fig. 2. The small dust enhancement exterior to the inner planet is still seen in the two-planet case, but this degree of enhancement is much smaller than the overall level of depletion in the interplanetary region and inner disc. Interestingly, a small enhancement in the surface density of the smallest species of dust (⁠|$\rm St = 0.002$|⁠) is seen in the inner pressure bump (see second plot in the third column of Fig. 2), which may well be due to the outer planet pushing more gas into the interplanetary region.

3.4 Big–Big (two pressure maxima): sandwiched planet formation

As in Section 3.3, the outer planet causes a significant dust depletion in the region interior to its orbit due to its pressure maximum stopping the dust in the outer disc from drifting inwards. This also causes the amount of dust collected in the trap created by the inner planet to decrease compared to the single-planet case for the largest two dust sizes (and especially for the largest dust size; also see Fig. 3). We again see a minor enhancement in the density profile of St = 0.002 dust in the inner dust trap compared to the inner-only case which again, may be attributed to the outer planet pushing more gas (and hence the small-sized dust since it is more coupled to the gas) into the interplanetary region.

We also note that the restricted amount of dust in between the two planets is evident at all times. Fig. 4 shows that the amount of large dust in the pressure maximum increases over time for the single-planet case, whereas with an additional planet the maximum density in the ring levels off.

4 DISCUSSION

4.1 Potential consequences for sequential planet formation

Section 3 shows that when the outer planet is small enough such that the perturbed pressure profile only creates a traffic jam, the impact on the inner dust collection is largely unaffected. As expected, the impact on the inner dust collection happens only when the pressure profile in the outer disc due to the second planet is such that a pressure maximum forms. In this case, if the inner planet is low in mass such that it cannot form a pressure maximum, the inner disc is depleted of its larger sized particles, consistent with past findings that explore dust filtration (Zhu et al. 2012; Rice et al. 2006).

However the more interesting aspect related to planet formation is when both planets are massive enough to create a pressure maximum. In this case, dust cannot be lost into the inner disc, but instead a depletion in the amount of dust in the inner ring occurs. Depending on the amount of dust that is collected in this ring, there could potentially be enough material to form a low-mass planet in between two higher mass planets. This therefore gives a pathway for the formation of planetary systems with planets in a Big–Small–Big mass pattern – we refer to this as ‘sandwiched planet formation’.

This is in contrast to the expectation that higher mass planets should form in the outer disc due to a larger isolation mass, a larger Hill radius, or the inside–out planet formation model (Chatterjee & Tan 2014), but does describe a subset of observed planets whose masses have been determined that do not follow the typical trend of increasing planet mass with radial location.

4.2 Comparison with exoplanet observations

To compare with observations, we only consider planetary systems from the NASA exoplanet archive² that satisfy certain criteria. Any planets that have been flagged as ‘controversial’ have not been included. Planets without measurements of both the orbital period and mass are excluded. To avoid complications arising from degeneracies between planet mass and inclination, we start by only considering systems with at least three transiting planets with known radii which we call the ‘transiting requirement’, and with masses determined to at least 3σ precision. As of 25 April 2023, there are only 17 such planetary systems known. A planet is considered to be small if the mass difference relative to both its neighbours satisfies the following condition,

where M is the planet’s mass. The subscripts i, and i ± 1 refer to the middle planet, and the planets exterior and interior to the middle planet, respectively. The denominator is the lower uncertainty on the mass difference of the middle planet and its neighbour. Equation (3) is calculated using the method described in appendix B in Laursen et al. (2019), which takes into account the asymmetric uncertainties of the planets’ mass measurements. Note that the denominator involves using the lower uncertainty on the mass difference as this is the one that determines whether the middle planet is larger or smaller than the outer planets. Of the 17 planetary systems, 6 of these have a clearly lower mass planet in the middle of the chain surrounded by higher mass planets as shown in Fig. 5. These systems are HD 219134 (Vogt et al. 2015; Gillon et al. 2017), Kepler-80 (MacDonald et al. 2016; Shallue & Vanderburg 2018), Kepler-223 (Mills et al. 2016), TOI-125 (Nielsen et al. 2020), TOI-1246 (Turtelboom et al. 2022), and TRAPPIST-1 (Agol et al. 2021). The blue markers represent the smaller planet sandwiched by two bigger planets. The size of the markers represents the planet’s mass, as does the colour of the two bigger planets. The black hatched markers show the other planets in the system. If the ‘transiting requirement’ is ignored, 19 out of 72 planetary systems contain a lower mass planet sandwiched between two more massive planets.

We now consider the systems with only two detected transiting planets. As observations are biased towards detecting the innermost planets, it is possible that there is an additional undetected outer planet that could potentially result in a sandwiched planetary system. Using the same criteria as above, there are 34 two-planet systems. Of these, six systems are in the Big–Small configuration, that is, where the inner planet is more massive. If this unknown outer planet is larger than the small planet, then these six systems could follow the Big–Small–Big mass ordering. If the ‘transiting requirement’ is relaxed, there are 33 systems out of 193 that could contain a sandwiched smaller planet. Similarly, if we now consider 3+ planet systems where the two outermost detected planets are in the Big–Small configuration, 15 out of 43 such systems could contain a sandwiched smaller planet. Relaxing the ‘transiting requirement’ makes it 34 potential systems out of 98.

The population of known planets is subject to extensive biases, as a result of instrument capabilities, observer choices, and stellar properties. Additionally, multiple detectable planets in-system produces complications, encouraging further observations due to the increased scientific value while making the modelling of such systems more challenging and increasing the amount of data required. Furthermore, the detection method will create biases in the data: observing planets further out by the radial velocity or transit method is hard, and to observe them by direct imaging requires the planet to be bright and massive. Given these biases and the small numbers available, it is premature to draw strong statistical conclusions about the mass ordering of planetary systems. However, the Big–Small–Big mass ordering as considered in this paper is present in a substantial fraction of the current sample, supporting research into the theoretical framework behind such systems. It is also worth noting that our own Solar System also contains sandwiched planets (Mars and Uranus).

4.3 Caveats

In this paper, we introduce the idea of sequential planet formation happening in between two already formed planets, rather than occurring from the inner to the outer disc. In this preliminary study, we focus on the qualitative results but there are many simplifications that have been made that would need to be addressed in future studies to produce more quantitative findings.

Dust coagulation and fragmentation:This is expected to alter the dust size distribution (Drążkowska et al. 2019; Laune et al. 2020; Lau et al. 2022), with the smallest dust in general coagulating to form larger dust particles. In our two-planet simulations, while we see large-scale depletions in the largest dust in the interplanetary regions, the depletion in the smallest dust is far less. This is illustrated in Fig. 6, where the dust density profiles of the three dust species are plotted on the same axis for the 20 and 35 M_⊕ system (and are multiplied by 100 for ease of comparison with the gas density). We expect that the small dust could coagulate to form larger dust that is susceptible to radial drift and trapping in the inner pressure bump. It is therefore expected that accounting for coagulation could result in a larger amount of material in the interplanetary region than that presented here. Further simulations modelling dust coagulation and fragmentation would be needed to check to what extent this process would alter our conclusions.

Dust backreaction: This is expected to result in the deformation of pressure bumps when the dust-to-gas ratio exceeds unity (Taki, Fujimoto & Ida 2016). This deformation acts to smooth out pressure bumps, spreading the dust ring and reducing dust trapping effects (Kanagawa et al. 2018), but does not completely erase the pressure bump (Onishi & Sekiya 2017) and can still lead to planetesimal formation (Carrera et al. 2021). In our simulations, the dust-to-gas ratio does reach the levels required for backreaction to occur. Therefore, we would expect a different disc dust profile and likely a lesser degree of dust accumulation. The general effects of radial drift and dust trapped in the outer pressure bump being unable to get to the inner bump are, however, expected to still hold. As a result, we expect the general conclusion that the amount of dust trapped in the pressure maximum in between the two planets being reduced, and thus the possibility of a low-mass planet being sandwiched in between two more massive planets, to still hold. Again, further simulations would be needed to confirm this.

Planetary dynamics: In this investigation, we adopt highly simplified planetary dynamics. We have not included planet migration nor planet–planet interactions, and it is certainly known that the effect of the former causes complex dust dynamics to take place (Dong et al. 2017; Meru et al. 2019; Wafflard-Fernandez & Baruteau 2020; Chametla & Chrenko 2022). Furthermore, as planets migrate they can pass through strong resonances which may alter their dynamics and stability (e.g. Papaloizou & Szuszkiewicz 2005; Quillen, Bodman & Moore 2013) and it is also important to consider the long-term dynamically stability of the resulting planetary system (after the disc has dispersed). While it is clear that these effects have a significant influence on how a planetary system evolves, it is unclear how these could affect the behaviour we isolate here, if at all. The large degree of complexity these effects introduce means that a careful selection of simulations would be needed to investigate this.

Planet growth: We make no assumptions about the growth of the two planets simulated. That is, we effectively assume that they both form quickly and simultaneously. In reality, a growing planet would pass from the ‘Small’ to the ‘Big’ regime, that is, there would be a period where the growth of two planets in a disc may start off in the regime described by the leftmost panel of Fig. 1, then it would transition to either of the middle panels of Fig. 1, before potentially ending up in that described by the right panel. Therefore, there would be a period where dust is able to flow past the outer planet. The degree to which this would collect in between the two planets would depend on the state of the inner planet at that time; for instance, if the inner planet had grown large enough to form a pressure maximum by this time, we would expect a greater degree of dust collection in the inter planetary region. In addition we note that whilst we perform simulations of planets for 3000 orbits (which is equivalent to ≈9×10⁴ yr at 10 au), there would continue to be a supply of material at least to the outer planet from the outer disc – the longer term impact of this on the growth of the outer planet would need to be considered. To determine the degree to which our conclusions would change when accounting for planetary formation and prior disc evolution, one would need to self-consistently account for when each planet would start to form as well as the formation time-scales at each location.

Disc viscosity: Our main set of simulations consider α = 10⁻³. This is a commonly chosen value for hydrodynamical protoplanetary disc simulations and a reasonable estimate for typical protoplanetary discs. Whilst in principle α may be as large as ∼10⁻², there is growing evidence for lower disc viscosities (α < 10⁻⁴, Fedele et al. 2018; also see Rosotti 2023 for a recent review). In Appendix  B, we also consider the case of α = 10⁻⁴. We draw the same general conclusions from these simulations, that is, that a smaller mass object may form in between two more massive planets, except we find that the planet masses for which this sandwiched effect occurs are smaller.

Stellar and planet masses: We note that here we have given planetary masses when assuming that the central star mass is 1 M_⊙. While most stars are lower in mass, planets have been observed around both more and less massive stars. Since fargo3d is scale-free and works with planet–star mass ratios rather than masses themselves, one can easily rescale for other mass stars, for example, a 20 M_⊕ planet here would correspond to a 40 M_⊕ planet with a 2 M_⊙ host star. However, we also note that for low-mass stars dust is expected to settle to the disc mid-plane faster (Mulders & Dominik 2012) and to move towards the star at a higher radial drift velocity (Pinilla et al. 2013).

Furthermore, the masses we have chosen all correspond to relatively small planets; it is natural to ask how the results would differ in the case of giant, gap-opening planets, such as those in the PDS 70 system (Keppler et al. 2018; Haffert et al. 2019). In this case, the planets would have a much stronger influence on the gas profile, producing much wider and deeper gaps. We note that in such cases there would also be the added effects of gap merging/overlap to consider in addition to the effects of radial drift and dust supply cut-off.

Sandwiched planet, dwarf planet or moon mass object?: It is important to point out that whilst we have presented the idea that a smaller mass planet may form in between two larger planets as a way to potentially describe a Big–Small–Big architecture of planetary systems, we note that we have not quantified the amount of dust that actually accumulates in the interplanetary region. More accurate simulations are needed that take into account some of the above caveats to accurately determine how small or massive the sandwiched object would be.

5 CONCLUSIONS

Planets embedded in protoplanetary discs are known to cause rings of dust, and can have a significant impact on the amount of dust flowing into the inner disc. One key aspect is whether such rings can be the sites of subsequent planet formation. In a two-planet system, the outer planet may influence the amount of solid material in the ring just exterior to the inner planet. When the outer planet is small (i.e. is below the pebble isolation mass at its location), the point of inflection it produces in the disc’s gas pressure profile has only a limited effect on dust inflow. When the outer planet is large enough to generate a pressure maximum the inflow of (larger) dust is halted, as expected. Of particular interest is the case when both planets are massive enough to generate pressure maxima in the disc. In this scenario, the amount of dust trapped exterior to the inner planet is greatly reduced compared to the single-planet case. As a result, if subsequent planet formation is to occur at this location, the architecture of the resulting planetary system is expected to show a small planet sandwiched in between two more massive planets. This is not the architecture that is typically expected through conventional planet formation processes but may describe a subset of planetary systems that are observed to have the architecture where the middle planet is smaller than the inner and outer planets.

ACKNOWLEDGEMENTS

We thank the referee whose insights improved this paper. MP acknowledges support from The Royal Society Enhancement Award. FM acknowledges support from The Royal Society Dorothy Hodgkin Fellowship. SR acknowledges support from a Royal Society Enhancement Award, and funding from the Science & Technology Facilities Council (STFC) through consolidated grant ST/W000857/1. DJA is supported by UK Research and Innovation through the STFC (ST/R00384X/1) and the Engineering & Physical Sciences Research Council (EP/X027562/1). We also thank Richard Booth and Dimitri Veras for helpful discussions. This work was performed using Orac and Avon, the High Performance Computing clusters at the University of Warwick. This work made use of the following software: fargo3d (Benítez-Llambay & Masset 2016), matplotlib (Hunter 2007), and numpy (van der Walt, Colbert & Varoquaux 2011). This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program (DOI 10.26133/NEA12, NASA Exoplanet Science Institute 2020).

DATA AVAILABILITY

The data from the simulations used to create all plots in this article are available on reasonable request to the corresponding author.

Footnotes

References

APPENDIX A: NUMERICAL TESTS

Our simulations are run with an inner boundary at R = 0.2. We note that there is a moderate overdensity of gas, rising sharply above the initial distribution close to the inner boundary. Fig. A1 shows the impact of changing the inner boundary to R = 0.5 for the four fiducial test cases. Moving the boundary outwards leads to an increase in the gas surface mass density and hence pressure gradient (at R ≲ 0.75), which changes the radial drift of larger dust in this region and results in an artificial buildup of dust (not shown here). The effect is minimized by moving the inner boundary as close to the star as is feasible. A logarithmic grid means that, for a given number of azimuthal grid cells, the number of radial grid cells quickly increases upon moving the inner boundary in. For this reason, the computational cost of simulations significantly increases for inner boundaries at smaller radii. Our choice of an inner boundary at R = 0.2 avoids the artificial dust buildup at a reasonable computational cost.

A very small gradual dip in the gas pressure profile is associated with the outer boundary. This effect only becomes noticeable in and near the outer damping zone. With the outer damping zone boundary defined as in McNally et al. (2019), this means that this effect will only be seen in the region local to the outer planet if the outer boundary is located at R ≲ 4. Hence, for the simulations presented in the main text where we set the the outer boundary to R = 10, this outer boundary effect has no influence on our conclusions. Changing the outer boundary also affects the computational cost of a simulation, but to a lesser degree than the inner boundary.

To test the effect of numerical resolution, we perform an additional simulation of the 20 and 35 M_⊕ system with double the resolution (N_ϕ = 2048 and N_R = 1276). The differences in the disc profiles between this case and that in the main text are minimal (see Fig. A2). This increases our confidence in the validity of results of our main suite of simulations. Note that simulations at lower viscosities seem to be less robust to changes in resolution (see Appendix  B, Fig. B2).

When using similar experimental setups, some authors (e.g Rosotti et al. 2016) aim to reduce numerical artefacts by gradually increasing the masses of embedded planets from zero to the full planet mass over the first few tens of orbits. We did not do this for the simulations presented here. An additional test of the fiducial Big–Big case including such a mass tapering (over 20 orbits) showed that the disc profiles with and without tapering become indistinguishable beyond a few hundred orbits.

APPENDIX B: DISC VISCOSITY

Reducing the viscosity parameter, α, has a very large effect on the disc profile both in the the one- and two-planet cases. First, decreasing α reduces the isolation mass (Bitsch et al. 2018; Ataiee et al. 2018). As a result the planet masses used for the α = 1 × 10⁻³ disc are now able to generate pressure maxima. Hence, less massive planets are needed to test the four cases described in Fig. 1. We therefore choose 6 and 8 M_⊕ inner planets and 13 and 16 M_⊕ outer planets, again located at R = 1 and 2, respectively (see Table B1) and run these for 9000 orbits. The qualitative results of these simulations are similar to those described in Section 3, and are displayed in Fig. B1. We note that the sandwiched planet formation effect discussed in Section 4.1 still occurs for lower viscosity discs, albeit for lower planet masses.

We note that these simulations have been run for less than the viscous time. This is purely to illustrate that the sandwiched architecture can potentially occur even for lower mass planets. For a more quantitative study, the simulations would need to be run for longer given that the viscous time is longer in low viscosity discs.

Decreasing α also has an influence on the numerical artefacts discussed in Appendix  A. First, the artefact associated with the inner boundary is enhanced with reduced viscosity. With the inner boundary at R = 0.5 for α = 1 × 10⁻³ a large pressure peak is seen interior to the inner planet. In the case of reduced viscosity, the impact of this on the disc dust profile is much larger (not shown here). Secondly, the effect of changing resolution is also more significant for these reduced α simulations. Whereas the effect of doubling the resolution is minimal when α = 10⁻³ (see Fig. A2), this is not the case for α = 10⁻⁴; with this reduced viscosity, though most qualitative features are preserved, Fig. B2 shows how there are significant quantitative effects when the resolution is changed. We do not closely consider the quantitative results of simulations in this paper. Should later work wish to do so, the choice of resolution would need to be carefully considered.