© ESO 2022

1 Introduction

During the past decade, the search of planets around very low mass stars (VLMSs) has increased significantly even with the confirmation of Earth-like planets around brown dwarfs (BDs) through transit and radial velocity observations mainly from the Keppler/K2 missions, the HARPS (High Accuracy Radial velocity Planet Searcher) and CARMENES (Calar Alto high-Resolution search for M dwarfs with Exoearths with Near-infrared and optical Échelle Spectrographs) spectrographs, and the Trappist (Transiting Planets and Planetesimals Small Telescope) telescope (Muirhead et al. 2012; Gillon et al. 2017; Astudillo-Defru et al. 2017; Grimm et al. 2018; Crossfield et al. 2019; Zechmeister et al. 2019; Dreizler et al. 2020; Sabotta et al. 2021). The close-in rocky planets that are hosted by such low-mass objects inside the habitable zones are ideal targets for the search of life in the solar neighborhood. Therefore, future space missions such as PLATO (Planetary Transits and Oscillations of stars) and JWST (James Webb Space Telescope) will be able to detect and atmospherically characterize Earth-like planets in the habitable zones around M dwarfs (Rauer et al. 2014). From a theoretical perspective, only few models have been developed to study the rocky planet formation around VLMSs. Initially Payne & Lodato (2007); Raymond et al. (2007); Ciesla et al. (2015); Liu et al. (2020) predicted planetary systems with more compact orbits than the systems that revolve around more massive stars. Moreover, they found that the mass of the planets increases with the mass of the host. Further results showed the relevance of tidal and general relativistic effects for the orbital changes of these compact systems. The contraction and spin-up of VLMSs during the pre-main sequence phase also allows the planet population to follow different dynamical paths (e.g., Barnes et al. 2010; Heller et al. 2010; Bolmont et al. 2011, 2013; Sánchez et al. 2020). Recently Coleman et al. (2019) studied the formation of the Trappist-1 system by incorporating the migration and orbital damping of the planetary embryos that is caused by torques exerted by the gas disk (e.g., Papaloizou & Larwood 2000; Tanaka & Ward 2004; Paardekooper et al. 2010, 2011; Cresswell & Nelson 2008; Ida et al. 2020) and showed its relevance in the early dynamic evolution of the system.

The standard initial mass function predicts that VLMSs and BDs are the most abundant objects in the solar neighborhood. The planetary formation around these objects is therefore essential for estimating the probability of Earth-like planets in habitable zones. Their proximity and number makes them relevant targets in the search for potentially habitable planets (e.g. Kasting et al. 1993; Selsis et al. 2007; Kopparapu et al. 2013; Barnes et al. 2013).

We study the rocky planet formation around a star with a mass of 0.08 M_⨀, which is close to the substellar mass limit. We conducted N-body simulations considering an embryo population orbiting the star. We included tidal and general relativistic effects and the contraction and spin-up of the star during the first 100 Myr of its evolution. We incorporated the gas-disk interactions onto the embryos during the first 10 Myr by assuming a standard disk model and using two different prescriptions for the corresponding torques: the classic formulas from Cresswell & Nelson (2008), and the recent results from Ida et al. (2020). Our aim is to make a comparative analysis of the two prescriptions for the gas-disk torques through a dynamical analysis along the gas and post-gas stages and in the final architectures of the resulting planetary systems. We also compare these systems with observed counterparts.

In Sect. 2, we describe the standard disk model. In Sect. 3 we explain the implementation of the two prescriptions for the gas torques and the set of test simulations we conducted to guarantee the agreement between the numerical and analytical models. In Sect. 4 we characterize the initial parameters of our N-body simulations. In Sect. 5, we develop a detailed analysis of the planet dynamics and architectures during the gas and post-gas stages. Finally, Sect. 6 summarizes our conclusions.

2 Disk model

We adopted the disk model from Ida et al. (2016) in which the structure is described by the gas surface density profile Σ_g, the disk temperature T_g, and the gas-disk aspect ratio h_g = H_g/r, where r is the radial coordinate in the midplane of the disk, and H_g is the gas scale height, which depends on the heating process. The model includes two different dominant heating mechanisms: the internal viscous dissipation for the inner disk and the irradiation from the central star for the outer disk. For the inner region of the disk Σ_g, T_g and h_g profiles are given by: $$\begin{array}{ll}{\displaystyle {\sum }_{\text{g,vis}}{=}^{-4/5}}\hfill & 2100{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot {\text{yr}}^{-1}}}\right)}^{3/5}{\left(\frac{M_{*}}{M_{\odot }}\right)}^{1/5}{\left(\frac{{\alpha }_{\text{g}}}{{10}^{-3}}\right)}^{-4/5}\hfill \\ \hfill & {\left(\frac{r}{\text{au}}\right)}^{-3/5}\text{g}{\text{cm}}^{-2}\hfill \end{array}$$(1) $$\begin{array}{ll}{T}_{\text{g,vis}}=\hfill & 200{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{2/5}{\left(\frac{M_{*}}{M_{\odot }}\right)}^{3/10}{\left(\frac{{\alpha }_{\text{g}}}{{10}^{-3}}\right)}^{-1/5}\hfill \\ \hfill & {\left(\frac{r}{\text{au}}\right)}^{-9/10}\text{K}\hfill \end{array}$$(2) $$\begin{array}{ll}{h}_{\text{g,vis}}=\hfill & 0.027{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{1/5}{\left(\frac{M_{*}}{M_{\odot }}\right)}^{-7/20}{\left(\frac{{\alpha }_{\text{g}}}{{10}^{-3}}\right)}^{-1/10}\hfill \\ \hfill & {\left(\frac{r}{\text{au}}\right)}^{1/20}\hfill \end{array}$$(3)

where M_⋆ is the mass of the central object, M_g is the gas accretion rate and α_g is the viscous coefficient related to the viscosity v = α_gc_sT_gH_g, where c_s is the sound speed at the temperature of the disk midplane at a given radial distance (Shakura & Sunyaev 1973). We assumed that the inner disk is optically thick with an average opacity K = 1 cm² g⁻¹.

In order to smooth Σ_g,vis at the inner edge of the disk we multiplied it with the term tanh[(r − r₀)/(r₀ho)], where r₀ and h₀ are the radius and aspect radius at the inner edge respectively, as was suggested by Cossou et al. (2014), Matsumura et al. (2017), and Brasser et al. (2018). Thus, the surface density profile in the viscous region can be expressed as follows: $${\displaystyle {\sum }_{\text{g,vis}}={\displaystyle {\sum }_{\text{g,vis}}\tanh \left[\left(r-r_0\right)/\left(r_0{h}_0\right)\right]}.}$$(4)

For the outer region of the disk the corresponding profiles are given by $$\begin{array}{c}{\displaystyle {\sum }_{\text{g,irr}}=2700\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right){\left(\frac{M_{*}}{M_{\odot }}\right)}^{9/14}{\left(\frac{{\alpha }_{\text{g}}}{{10}^{-3}}\right)}^{-1}}\\ {\left(\frac{L_{*}}{L_{\odot }}\right)}^{-2/7}{\left(\frac{r}{\text{au}}\right)}^{-15/14}\text{g}{\text{cm}}^{-2}\end{array}$$(5) $$T_{\text{g,irr}}=150{\left(\frac{M_{*}}{M_{\odot }}\right)}^{-1/7}{\left(\frac{L_{*}}{{\text{L}}_{{}_{\odot }}}\right)}^{2/7}{\left(\frac{r}{\text{au}}\right)}^{-3/7}\text{K}$$(6) $$h_{\text{g,irr}}=0.024{\left(\frac{M_{*}}{M_{\odot }}\right)}^{-4/7}{\left(\frac{L_{*}}{{\text{L}}_{\odot }}\right)}^{1/7}{\left(\frac{r}{\text{au}}\right)}^{2/7}$$(7)

where L_⋆ is the luminosity of the central object, and for simplicity, the disk was assumed to be vertical optically thin.

The boundary r_tran that separates the viscous region from the irradiated region is given by $$\begin{array}{ll}{r}_{\text{tran}}=\hfill & 1.8{\left(\frac{L_{*}}{L_{\odot }}\right)}^{-20/23}{\left(\frac{M_{*}}{M_{\odot }}\right)}^{31/33}{\left(\frac{\alpha }{{10}^{-3}}\right)}^{-14/33}\hfill \\ \hfill & {\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot {\text{yr}}^{-1}}}\right)}^{28/33}\text{au}\hfill \end{array}$$(8)

All the initial parameters were set at 1 Myr. We considered a central object with a mass M_⋆ = 0.08 M_⊙, which is close to the substellar mass limit, and an evolving luminosity as predicted by the evolutionary models from Baraffe et al. (2015).

The inner disk edge was set at r₀ = 0.015 au ≃ 3R_⋆‚0 with R_⋆,0 the initial radius of the central object from the Baraffe et al. (2015) models and a corresponding aspect ratio h0 = 0.03, calculated using h_g,vis as r_tran = 0.086 au at 1 Myr. For the viscous coefficient, we used α_g = 0.001, which is a commonly used value for disks around BDs (Adame et al. 2011). The gas accretion rate M_g was obtained from the fit proposed by Manara et al. (2012) which is based on a large sample of accreting stars and BDs in the Orion nebular cluster. This approach was also adopted in Liu et al. (2020) to treat planet formation down to the substellar mass regime. Then, the evolution of M_g in time t is given by $$\begin{array}{ll}\log \left(\frac{{\dot{M}}_{\text{g}}}{M_{\odot }{\text{yr}}^{-1}}\right)=\hfill & -5.12-0.46\log \left(\frac{t}{\text{yr}}\right)-5.75\log \left(\frac{M_{*}}{M_{\odot }}\right)\hfill \\ \hfill & +1.17\log \left(\frac{t}{\text{yr}}\right)\log \left(\frac{M_{*}}{M_{\odot }}\right).\hfill \end{array}$$(9)

Several observational results for young stellar populations show that for a given age, the number fraction of BDs harboring disks is higher than for low-mass stars (LMSs), which in turn is higher than the faction for more massive stars. For instance, at ages of ~7 Myr, ~30% of the BDs still retain their disks. For low-mass stars, this fraction is 5% (e.g. Damjanov et al. 2007; Bayo et al. 2012; Riaz et al. 2012; Dawson et al. 2013; Downes et al. 2015; Manzo-Martínez et al. 2020). These results suggest that VLMSs and BDs could retain their primordial disks for longer times than more massive stars; the times may be up to several tens of million years. Then, we set a disk lifetime of 10 Myr for our disk model. Although photoevaporation might cause some gas dispersion, the dominant dispersion process is the accretion onto the central object. As explained by Stamatellos & Herczeg (2015), current uncertainties of how UV and X-ray emissions from BDs would affect their disks prevents an accurate estimate of the contribution of photoevaporation to the disk dispersion. Although is well known that this effect could occur in disks surrounding VLMSs (e.g. Alexander et al. 2006), their photoevaporation rates are estimated to be as low as ~10−¹¹ M_⊙ yr−¹ (Herczeg 2007). Even in the rough approximation of a constant photoevaporation rate of 10−¹¹ M_⊙ yr−¹ during the first 10 Myr, the accretion onto the central object is always higher than this value in this period of time. We therefore ignored the effect of disk photoevaporation produced by radiation from the central object.

3 Treatment of the disk-embryo interaction

In this section, we describe and compare two different analytical prescriptions for the torques exerted by the gas-disk onto the embedded embryos. We compare them using test simulations.

3.1 IDA20 and CN08 torque prescriptions

We analyzed the migration, and the eccentricity and inclination decay that embryos experience through their interaction with the gas component of the disk following the new prescription from Ida et al. (2020, from now on IDA20) and the classic formulas from Cresswell & Nelson (2008, from now on CN08). In both cases, the torques that the gas exerts on the embryos were computed using the nonisothermal disk model from Paardekooper et al. (2010, 2011), which includes thermal and viscous diffusion. When a gravitational smoothing length of b = 0.4h_g is assumed, the total torque over each embryo is given by $${\Gamma }_{\text{total}}={\Delta }_{\text{L}}{\Gamma }_{\text{L}}+{\Delta }_{\text{C}}{\Gamma }_{\text{C}},$$(10)

where ∆_L and ∆_C are the reduction factors for noncircular or coplanar planetary orbits. The reduction factors differ in the CN08 and IDA20 prescriptions, as we discuss in Sects. 3.2 and 3.3. The factors Γ_L and Γ_C represent the Lindblad and coro-tation torques for a circular and coplanar motion, respectively, given by $${\Gamma }_{\text{L}}=\left(-2.5-1.7\beta +0.\delta \right)\frac{{\Gamma }_0}{\gamma \text{eff}}$$(11)

and $$\begin{array}{ll}{\Gamma }_{\text{C}}=\hfill & {\Gamma }_{\text{c,hs,baro}}F\left(P_v\right)G\left(P_v\right)+\left(1-K\left(p_v\right)\right){\Gamma }_{\text{c,lin,baro}}\hfill \\ \hfill & +{\Gamma }_{\text{c,hs,ent}}F\left(P_v\right)F\left(P_x\right)\sqrt{G\left(P_v\right)G\left(P_x\right)}\hfill \\ \hfill & +\sqrt{\left(1-K\left(P_v\right)\right)\left(1-K\left(P_x\right)\right){\Gamma }_{\text{c,lin,ent}}.}\hfill \end{array}$$(12)

Where Γ_c,hs,baro and Γ_c,lin,baro are barotropic terms related to the horseshoe drag and the linear corotation torque, respectively, and Γ_c,hs,baro and Γ_c,lin,baro are their corresponding nonbarotropic entropy counterparts. These terms are given by $${\Gamma }_{\text{c,hs,baro}}=1.1\left(1.5-\delta \right)\frac{{\Gamma }_0}{\gamma \text{eff}},$$(13) $${\Gamma }_{\text{c,lin,baro}}=0.7\left(1.5-\delta \right)\frac{{\Gamma }_0}{\gamma \text{eff}},$$(14) $${\Gamma }_{\text{c,lin,ent}}=7.9ϵ\frac{{\Gamma }_0}{{\gamma }^2{}_{\text{eff}}},$$(15) $${\Gamma }_{\text{c,lin,ent}}=\left(2.2-\frac{1.4}{\gamma e\text{ff}}\right)ϵ\frac{{\Gamma }_0}{\gamma e\text{ff}},$$(16)

where the scaling torque is ${\text{Γ}}_{\text{0}}={\left(M_{\text{p}}/M_{*}\right)}^2{\displaystyle {\sum }_{\text{g}}{r}^4{h}_{\text{g}}^{-2}{\Omega }_{\text{k}}^2}$ with the angular Keplerian velocity Ω_k. The negative of the entropy slope is ϵ = β−(γ - 1)δ, with δ = −d ln Σ_g/d ln r,β = −d ln T_g/d ln r and γ = 1.4 the adiabatic index. In the viscous region, the parameter δ was calculated by considering the surface density profile for the viscous region from Eq. (4). The multiplied factor added to the original surface density viscous profile becomes relevant only for a region close to the inner edge of the disk up to r ~ 0.016 au. On the other hand, we note that in the irradiated region of the disk, δ and β take values that make ϵ = 0. In the irradiated region, the entropy contributions to the corotation torques are therefore null. The effective γ_eff is given by $$\gamma \text{eff=}\frac{2Q\gamma }{\gamma Q+0.5\sqrt{2\sqrt{{\left({\gamma }^2{Q}^2+1\right)}^2-16Q^2\left(\gamma -1\right)}+2{\gamma }^2{Q}^2-2}}$$(17)

which is related to thermal diffusion by the coefficients $Q=2\chi /3h_{\text{g}}^3{r}^2{\Omega }_{\text{k}}$ and$\chi =16\gamma \left(\gamma -1\right)\sigma T_{\text{g}}^4/\left[3\kappa {\left({\rho }_{\text{g}}{h}_{\text{g}}r{\Omega }_{\text{k}}\right)}^2\right]$ with σ the Stefan-Boltzmann constant, K the gas opacity, and ρ_g the volumetric gas density ${\rho }_{\text{g}}={\displaystyle {\sum }_{\text{g}}/\left(H_{\text{g}}\sqrt{2\pi }\right)}$.

Additionally, the functions F(p), G(p), and K(p) from Eq. (12) are given by $$F\left(P\right)=\frac{1}{1+{\left(\frac{P}{1.3}\right)}^2},$$(18) $$G\left(P\right)=\{\begin{array}{ll}\frac{16}{25}{\left(\frac{45\pi }{8}\right)}^{3/4}{P}^{3/2}\hfill & \text{if}p<\sqrt{\frac{8}{45\pi }}\hfill \\ 1-\frac{9}{25}{\left(\frac{8}{45\pi }\right)}^{4/3}{p}^{-8/3}\hfill & \text{if}p\ge \sqrt{\frac{8}{45\pi }}\hfill \end{array}$$(19) $$k\left(p\right)=\{\begin{array}{ll}\frac{16}{25}{\left(\frac{45\pi }{8}\right)}^{3/2}{P}^{3/2}\hfill & \text{if}p\text{<}\sqrt{\frac{28}{45\pi }}\hfill \\ 1-\frac{9}{25}{\left(\frac{28}{45\pi }\right)}^{4/3}{p}^{-8/3}\hfill & \text{if}p\ge \sqrt{\frac{28}{45\pi }}\hfill \end{array}.$$(20)

The functions are evaluated in p, which takes the form of p_v, the saturation parameter associated with viscosity, or p_χ, the saturation parameter related to thermal diffusion. Both parameters are given by $$p_v=\frac{2}{3}\sqrt{\frac{r^2{\Omega }_{\text{k}}{x}_{\text{s}}^3}{2\pi \nu }},$$(21)

where $x_s=\left(1.1/{\gamma }_{\text{eff}}^{0.25}\right)\sqrt{M_p/\left(M_{*}{h}_{\text{g}}\right)}$ is the nondimensional half-width of the horseshoe region, and $$P_{\chi }=\sqrt{\frac{r^2{\Omega }_{\text{k}}{x}_{\text{s}}^3}{2\pi \chi }.}$$(22)

Following Coleman & Nelson (2014); Cossou et al. (2014); Izidoro et al. (2017); Carrera et al. (2018); Raymond et al. (2018), we adopted a unique value for α. Recent models (e.g., Ida et al. 2018; Matsumura et al. 2021) include both a widen-driven disk accretion α for the inner disk profiles and a turbulent α to calculate the local planet-disk interactions. The same values we adopted for the accretion and turbulent α = 10⁻³ also agree with one of the values proposed by Matsumura et al. (2021). An exploration of different values for the turbulent α is beyond the scope of this work.

All previous formulas that were used to calculate Lindblad and corotation torques were evaluated at the semimajor axis a of the orbit of the embryo. The eccentricity e and inclination i values of the orbit of the embryo were included when the reduction factors were calculated.

3.2 Reduction factors and acceleration from IDA20

The new prescription from IDA20 studies the gravitational interactions between the gas and the embryos on the basis of dynamical friction resulting in reduction factors given by $${\Delta }_{\text{L}}={\left(1+\frac{C_{\text{p}}}{C_{\text{M}}}\sqrt{e_{\text{rat}}^2+i_{\text{rat}}^2}\right)}^{-1},$$(23) $${\Delta }_{\text{C}}=\exp \left(-\frac{\sqrt{e^2+i^2}}{e_{\text{f}}}\right),$$(24)

where C_P = 2.5 − 0.1δ + 1.7β, C_M = 6(2δ − ß + 2), e_rat = e/h_g, i_rat = i/h_g and e_f = 0.5h_g + 0.0l.

In cylindrical coordinates (r, θ, z), the equations of motion for an embryo with a velocity u = (v_r, v_θ, v_z) are given by $$\frac{\text{d}\upsilon }{\text{d}t}=-\frac{\upsilon \text{r}}{t_{\text{e}}}{ê}_{\text{r}}-\frac{\left({\upsilon }_{\theta }-{\upsilon }_{\text{k}}\right)}{t_{\text{e}}}{ê}_{\theta }-\frac{{\upsilon }_{\text{k}}}{2t_{\text{a}}}{ê}_{\theta }-\frac{{\upsilon }_Z}{t_{\text{i}}}{ê}_Z,$$(25)

where ${\widehat{e}}_{\text{r}},{\widehat{e}}_{\theta }$, and ${\widehat{e}}_{\text{z}}$ are versors in the respective directions. The gas velocity is given by υ_g = (0, (1 - η)υ_k, 0), with υ_k the Keplerian velocity, $\eta ~h_{\text{g}}^2$ and η⁻¹t_e = t_a. The variables t_a, t_e, and t_i represent the damping timescales of the semimajor axis a, the eccentricity e, and the inclination i of the orbit of the embryo, respectively. Considering that the embryo migration due to its interaction with the gas is nonisothermal and assuming the condition i < h_g, we can express the damping timescales as $$t_{\text{a}}=-\frac{t_{\text{wave}}}{2h_{\text{g}}^2}\frac{{\Gamma }_0}{{\Gamma }_{\text{total}}},$$(26) $$t_{\text{e}}=\frac{t_{\text{wave}}}{0.78}\left[1+\frac{1}{15}{\left(e_{\text{rat}}^2+i_{\text{rat}}^2\right)}^{3/2}\right],$$(27) $$t_{\text{i}}=\frac{t_{\text{wave}}}{0.544}\left[1+\frac{1}{21.5}{\left(e_{\text{rat}}^2+i_{\text{rat}}^2\right)}^{3/2}\right],$$(28)

where $$t_{\text{wave}}=\left(\frac{M_{*}}{M_{\text{p}}}\right)\left(\frac{M_{*}}{{\displaystyle \sum {}_{\text{g}}{r}^2}}\right)h_{\text{g}}^4{\Omega }_{\text{k}}^{-1},$$(29)

is the timescale from Papaloizou & Larwood (2000) and Tanaka & Ward (2004), in which all physical parameters are evaluated at the semimajor axis of the orbit.

3.3 Reduction factors and acceleration from CN08

The CN08 prescriptions were obtained by fitting analytic formulas to hydrodynamic simulations of planets with eccentric and inclined orbits embedded in the gas disk. For the hydrodynamic simulations they treated the disk as a viscous fluid and preserved its mass by applying reflecting boundary conditions at the inner and outer boundaries. In this prescription the reduction factors are given by $${\Delta }_{\text{L}}={\left[P_{\text{e}}+\frac{P_{\text{e}}}{\left|P_{\text{e}}\right|}\left\{0.07i_{\text{rat}}+0.085i_{\text{rat}}^4-0.08e_{\text{rat}}{i}_{\text{rat}}^2\right\}\right]}^{-1},$$(30)

where $$P_{\text{e}}=\frac{1+{\left(\frac{e}{2.25h_{\text{g}}}\right)}^{6/5}+{\left(\frac{e}{2.84h_{\text{g}}}\right)}^6}{1-\left(\frac{e}{2.02h_{\text{g}}}\right)}$$(31) $${\Delta }_{\text{C}}=\exp \left(-\frac{e}{e_{\text{f}}}\right)\left[1-\tanh \left(i_{\text{rat}}\right)\right].$$(32)

From the comparison of the hydrodynamic simulations with N-body simulations, CN08 found that the acceleration for the embryos is given by $$\frac{\text{d}\upsilon }{\text{d}t}=-\frac{\upsilon }{t_{\text{m}}}-2\frac{\left(\upsilon .r\right)}{r^2{t}_{\text{e}}}-2\frac{{\upsilon }_{\text{z}}}{t_{\text{i}}}k,$$(33)

where r and u are the position and velocity vectors of the embryo in Cartesian coordinates and k is the versor in the z-direction. Additionally, the migration, eccentricity and inclination damping timescales are given by $$t_{\text{m}}=-L\frac{{\Gamma }_0}{{\Gamma }_{\text{total}}},$$(34) $$t_{\text{i}}=\frac{t_{\text{wave}}}{0.78}\left(1-0.14e_{\text{rat}}^2+0.06e_{\text{rat}}^3+0.18e_{\text{rat}}{i}_{\text{rat}}^2\right),$$(35) $$t_{\text{i}}=\frac{t_{\text{wave}}}{0.544}\left(1-0.3i_{\text{rat}}^2+0.24i_{\text{rat}}^3+0.18e_{\text{rat}}^2{i}_{\text{rat}}\right),$$(36)

where $L=M_{\text{p}}\sqrt{G{M}_{*}a\left(1-e^2\right)}$ is the orbital angular momentum of the embryo, M_p is its mass, G is the gravitational constant, t_wave is given by Eq. (29) and the total torque Γ_total includes the reduction factors given by Eq. (30) and Eq. (32).

It is important to point out that in CN08, the associated accelerations are related with t_m, t_e, and t_i, while in IDA20, they are associated with ta, t_e, and t_i. We clarify that the timescales are related by $$t_{\text{m}}={\left[\frac{1}{2}{t}_a^{-1}-\frac{e^2}{1-e^2}{t}_e^{-1}-{\text{i}}^2{t}_{\text{i}}^{-1}\right]}^{-1}.$$(37)

The migration timescale tm does not represent the actual evolution of the semimajor axis. Thus a change in the orbital angular momentum L could happen both with a change in semimajor axis or with a change in the orbital eccentricity of a planet. This indicates that an inward migration is not always associated with a decay in semimajor axis if the orbit is noncircular.

Fig. 1 Maps of the total normalized torque Gtotal/G0 exerted by the gas onto embryos in a 1 Myr old disk for: coplanar and circular orbits (top panels), orbits with e = 0.1 following the prescriptions from IDA20 (middle panels), and orbits with e = 0.1 following the prescriptions from CN08 (bottom panels). The right panels show zoom-ins of the inner part of the disk.

3.4 Comparison of type I migrations from IDA20 and CN08

The migration and orbital decay experienced by the embryos due to the gas were modeled by IDA20 and CN08. Both teams considered the same estimates for the Lindblad and corotation torques, but differed in the reduction factors when considering a noncircular or coplanar orbit, the damping timescales for a, e, and i and the accelerations terms. We evaluated the consistency of the two prescriptions by comparing the corresponding total torques for an embryo in two different scenarios: when its orbit is circular and coplanar with the disk midplane, and when it is not. In both cases we used the disk model described in Sect. 2.

Figure 1 shows the maps of Γ_total normalized by Γ₀ for embryos with masses M_emb and semimajor axis a within the ranges $M<M_{emb<10M_{\oplus }}$ and 0.015 < a/au < 5, and setting the disk parameters for an age of 1 Myr. The left panel shows the complete range of a, and the right panel focuses on the inner part 0.015 < a/au < 0.02.

The top panels show the circular and coplanar case, in which the reduction factors ∆_L and ∆_C equal unity and both prescriptions result in the same Γ_total. In this case, the total torque becomes positive for embryos closer to the inner edge and for masses M_emb < 3 M_⊕, while Γ_total ~ 0 only for embryos at some particular distances close to the inner edge of the disk, which depends on the mass of the embryo. For the remaining combinations of a and M_emb the total torque Γ_total remains negative. However, it became even more negative for r > r_tran and M_emb > 0.3 M_⊕. The middle panels shows the results considering coplanar orbits with e = 0.1 and Γ_total calculated using ∆_L and ∆_C from the IDA20 prescription. In this case Γ_total is always negative (even more negative for r > r_tran), except for a narrow region at a ~ 0.016 au where it becomes positive even close to zero. This small region of positive torque values is produced by the maximum of the density profile of the gas and the inner negative values for the torques are due to the trap considered close to the inner edge. As we described in Sect. 3.2, ∆_L under the IDA20 prescription includes the parameters δ and β, which are related to the disk density profiles and lead to differences when either the eccentricity or the inclination of the planet is non-negligible. The bottom panel shows the Γ_total calculated for a coplanar orbit with e = 0.1 using ∆_L and ∆_C from CN08. In this case, Γ_total changes from positive values in the inner region of the disk to negative values in the outer disk. The values are almost zero around a ~ 0.4 au. The region closest to the inner edge has the highest positive torque values, and it decreases with r even more after r > r_tran. The increase in e produces positive values of Γ_total for a wide range of radial distances under the CN08 treatment. However, there is no particular treatment for the trap close to the inner edge because the factors ∆_L and ∆_C are not related with the disk profiles, just with the value of the eccentricity and inclination.

We found that the addition of an orbital inclination i ≤ h_g results in similar torque patterns, but they are closer to zero than those in the previous analysis. As i increases, Γ_total becomes even closer to zero, as expected because the embryo orbits are further from the midplane of the disk most of the time.

Finally, we explored how the Γ_total maps change throughout the gas-disk lifetime. We found that at t ~ 3 Myr and later, the Γ_total remained negative for the complete ranges of M_emb and a for both circular and coplanar orbits and also for orbits with e = 0.1 if it follows the IDA20 prescription. On the other hand, for an orbit with e = 0.1 under the CN08 treatment, an analogous pattern in the Γ_total map remained, but the inner positive torques decreased in absolute value and the outer positive torques were extended up to higher values of a until they reached a ~ 0.8 ua at t = 10 Myr.

To summarize, the Γ_total maps from the two prescriptions differ significantly in the case of noncircular orbits. We therefore searched for the differences in the evolution history paths of the embryos and the final architectures of the planetary systems when we applied one or the other prescription.

3.5 Test of numerical simulations

We modified the well-known MERCURY code (Chambers 1999) by adding the torque prescriptions from IDA20 and CN08 and the disk model discussed in Sect. 2 to our previous modification including tidal and general relativistic effects (Sánchez et al. 2020). Then, we tested the agreement within the external forces and the damping timescales by analyzing the orbital evolution from a set of N-body simulations of a planet with different initial parameters that followed either IDA20 or CN08.

Figure 2 shows the first 1 Myr evolution of a, e, and the absolute values of t_a, t_e and t_i for an inner planet initially located in an orbit with a = 0.016 au and e = 0.1 and an outer planet initially located in an orbit with a = 0.1 au and e = 0.5. In both cases the planets have masses M_P = M_Mars and the simulations were performed twice following the prescriptions from IDA20 and CN08. We indicate the separation between the subsonic (e < h_g) and supersonic (e > h_g) regimes.

Regarding the inner planet evolution, if we follow IDA20, while e < h_g, the embryo first moves away from the star and then moves slowly inward until e = h_g, from where we can separate the direction of migration into two: first the planet moves outward until Γ_total = 0 when the planet starts moving inward because of the change in the torque sign. On the other hand, when we follow CN08, while e < h_g, the planet moves inward until it reaches e = h_g, at which time it reproduces the same migration direction as for IDA20. When e is non-negligible, Γ_total does not give the direction of evolution of the semimayor axis when the planet follows CN08, but it does when the planet follows IDA20 (see the torques in Fig. 1). This means that the CN08 prescription involves tm instead of t_a in the acceleration expressions of the planets, which can differ in sign as t_m is related to both t_a and t_e (see Eq. (37)). On the other hand, when the embryo orbit starts to circularize, an inward migration is represented by Γ_total < 0 and an outward migration by Γ_total > 0, in both prescriptions, thus t_m and t_a preserve their sign.

By following either of the two prescriptions, the outer planets moves inward. As in the case of the inner planet, the sign of Γ_total does not agree with the direction in which a evolves for a non-circular orbit when we use the CN08 prescription. On the other hand, the planet migration under the IDA20 treatment coincides at any position with a Γ_total < 0, in agreement with an inward migration. For the outer planets, the semimajor axis decreased considerably up to one order of magnitude during the integration time. The sinks in the evolution of a from both prescriptions when the planet reached e = h_g came from the variation in t_a that the planet experienced while migrating inward. Figure 3 shows the initial values of t_a from CN08 and IDA20, which were lower than the values they took when the orbit became circular. In the case of a circular orbit, t_a presents two peaks, one related to a distance that is close to the inner edge of the disk (a ~ 0.0165 au) where Γ_total changes its sign, and the other related to the distance a = r_trans ~ 0.8 au that separates the viscous and irradiated zones of the disk where Γ_total starts to decrease as shown in Fig. 1.

For the inner and outer planet simulations, t_a shows the larger difference within the two prescriptions in the supersonic regime. Oscillations in their absolute values are visible that are associated with the calculation of Γ_total, which differs from one prescription to the other. On the other hand, when the orbit is quasi-circular, both timescales are similar.

When the orbit is quasi-circular, the timescales te and ti are equivalent in the two prescriptions. When the orbit is eccentric while t_i remains constant for CN08, it follows the same increment as t_e for IDA20, because one depends on the other, as shown in Eqs. (27) and (28).

Finally, we explored the effect of orbital inclination through a set of new simulations for an orbit initially inclined by an angle i ≤ h_g. We found no modification in the evolution for a, e, i, t_i, and t_a while the timescale t_e differed slightly within the two prescriptions, although these differences did not produce any considerable change in the orbital evolution or migration directions.

Our analysis shows important differences in the orbital evolution of a planet for CN08 or IDA20. We therefore decided to develop two sets of simulations for a detailed study of the impact of the two prescriptions on the formation and evolution of a planetary system around an object close to the substellar mass limit.

Fig. 2 First 1 Myr evolution of a and e, and the absolute values of ta, te and ti for an inner planet (left panels) initially located at a = 0.016 au with e = 0.1 and an outer planet (right panels) initially located at a = 0.1 au and e = 0.5. The planet masses are MP = MMars. Solid lines indicate the results following the prescriptions from IDA20, and dotted lines show the corresponding results from CN08. The vertical lines indicate the moment at which e = hg, which separates the subsonic (e < hg) and the supersonic (e > hg) regimes.

Fig. 3 Timescale ta as a function of a for a planet with a mass MP = Mars. The black line indicates the case of a circular and copla-nar orbit. The red and blue lines indicate a planet initially located at a = 0.1 au with e = 0.5 following the IDA20 and CN08 prescriptions, respectively. All cases consider disk parameters set to 1 Myr.

4 Simulations of planetary system formation

We present the scenario of planetary system formation starting with a sample of embryos orbiting an evolving central object with a mass of 0.08 M_⨀. We developed a set of 23 N-body simulations following the prescription of IDA20 and another 23 simulations following CN08 with our modified version of the MERCURY code. As external forces, we incorporated the acceleration corrections generated by the interaction within the gas-disk and the embryos of the two prescriptions, as well as those produced by tidal and general relativistic effects which are relevant during rocky planet formation at the substellar mass limit, as we showed in Sánchez et al. (2020). We also included the central object contraction and rotational period evolution as well as a fixed pseudo-synchronization period for embryos.

4.1 Tidal and general relativistic effects

We incorporated tidal effects following the equilibrium tide model from Hut (1981) and Eggleton et al. (1998). We included tidal distortions and dissipation terms, considering the tide raised by the central object on each embryo and by each embryo on the central object and neglected the tide between embryos as follows, $$f_{\omega }=-3\frac{\mu }{r^8}\left[k_{2,*}\left(\frac{M_p}{M_{*}}\right)R_{*}^5+k_{2,\text{p}}\left(\frac{M_{*}}{M_{\text{p}}}\right)R_{\text{p}}^5\right]r,$$(38) $$\begin{array}{cc}{f}_{\text{ae}}& =-3\frac{\mu }{r^{10}}\left[\frac{M_{\text{p}}}{M_{*}}{k}_{2,*}\Delta {\text{t}}_{*}{R}_{*}^5\left(2r\left(r\cdot \upsilon \right)+r^2\left(r\times {\Omega }_{*}+\upsilon \right)\right)\right]\\ & -3\frac{\mu }{r^{10}}\left[\frac{M_{*}}{M_{\text{p}}}{k}_{2,\text{p}}\Delta {\text{t}}_{\text{p}}{R}_{\text{p}}^5\left(2r\left(r\cdot \upsilon \right)+r^2\left(r\times {\Omega }_{\text{p}}+\upsilon \right)\right)\right],\end{array}$$(39)

where k_2,* = 0.307 and k_2,p = 0.305 are the potential Love numbers of degree 2 of the star and the embryos, respectively. For the star we assumed the Love number of an object with a mass at the substellar mass limit, and for the embryos, we assumed the Love number estimated for the Earth (Bolmont et al. 2015). The variable r is the position vector of the embryo with respect to the central object, µ = G(M_⋆ + M_p), G is the gravitational constant, and M_⋆, R_⋆, M_p and R_p are the masses and radius of the star and the protoplanetary embryo, respectively, under the approximation that these objects can instantaneously adjust their equilibrium shapes to the tidal force and considering only up to second-order harmonic distortions (Darwin 1908). The variable υ is the velocity vector of the embryo with respect to the central star, ∆t_⋆ and ∆t_p are the time-lag model constants for the star and the protoplanetary embryo, respectively. The factors k_2,⋆∆t_⋆ and k_2‚p∆t_p, are related to the dissipation factors by $$\begin{array}{cc}{k}_{2,\text{p}}\Delta {\text{t}}_{\text{p}}=\frac{3R_{\text{p}}^5{\sigma }_{\text{p}}}{2G}& k_{2,*}\Delta {\text{t}}_{*}=\frac{3R_{*}^5{\sigma }_{*}}{2G}\end{array}$$(40)

with the dissipation factor for each protoplanetary embryo σ_p = 8.577 × 10⁻⁴³ k⁻¹m⁻²s⁻¹, which is the same dissipation factor as estimated for the Earth (Neron de Surgy & Laskar 1997), and the dissipation factor of the central object is σ_⋆ = 2.006 × 10⁻⁵³ k⁻¹ m⁻² s⁻¹ (Hansen 2010).

We included the rotational evolution (Ω_⋆) and contraction of the central object (R_⋆) from Bolmont et al. (2011) and Baraffe et al. (2015), and fixed each embryo at pseudo-synchronization (Ω_p) following Hut (1981). Thus, in a heliocentric reference frame, tidal interactions lead the precession of the argument of periastron ω, as well as the a and e decays.

We also incorporated the acceleration corrections associated with the precession of periastron caused by the central object as derived from the General Relativity Theory (GRT) (Einstein 1916) as follows, $$f_{\text{GR}}=\left[\left(\frac{G{M}_{*}}{r}-{\upsilon }^2\right)r+4\left(\upsilon .r\right)\upsilon \right],$$(41)

with c the speed of light. Equation (41) was proposed by Anderson et al. (1975), who used the parameterized post-Newtonian theories and reported a relative correction associated with two parameters βand γ, which are equal to unity in the GRT case. We refer to Sánchez et al. (2020) for a detailed description of tidal and relativistic corrections and their associated orbital decay timescales.

4.2 Initial distribution of embryos

The initial spatial distribution of embryos extends from r_ice < r < r_final, where r_ice = 0.23 au is the location of the snow line at 1 Myr and r_final = 5 au as in Coleman et al. (2019). The location of the snow line was computed following the parameterization from Ida et al. (2016) given by r_ice ~ max(r_ic,vis, r_ice,irr), with $$r_{\text{ice,vis}}=1.2{\left(\frac{M_{*}}{m_{\odot }}\right)}^{1/3}{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)}^{4/9}{\left(\frac{{\alpha }_{\text{g}}}{{10}^{-3}}\right)}^{-2/9}au,$$(42)

and $$r_{\text{ice,irr}}=0.75{\left(\frac{M_{*}}{m_{\odot }}\right)}^{-1/3}{\left(\frac{L_{*}}{L_{\odot }}\right)}^{2/3}\text{au},$$(43)

where M_g and L_⋆ are evaluated at 1 Myr. The first embryo was located at a = rice while the location of the remaining consecutively embryos was calculated with a_i+1 = a_i + ∆R_Hill, assuming ∆ to be a randomly integer number between 5 and 10 and R_um = a(2M_emb/3 M_⋆)^1/3, with i = 1, 2, etc.

Our simulations are intended to explore the rocky planet formation from an embryo population and do not include earlier formation stages such as pebble accretion or planetesimals. We therefore set an initial sample of already formed embryos with masses M_emb = 0.16 M_⊕ (~1.5 M_Mars) comparable to those used in previous work on planet formation at the substellar mass limit (e.g., Coleman et al. 2019).

The number of embryos was computed from the total mass of solids M_solid as $$M_{\text{solid}}=2\pi {\displaystyle {\int }_{r_{\text{ice}}}^{r_{\text{final}}}r{\displaystyle {\sum }_{\text{solid}}\text{d}r}},$$(44)

where $${\displaystyle {\sum }_{\text{solid}}=\{\begin{array}{ll}{\displaystyle {\sum }_{\text{g,vis}}{z}_0{\eta }_{\text{ice}}}\hfill & \text{if}r<r_{\text{tran}},\hfill \\ {\displaystyle {\sum }_{\text{g,irr}}{z}_0{\eta }_{\text{ice}}}\hfill & \text{if}r>r_{\text{tran}},\hfill \end{array}}$$(45)

where ∑_solid is the solid surface density profile, z₀ = 0.0153 is the primordial solar abundance of heavy elements (Lodders et al. 2009), and η_ice represents the increase in the amount of solids due to the condensation of water at r > r_ice whose values from Lodders (2003) are η_ice = 1 if r < r_ice and η_ice = 2 if r > r_ice. The factor 2 in Σ_solid is related to the water radial distribution, so that all the bodies located outside the snow line have 50% of water in mass. Then, by solving Eq. (44) we obtained M_solid ~ 7 M_⊕. We considered that ~30% of the mass is the dust that causes the opacity in the disk, which is representative of the typical values of dust masses that have been found in disks around low-mass objects at early ages (e.g., Ward-Duong et al. 2018). Then the initial mass in embryos is ~4.5 M_⊕ and the total number of initial embryos is 28. We neglected the mass of solids between r_init and rice because is negligible in comparison with the mass obtained from solving Eq. (44).

We considered that all embryos initially have e < 0.02 and i < 0.5°, and the initial values of the argument of perias-tron ω, longitude of the ascending node Ω, and mean anomaly M, were randomly determined for each embryo from uniform distributions between 0° and 360°.

Fig. 4 Semimajor axis as a function of eccentricity for embryos at different times before the complete dissipation of the gas at 10 Myr. The color scale represents the initial a of each embryo, and the sizes indicate their masses. The top and bottom panels show the results from simulations including the prescriptions from IDA20 and CN08, respectively. The black lines indicates e = hg, and the vertical lines the inner edge of the disk (rinner = 0.015 au).

4.3 Characterization of N-body simulations

To develop our simulations, we chose the hybrid integrator, which uses a second-order symplectic algorithm to treat interactions between objects with separations greater than 3R_Hill and the Bulirsch-Stöer method to solve close encounters. The time step we adopted corresponds to 1/30th of the orbital period of the innermost body of the simulations, which is 0.08 days. We considered an embryo as ejected from the system when it reached a distance r > 100 au, and we considered that the embryo had collided with the central star when it was closer than 0.0045 au, which corresponds to the maximum radius of the star. We fixed this value for the entire simulation in order to avoid any numerical error for small-perihelion orbits. All simulations ran for 100 Myr in order to analyze the dynamic of planetary systems well after the gas has dissipated from the disk at 10 Myr.

5 Results

In this section, we analyze the gas effects on the dynamical evolution of embryos during the gas-disk lifetime regarding the prescriptions from IDA20 and CN08 independently, and also the dynamics of planetary systems after the gas has been dissipated from the disk. We show the final architectures of planetary systems, focusing on the close-in planet population, and compare it with observational results.

5.1 Gas stage

As discussed in Sect. 3, the gas disk exerts torques on the embryos, allowing them to migrate. The direction of the migration depends on each prescription for the disk model used, as well as on the physical and orbital parameters of the embryos. The timescales for orbital decay in both treatments are comparable when the planetary orbit is quasi-circular and coplanar with the midplane of the disk, but they differ from each other when the orbit is eccentric. This discrepancy has an important effect on the dynamic of embryos, as we discuss below.

Figure 4 shows e as a function of a for the embryos that have survived in the disk at 1 Myr, 5 Myr, and 10 Myr during the gas-disk lifetime, distinguished by their initial semimajor axis and mass. In the close-in population (a < 0.1 au), the embryos decreased their e and migrated inward faster under the CN08 prescription than under the IDA20 treatment. As an example of this behavior, Fig. 5 shows the t_a and t_e damping timescales for an embryo with a mass of M_Mars located at a = 0.02 au for the IDA20 and CN08 prescriptions.

The eccentricity damping timescales are shorter in the CN08 regime. Thus the planets involved have lower eccentricities than those ones under the IDA20 treatment (see Fig. 4). This difference in eccentricity in the supersonic regime causes that for a given time, at a particular location of the planet, the semimajor axis damping timescale is higher for an IDA20 embryo than for a CN08 embryo, as shown in Fig. 5. This explains why the inner embryo population migrated faster and left the region 0.015 < a/au < 0.1. Only one embryo survived in ~17% of the CN08 simulations in quasi-circular and coplanar orbits, but more than one embryo survived in this region in all IDA20 simulations, in some cases, in more eccentric orbits (see Fig. 4).

Figure 6 shows the fraction over the initial number of embryos that ended up inside the cavity (a < 0.015 au) during the gas stage for both prescriptions. In all simulations, between ~45% and 70% of the total number of initial embryos entered the cavity under the CN08 treatment, while just between ~5% and 30% of them did so under the IDA20 prescription. A greater number of embryos entered the cavity under the CN08 treatment for the fast migration and the strong total torque at the inner edge of the disk that the embryos experienced, especially during the first 5 Myr of the disk lifetime.

Approximately 10% of the embryos that entered the cavity under the CN08 prescription collided, and then the remaining 90% collided with the central object. Under the IDA20 treatment, ~10% of the embryos that entered the cavity collided in only ~30% of the simulations, while in the remaining ~70% of the simulations, all the embryos that entered the cavity collided directly with the central object during the integration time.

Figure 7 shows the percentage of the initial number of embryos that collided during the disk lifetime (10 Myr) for every simulation of each of the two prescriptions. The collisions were more frequent for embryos in simulations including the IDA20 prescription (~78%) than those in CN08 (~22%). In both cases, about 98% of the collisions occurred among embryos with a < 0.1 au. Under the CN08 prescription, all the embryos that collided inside the disk entered the cavity. On the other hand, under the IDA20 prescription, this occurred for all of the embryos that collided inside the disk, in later stages during the gas-disk lifetime, given the slow migration rate under this treatment, and thus they survived in the planetary systems.

For a detailed analysis of the collision history of embryos during the gas-disk lifetime, Fig. 8 shows the medians of the accumulated collisions of embryos and with the central object, during and after the gas dissipated from the disk, considering all the simulations for each prescription. During the first 3 Myr, no collisions with the central object occurred in the IDA20 simulations, while less than 10% of these collisions occurred in ~17% of the CN08 simulations, even though the median value is equal to zero. Furthermore, a median of ~10% of the embryos collided in the simulations of each treatment. Later, between 3 and 5 Myr, the IDA20 simulations have a median of collisions among embryos of ~32% and still a median of ~0% of collisions with the central object, while the CN08 simulations present a median of ~21% and ~32%, respectively. Finally, between 5 and 10 Myr, the median of collisions between embryos increase up to ~40% in the IDA20 simulations, while the median remained the same for the CN08 simulations. For the collisions with the central object, IDA20 simulations have a median of less than 10%, while CN08 simulations median value increased up to ~40%.

After this analysis, we point out that in a standard disk model, the prescriptions from CN08 and IDA20 lead to differences within the dynamics of the systems that produces different inner surviving embryo populations not only at the age of disk dissipation, but throughout the entire gas stage.

Fig. 5 Semimajor axis migration timescales ta (top panel) and eccentricity damping timescales te (bottom panel) as a function of e/hg. The red and blue lines indicates the prescriptions from IDA20 and CN08, respectively. Both cases consider an inner planet with mass MMars located at a = 0.02 au. The red dot indicates e = 0.2 from the IDA20 prescription, and the blue dot shows e = 0.05 from CN08. The solid black line represents the limit between subsonic and supersonic regime.

Fig. 6 Percentage of the initial number of embryos that entered the cavity (a < 0.015 au) during the gas stage. Each bin indicates a single simulation. Red bars represent simulations using IDA20, and blue bars show simulations using CN08.

Fig. 7 Percentage of the initial number of embryos that collided during the gas stage. Colors are as in Fig. 6.

Fig. 8 Medians of accumulated collisions of embryos (top panel) and of embryos that collided with the central object (bottom panel) of all IDA20 (red lines) and CN08 simulations (blue lines) before and after the dissipation of the gas-disk at 10 Myr (vertical line).

5.2 Resonances

During the gas stage, all simulations from both prescriptions show mean motion resonances in the inner population (a < 0.1 au). The difference between the two treatments lies in the duration and rupture of these resonances.

Two planets whose semimajor axis a_i and a_j satisfy a_i < a_j, are in commensurable orbits if they follow the relation $$a_{\text{i}}\approx {\left(\frac{q}{p}\right)}^{2/3}{a}_{\text{j}},$$(46)

where p and q give the order of the commensurability p: q.

On the other hand, if a pair of planets satisfies this relation and also presents a libration critical angle, then this pair of planets is in mean motion resonance. In all our simulations, the critical angles that librated were $${\theta }_1=p\overline{{\lambda }_{\text{i}}}-q\overline{{\lambda }_{\text{j}}}-\left(p-q\right)\overline{{\omega }_{\text{i}}},$$(47) $${\theta }_2=p\overline{{\lambda }_{\text{i}}}-q\overline{{\lambda }_{\text{j}}}-\left(p-q\right)\overline{{\omega }_{\text{j}}},$$(48)

where $\overline{{\lambda }_{\text{i,j}}}={\Omega }_{\text{i,j}}+{\omega }_{\text{i,j}}+M_{\text{i,j}}$ and $\overline{{\omega }_{\text{i,j}}}={\Omega }_{\text{i,j}}+{\omega }_{\text{i,j}}$, with Ω_i,j the node longitude, ω_i,j the argument of periastron and M_i,j the mean anomaly of each body involved. In Fig. 9, we show one example of a pair of planets under the IDA20 treatment in mean motion resonance, which fulfilled the requirements of a libration angle and being in commensurable orbit in an interval of integration time between t ~ 4 Myr and t ~ 10 Myr. This pair of planets in mean motion resonance was part of a resonant chain. The pair migrated inward at t ~ 4 Myr after the outer planets involved in the resonant chain collided and caused the remaining innermost planets to migrate inward. When the gas dissipated at t ~ 10 Myr, the bodies migrated toward the star at different times, and the mean motion resonance broke apart.

For both treatments, the most common resonances are on the order of 2:1, 3:2, and 4:3. However, resonances of order 3:1, 4:1, 5:1, 5:2, and 5:3 were found with lower frequency. We show in Fig. 10 the number of resonances of different orders that we found taking into account all the simulations made under the two prescriptions IDA20 and CN08. We found resonant chains of resonances with different orders throughout the simulations.

Planets under the CN08 treatment entered in their resonance during the first 1 Myr and usually remained there for 1.5–2 Myr. When the innermost embryo entered the cavity, all the embryos started to move inward, maintaining their respective resonances. The fast migration was mainly caused by the strong gas-disk torques at early stages at distances close to the inner edge of the disk, and second by tidal effects when the embryo was getting closer to the star. The resonances broke apart when the embryos collided with the star. However, we cannot confirm that such a collision caused the breaking as the inward migration occurred in a short period of time (~40 000 yr). Then, at the age of t = 5 Myr, no resonant chain can be found in any of the systems because all the inner bodies have already migrated inward and collided with the central object. One example of this is shown in the left panel of Fig. 11, which shows different pairs of planets in a resonant chain that migrated inward after the innermost body entered the cavity. They migrated inward, maintaining their resonances, and finally collided with the star. In most of the cases, no resonant chain was formed later due to a lack of inner planets (a < 0.1 au). In less than ~10% of the simulations, one pair of resonant planets could be found from 8 Myr to 10 Myr, which also finally broke up after the gas dissipated from the disk.

For planets under the IDA20 treatment, we can distinguish different moments associated with the formation of resonant pairs. In all simulations, some pair of planets entered their resonance during the first 1 Myr. Some of them finally broke apart because one or more of the embryos involved entered the cavity or because they experienced collisions with the embryos involved in the resonant chain, while others remained in resonance until the gas dissipated. On the other hand, some other resonant chains were formed after 4 Myr of integration time, which usually remained during the entire gas stage. One example of this behavior is shown in the right panel of Fig. 11, which shows different pairs of planets in a resonant chain that entered in their resonance at different times. Some resonances broke apart when one pair of embryos collided and others because the innermost body entered the cavity. Thousands of years after the innermost planet started to migrate inward, the resonances broke apart and the planets involved started to migrate inward. Then, right after t ~ 10Myr, when the gas dissipated, all the surviving resonant chains broke up because the innermost body entered the cavity and started to migrate inward, breaking the other surviving resonances. The migration was mainly caused by the total gas-disk torques and second, by tidal effects when the planet approached the star. In this case, the planets experienced a slower migration than in the CN08 example (in which the bodies entered the cavity at t ~ 4 Myr), as at t = 10 Myr the strength of the negative gas-disk torques is lower. Even though no resonant chain was formed after the gas-disk dissipated, some pairs of surviving planets remained close to commensurabilities throughout the integration time.

Fig. 9 Example of one mean motion resonance under the IDA20 treatment. The critical angle libration (top panel) and commensurability 3:2 (bottom panel) of a pair of planets in resonance are shown during the integration time. The vertical lines represent t = 10 Myr.

Fig. 10 Histogram with the total number of resonances of different orders taking into account the total number of simulations made for each prescription, IDA20 or CN08. An offset was added to the IDA20 resonances for clarity.

Fig. 11 Examples of planets in commensurable orbits that represent typical resonant chains from a simulation under the CN08 treatment (left) and the IDA20 treatment (right).

5.3 Post-gas stage

The main difference in the post-gas evolution in both prescriptions occurs at distances in the range 0.015 < a/au < 0.1. There, a large population of embryos remains for the IDA20 prescription with eccentricities up to ~0.08, while only a few embryos survive in the CN08 treatment and present quasi-circular orbits. At larger distances (a > 0.1 au), the effects of the gas are less relevant and embryos retain their original quasi-circular orbits. These results are summarized in Fig. 12, which shows the surviving embryos, distinguished by their initial semimajor axis and mass, in all simulations for the two prescriptions at ages of 10 Myr, 50 Myr, and 100 Myr.

After the dissipation of the gas, no other embryo entered the cavity in simulations under the CN08 prescription, while many of them did in simulations under the IDA20 treatment. All the remaining embryos inside the cavity (a < 0.015 au) eventually collided with the central object. On the one hand, in the simulations under the CN08 treatment, the embryos that entered the cavity just before the gas-disk dissipation were not part of a resonant chain and migrated inward slowly given the weak intensity of the total torque at the inner edge of the disk and the stellar tide. For the simulations under the IDA20 treatment, their evolution toward the central object was mainly caused by gravitational interactions with outer embryos that were part of a resonant chain and were immersed in the gas-disk just before its dissipation. The evolution was also affected by stellar tides when they approached the star (a < 0.008 au).

As explained in Fig. 8 in Sect. 5.1, most of the collisions occured before the gas dissipation. After this, embryos under the CN08 treatment experienced no collisions among themselves, but with the central object, reaching a median value of almost 50% considering all the simulations. On the other hand, embryos under the IDA20 prescription continued to experience collisions among themselves and with the central object up to 50 Myr, increasing the median value up to ~45% and ~15%, respectively. We note that after 50 Myr, no simulation showed a collision between embryos in either of the two prescriptions, while ~22% and ~40% of them presented one collision with the central object under the CN08 and IDA20 treatments, respectively, even though the median values are the same as at 50 Myr.

In Fig. 13, the final architectures of the surviving planets at 100 Myr are shown for IDA20 and CN08. Planets are distinguished by size according to their final masses. Most of the CN08 planets remained at their initial mass, while IDA20 planets present a wide range of masses of up to ~1 M_⊕ exclusively in the inner region with a < 0.1 au. The isolated habitable zone (IHZ) of the systems at 100 Myr and at 1 Gyr are also included. We calculated the IHZ according to Selsis et al. (2007) and Barnes et al. (2013) and adapted it for our central object of 0.08 M_⊕ as explained in Sánchez et al. (2020).

Finally, we compared the spatial distribution of the resulting embryos with the position of the IHZ. Because of the expected evolution of the IHZ toward the central object, we considered its corresponding positions at 100 Myr and 1 Gyr assuming that in this period, the migration and e damping timescales due to tidal effects are long enough to avoid a change in the location of the embryos. The resulting systems from the IDA20 treatment present several embryos inside the IHZ for all the simulations, while from the CN08 prescription only ~10% of the simulations present just one single embryo inside the IHZ. On the one hand, we note that the embryos inside the 100 Myr IHZ will not migrate inward following the evolving IHZ. Thus, after few million years, these planets will be located farther away from the external limit of the habitable zone. On the other hand, the embryos located inside the 1 Gyr IHZ were unprotected for a considerable amount of time until the evolving IHZ reached them. The fact that the planets are exposed to stellar radiation inside the inner edge of the IHZ for a long period of time makes the survival of water among other volatiles on the surface of the planets challenging because they experience greenhouse runaway (e.g., Luger & Barnes 2015). Thus, the potential habitability of these planets should be carefully evaluated.

Fig. 12 Semimajor axis vs. eccentricity for the surviving embryos at ages of 10 Myr, 50 Myr, and 100 Myr for the IDA20 and CN08 prescriptions (top and bottom panels, respectively). The colored scale indicates the initial values of the semimajor axis, and the points sizes represent the embryo mass according to the upper labeled black dots. The vertical line indicates the position of the inner radii of the gas-disk.

Fig. 13 Final planetary architectures for each simulation from the IDA20 prescription (leftpanel) and CN08 prescription (rightpanel). The color and size of the points represent the final masses of the planets as shown in the upper labeled panels. The cream band represents the IHZ at 100 Myr, and the pink band shows the IHZ at 1 Gyr.

Fig. 14 Comparison of cumulative distributions of survival planets in the IDA20 (red line) and CN08 (blue line) simulations at 10 Myr (top panels) and at 100 Myr (bottom panels), and the confirmed Earth-like planets around stars with 0.08 < M⋆/M⊙ < 0.14 M⊙ (solid black line) and 0.14 < M⋆/M⊙ < 0.5 M⊙ (dotted black line). All planets are located at a < 0.1 au. The color shadows correspond to the Poissonian errors in IDA20 (red), CN08 (blue) and confirmed exoplanets around M⋆ < 0.14 M⊙ (orange) and M⋆ < 0.5 M⊙ (gray).

5.4 Comparison with confirmed exoplanets

In this section, we compare the period distributions of the innermost surviving planets (a < 0.1 au) with the corresponding confirmed Earth-like exoplanets from the Exoplanet Archive¹. We assumed a mass/radius cutoff of 2 M_⊕/2 R_⊕ and considered the sample of exoplanets detected through transit and radial velocity techniques.

As in previous works that studied gas-disk interactions down to the substellar limit (e.g., Coleman et al. 2019; Miguel et al. 2020) we compared our results at the end of the gas stage. The top panels of Fig. 14 show the cumulative distribution of the period ratio of adjacent planets and the distribution of periods of the innermost planet from each simulation under IDA20 and CN08 at 10 Myr together with the corresponding distributions of confirmed terrestrial exoplanets orbiting stars with masses 0.08 < M_⋆/M_⊙ < 0.14 and 0.14 < M_⋆/M_⊙ < 0.5. We included the Poissonian errors for the numerical simulations and for the exoplanet population around the less and the more massive stars. These errors were calculated with ± the square of the cumulative planets discretized by period and period ratio, respectively. The close-in planetary population obtained following the CN08 prescription around a star of 0.08 M_⊙ agrees better in the distribution of the period ratio of adjacent planets with the corresponding population from confirmed exoplanets around stars with masses M_⋆ < 0.5 M_⊙. Despite this similarity, we caution that this comparison was performed against observations of field stars that are older than 1 Gyr and also span a wider mass range than the central star of our simulations. On the other hand, IDA20 and CN08 both overestimate the number of planets with innermost periods in comparison with the observations.

We aim to show how the comparison with the same observations changes at the end of our simulations. The bottom panels of Fig. 14 show the distribution of the period ratio of adjacent planets and the distribution of periods of the innermost planet from each IDA20 simulation at 100 Myr, together with the corresponding distributions of confirmed Earth-like exoplanets. In this case, the close-in planetary population obtained following the IDA20 prescription shows a distribution of the period ratio of adjacent planets that agrees with the distribution of the confirmed exoplanets around stars with masses M_⋆ < 0.14 M_⊙. A K-S (Kolmogórov-Smirnov) test was performed by applying 500 bootstrap realizations. We found that in only ~15% of the realizations we can reject that both distributions came from the same distribution with a 99% probability. A corresponding analysis is not possible for the simulations following the CN08 prescriptions because, as we showed in Fig. 13 in Sect. 5.3, not more than one planet at a < 0.1 au survived at the end of all these simulations. We note that considering gravitational, tidal, and general relativistic effects, the planetary configurations are not expected to change significantly after 100 Myr. Thus, it is valid to compare the simulations with the known exoplanet population around old stars.

We were unable to find a good agreement within the innermost planet period distributions regarding the surviving planets in our simulations under either of the prescriptions (CN08 and IDA20) and the confirmed exoplanets around either of the selected samples. On the one hand, ~80% of the surviving innermos planets in simulations under the CN08 treatment have periods longer than 10 days, while all the surviving innermost planets in simulations under the IDA20 prescription have periods between 4 days and 10 days. On the other hand, the innermost observed exoplanet samples have periods between 0.4 days and 10 days. Thus, the innermost planet period distribution taken from the surviving planets in simulations under IDA20 remains closer to the distribution taken from either of the two exoplanet samples, even though it presents a deficit of inner planets.

The deficit of inner planets (a < 0.025 au) in simulations under the IDA20 treatment would not change on a timescale of a billion years. Considering tidal effects, the resulting planets with a > 0.025 au would decrease their semimajor axis in more than ~ 10 Gyr. One possible explanation for this deficit might lie in the standard disk model. Changing the disk lifetime, the turbulent and accretion α, as well as including an evolving inner edge of the gas-disk might allow the survival of more planets located closer to the star. However, we expect that this consideration will not produce qualitative differences in the agreement found within the period ratio of adjacent planets for IDA20 planets and the confirmed Earth-like planets because with the incorporation of inner planets in the systems, the systems would remain in compact configurations.

We conclude that it is not correct to compare the resulting planetary configurations at the end of the gas stage with the confirmed exoplanet population. Moreover, this comparison should be made considering an analogous mass for the central object of the systems. We therefore ran the simulations for 100 Myr to be able to compare the planetary architectures with the confirmed exoplanet population. This is as a representative time of the final stage of the system because the number of surviving planets remained almost constant during the last 50 Myr, as shown in Fig. 8. Furthermore, regarding gravitational interactions and tidal effects, the planetary configuration would not change significantly in 1 Gyr given the low eccentricity values, the mass range, and the semimajor axis of the resulting planets.

6 Conclusions and discussions

We treated rocky planet formation around an evolving 0.08 M_⊙ star by using 100 Myr long N-body simulations that incorporate tides and general relativistic effects, as well as gas-disk interactions during the first 10 Myr. The disk was simulated according a standard disk model, and its effect over the embryos was simulated by independently following the torque prescriptions from IDA20 and CN08.

We found a resulting close-in (a < 0.1 au) planet population under IDA20 that did not survive under the CN08 treatment at the end of our simulations (100 Myr). Moreover, the IDA20 close-in planet population agrees with the period ratio of adjacent confirmed Earth-like exoplanets around stars with M_⋆ < 0.14 M_⊙. Furthermore, we note that to compare the simulated planetary architectures with the confirmed exoplanet population, it is important to properly constrain the comparison within samples in similar stellar mass ranges and using simulations whose integration times are long enough to be compared accordingly with the old observed planetary systems.

We also found a surviving planet population located in the habitable zone of IDA20 systems that is not present in simulations under the CN08 treatment. However, this population should be studied carefully, taking the dynamical evolution history of the planets and the location of the evolving IHZ into account. Considering the influence of gravitational, tidal, and general relativistic effects, no significant changes over time are expected in the planetary configurations after 100 Myr. Thus, our resulting planets are therefore exposed to stellar radiation due to their location inside the inner edge of the IHZ for 1 Gyr. Then, the planets located in the IHZ may not satisfy the conditions to be considered potentially habitable (e.g., Luger & Barnes 2015).

We obtained several mean motion resonant chains that survived during the entire gas stage only under the IDA20 prescription. After the gas dissipated, all these resonant chains broke apart. However, some pairs of planets remained close to com-mensurabilities. This allows the existence of compact planetary systems, which agrees with the period ratio of an adjacent close-in confirmed Earth-like population around less massive stars (M_⋆ < 0.14 M_⊙).

Even though the resonant breaking mechanism is still not well understood, we aim to discretize three different breaking scenarios that we found in our work:

• Collision with the central object. The resonant chains broke apart when the embryos collided with the central object. However, given the fast inward migration, we cannot distinguish whether the breaking was due to the collision with the central object or due to the planets entering the cavity. This is the case for all resonant breaking in simulations under the CN08 treatment;

• Collision among the outer embryos. The resonances broke apart because the outer embryos involved in the mean motion resonance chain experienced collisions with other embryos during the gas-disk lifetime. This was the case for some of the resonant breaking in simulations under the IDA20 treatment;

• Disk dispersal/disk absence. Either the resonances broke apart when the innermost planet entered the cavity while the outer planets were still immersed in the disk, or after disk dispersal, when the innermost planet involved in a resonance was located inside the cavity or close to the inner edge of the disk and the outer planet was still immersed in the disk just before the gas-disk dissipation. This was the case of the remaining breaking cases in simulations under the IDA20 treatment.

We highlight that the mean motion resonance breaking after the gas-disk dissipation was also found by other authors who studied planet formation around Sun-like stars, even when stellar tides were not considered (e.g., Izidoro et al. 2017, 2021).

We analyzed tidal and general relativistic effects as well as gas-disk interactions. During the gas stage, the gas-disk interactions play a primary role in the dynamical history of the planets. However, tidal and general relativistic effects give a more detailed model of the orbital dynamic of the planets both during the gas stage and after the gas disk dissipated (Sánchez et al. 2020).

In our model, we used a set of initial parameters for the standard gas-disk model. We aim to explore different scenarios for disk masses and lifetimes as well as the impact of different values for the accretion and turbulent α, as in Matsumura et al. (2021), onto the resulting planetary configurations in future works. We expect that for different initial conditions, we may reproduce the innermost exoplanet period distributions if an innermost planet population survives in our new numerical simulations.

This study allow us to better understand the rocky planet formation at the substellar mass limit. We conclude that within the framework of parameters we explored, it would be more appropriate to use the IDA20 prescriptions to treat the gas-disk interactions to study rocky planet formation at the substellar mass limit.

Acknowledgements

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. MS and GdE acknowledges the partial financial support by Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT) by PICT 201-0505, Argentina, and the partial financial support by Facultad de Ciencias Astronómicas y Geofísicas de la Universidad Nacional de La Plata, Argentina (FCAGLP-UNLP) and Instituto de Astrofísica de La Plata (IALP) for extensive use of their computing facilities. JJD acknowledges funding from the MIA program at Universidad de la República, Uruguay, and is grateful for the hospitality and support of FCAGLP-UNLP where part of this research was carried out. We thank also to the referee Brasser R. for his constructive comments.

References

All Figures

	Fig. 1 Maps of the total normalized torque Gtotal/G0 exerted by the gas onto embryos in a 1 Myr old disk for: coplanar and circular orbits (top panels), orbits with e = 0.1 following the prescriptions from IDA20 (middle panels), and orbits with e = 0.1 following the prescriptions from CN08 (bottom panels). The right panels show zoom-ins of the inner part of the disk.
In the text

Fig. 2 First 1 Myr evolution of a and e, and the absolute values of ta, te and ti for an inner planet (left panels) initially located at a = 0.016 au with e = 0.1 and an outer planet (right panels) initially located at a = 0.1 au and e = 0.5. The planet masses are MP = MMars. Solid lines indicate the results following the prescriptions from IDA20, and dotted lines show the corresponding results from CN08. The vertical lines indicate the moment at which e = hg, which separates the subsonic (e < hg) and the supersonic (e > hg) regimes.

In the text

	Fig. 3 Timescale ta as a function of a for a planet with a mass MP = Mars. The black line indicates the case of a circular and copla-nar orbit. The red and blue lines indicate a planet initially located at a = 0.1 au with e = 0.5 following the IDA20 and CN08 prescriptions, respectively. All cases consider disk parameters set to 1 Myr.
In the text

Fig. 4 Semimajor axis as a function of eccentricity for embryos at different times before the complete dissipation of the gas at 10 Myr. The color scale represents the initial a of each embryo, and the sizes indicate their masses. The top and bottom panels show the results from simulations including the prescriptions from IDA20 and CN08, respectively. The black lines indicates e = hg, and the vertical lines the inner edge of the disk (rinner = 0.015 au).

In the text

Fig. 5 Semimajor axis migration timescales ta (top panel) and eccentricity damping timescales te (bottom panel) as a function of e/hg. The red and blue lines indicates the prescriptions from IDA20 and CN08, respectively. Both cases consider an inner planet with mass MMars located at a = 0.02 au. The red dot indicates e = 0.2 from the IDA20 prescription, and the blue dot shows e = 0.05 from CN08. The solid black line represents the limit between subsonic and supersonic regime.

In the text

	Fig. 6 Percentage of the initial number of embryos that entered the cavity (a < 0.015 au) during the gas stage. Each bin indicates a single simulation. Red bars represent simulations using IDA20, and blue bars show simulations using CN08.
In the text

	Fig. 7 Percentage of the initial number of embryos that collided during the gas stage. Colors are as in Fig. 6.
In the text

	Fig. 8 Medians of accumulated collisions of embryos (top panel) and of embryos that collided with the central object (bottom panel) of all IDA20 (red lines) and CN08 simulations (blue lines) before and after the dissipation of the gas-disk at 10 Myr (vertical line).
In the text

	Fig. 9 Example of one mean motion resonance under the IDA20 treatment. The critical angle libration (top panel) and commensurability 3:2 (bottom panel) of a pair of planets in resonance are shown during the integration time. The vertical lines represent t = 10 Myr.
In the text

	Fig. 10 Histogram with the total number of resonances of different orders taking into account the total number of simulations made for each prescription, IDA20 or CN08. An offset was added to the IDA20 resonances for clarity.
In the text

	Fig. 11 Examples of planets in commensurable orbits that represent typical resonant chains from a simulation under the CN08 treatment (left) and the IDA20 treatment (right).
In the text

Fig. 12 Semimajor axis vs. eccentricity for the surviving embryos at ages of 10 Myr, 50 Myr, and 100 Myr for the IDA20 and CN08 prescriptions (top and bottom panels, respectively). The colored scale indicates the initial values of the semimajor axis, and the points sizes represent the embryo mass according to the upper labeled black dots. The vertical line indicates the position of the inner radii of the gas-disk.

In the text

	Fig. 13 Final planetary architectures for each simulation from the IDA20 prescription (leftpanel) and CN08 prescription (rightpanel). The color and size of the points represent the final masses of the planets as shown in the upper labeled panels. The cream band represents the IHZ at 100 Myr, and the pink band shows the IHZ at 1 Gyr.
In the text

Fig. 14 Comparison of cumulative distributions of survival planets in the IDA20 (red line) and CN08 (blue line) simulations at 10 Myr (top panels) and at 100 Myr (bottom panels), and the confirmed Earth-like planets around stars with 0.08 < M⋆/M⊙ < 0.14 M⊙ (solid black line) and 0.14 < M⋆/M⊙ < 0.5 M⊙ (dotted black line). All planets are located at a < 0.1 au. The color shadows correspond to the Poissonian errors in IDA20 (red), CN08 (blue) and confirmed exoplanets around M⋆ < 0.14 M⊙ (orange) and M⋆ < 0.5 M⊙ (gray).

In the text