© The Authors 2025

1 Introduction

Meteorites are divided into two broad categories: chondrites and achondrites (with iron meteorites included in the latter). The parent bodies of achondrites formed in the first megayear (Myr) of the Solar System’s existence. They contained a large abundance of short-lived radioactive elements, whose decay liberated enough energy to melt these bodies (Neumann et al. 2012). Under the effect of gravity, these bodies then differentiated into an iron core, a silicate mantle and a crust, akin to the structure of asteroid 4 Vesta. The parent bodies of chondrites formed later, so that most short-lived radioactive atoms had already decayed. Consequently, the energy liberated in their interior was not as high; temperatures increased significantly, up to several hundred degrees, but not enough to cause melting. Thus, these bodies preserved their original petrological structure, made of chondrules and matrix and refractory inclusions. From measurements of the ages of formation of the individual chondrules, the chondrite parent bodies formed approximately 2–4 My after the beginning of the Solar System (Fukuda et al. 2022; Piralla et al. 2023). The formation of planetesimals so late in the history of the protoplanetary disc is problematic, because dust particles (i.e. the individual chondrules, or the matrix grain aggregates) tend to drift radially towards to Sun on a timescale that is much shorter than a million years. We might consider that the dust lost by radial drift at a given location was substituted by dust drifting from larger radial distances in the disc; however, this simple view is negated by another crucial observation in meteorite science. The analysis of nucleosynthetic isotopic anomalies inherited from the interstellar medium reveals that both chondrites and achondrites are distinguished into two distinct groups, dubbed CC – as the measured anomalies are like those of carbonaceous chondrites – and NC – as the anomalies are like those of non-carbonaceous chondrites (Warren 2011; Kruijer et al. 2017). The existence of an isotopic dichotomy such as that described above implies that the objects of the NC and CC groups formed at different locations of the disc, incorporating distinct materials. Because CCs exhibit evidence of prominent water alteration, whereas the isotopic properties of NC meteorites are similar to those of Earth and Mars, the common interpretation is that the NC bodies (both achondrites and chondrites) formed in the inner Solar System and the CC bodies (again both achondrites and chondrites) formed in the outer Solar System, where water could be accreted in solid form (Kruijer et al. 2017; Nanne et al. 2019). An implication of these considerations is that the material that formed the NC chondrites in the inner Solar System at a later time could not have drifted from the outer disc; otherwise it would have exhibited CC isotopic properties. The formation of Jupiter, opening a gap in the disc (Kruijer et al. 2017), or the formation of a pressure bump near the snowline (Brasser & Mojzsis 2020) must have blocked the drift of particles from the outer disc to the inner disc. Again, meteorites provide evidence that only particles smaller than 200 microns could pass through such dynamical barriers, with negligible effect on the overall mass balance (Haugbølle et al. 2019).

Thus, in the absence of a refilling mechanism, the NC dust in the inner part of the protoplanetary disc had to survive for millions of years without drifting into the Sun. This would require the existence of a pressure bump in the disc where dust could be trapped. A pressure bump may be generated near 1 au if the disc evolves primarily under the effects of magnetised winds and if the mass removal rate in these winds increases with decreasing heliocentric distance (Suzuki et al. 2010; Ogihara et al. 2018). The existence of a unique pressure bump, however, is not completely satisfactory (Kleine et al. 2020; Schneider et al. 2020). Four kinds of NC chondrites exist, with somewhat different isotopic anomalies: the prominent ordinary and enstatite chondrite groups, followed by the Rumuruti (R) and Kakangari (K) groups, which are represented to a lesser extent in the meteorite collections. Together with the fact that Enstatite chondrites have essentially identical isotopic properties of a class of NC achondrites called Aubrites, this implies that the dust did not get completely mixed up, even on a timescale of millions of years. This strongly suggests that it was not one, but several pressure bumps, which had to form and persist in the inner part of the protosolar disc.

Multiple rings of dust are observed in protoplanetary discs (Andrews et al. 2018), suggestive of a multitude of pressure bumps. A long-debated question is whether the formation of these pressure bumps is due to the formation of massive planets opening gaps in the disc (Bae et al. 2017) or to spontaneous non-ideal magnetohydrodynamic (MHD) effects (Béthune et al. 2017; Riols & Lesur 2018; Cui & Bai 2021). We remark that spontaneous generation of rings has been obtained (by Béthune et al. 2017; Riols & Lesur 2018; Cui & Bai 2021) in cases where ampipolar diffusion is the dominant non-ideal effect. This result suggests that such rings would not not form in the inner disc because the MHD dynamics there is thought to be dominated by Hall and Ohmic dissipation, rather than ambipolar diffusion (Armitage 2011).

The rings are observed only in the outer part of protoplanetary discs and, due to resolution effects, it is not possible to assess observationally whether they are common also in the inner part. Also, none of the terrestrial planets was massive enough to open a gap in the disc, even if already present with its current mass, because their Hill sphere is much smaller than the vertical scale height of the disc. Therefore, Jupiter remains the main suspect for the potential formation of rings in the inner disc. Jupiter obviously forms a prominent pressure bump outside of the main gap that opens along its orbit, but this is of no help with respect to the preservation of dust in the inner system. It has been shown, however, that massive planets can also open secondary gaps, separated by density bumps, also inwards of the main gap (Bae et al. 2016, 2017; Miranda & Rafikov 2020a).

Most studies of the process of secondary gap formation have been conducted in 2D equatorial simulations (Bae et al. 2017; Zhang et al. 2018; Miranda & Rafikov 2019, 2020a; Ziampras et al. 2023) and show that three necessary conditions have to be met: (i) the disc must have a very small effective turbulent viscosity; (ii) the disc cooling time normalised by the local orbital frequency, denoted β_c, has to be either small (β_c < 0.01) or large (β_c > 10) (Miranda & Rafikov 2020a; Ziampras et al. 2023); and (iii) the planet must be massive enough (with the mass depending on the disc’s viscosity, β_c, and on the disc’s scale height.

Including a prescription of radiative effects gives the appropriate cooling time associated to disc properties and is a necessary step in capturing the dynamics of the formation of secondary gaps in specific discs. Taking into account the cooling from the disc surface (Ziampras et al. 2020) and adding radiation transport along the disc midplane (Miranda & Rafikov 2020a; Ziampras et al. 2023) have improved the treatment of disc thermodynamics in 2D simulations. However, a self-consistent treatment of radiative effects in 3D simulations of secondary gap opening has not yet been considered. To our knowledge, the only 3D simulations with an appropriate β_c from disc surface irradiation and a Jupiter-mass planet in a massive, inviscid protoplanetary disc were presented in Bae et al. (2016). This work indicated several waves in the disc’s surface density distribution inwards of Jupiter’s orbit. The authors, however, did not perform an analysis of the pressure profile in the disc to verify whether dust-trapping pressure bumps formed.

Given the relevance of this problem to understand the formation of different classes of NC chondrites, we present an exploration of the capability of Jupiter to form multiple dust-trapping rings in the inner Solar System using a self-consistent treatment of radiative effects in a disc with low bulk viscosity. Moreover, since low-viscosity discs transport gas towards the star at a very low rate with respect to typical accretion rates observed in young stars (Hartmann et al. 1998; Manara et al. 2016), we decided to extend our model by including non-ideal MHD effects. In this case, accretion flow onto the star is obtained by angular momentum removal by disc winds (Béthune et al. 2017; Lesur 2021b) and our aim is to explore the feedback of the magnetic field on secondary gap formation.

The paper is organised as follows. In Section 2, we present our simulation tools and methods. In Section 3, we reproduce the work of Bae et al. (2016) and compute the pressure profile on the midplane. In Section 4, we present a different simulation, where we adopt a density comparable to the minimum mass solar nebula at the planet location, instead of the massive disc used by Bae et al. (2016). In addition, we compute β_c self-consistently using flux-limited energy diffusion. In Section 5, we extend our model including non-ideal MHD effects. Finally, in Section 6, we progressively reduce the mass of the planet to see at which mass a planetary core might have started to form dust-trapping rings in the inner Solar System.

2 Models

2.1 Physical model

The protoplanetary disc (PPD hereafter) is treated as a non-self-gravitating gas whose motion is described by the Navier-Stokes equations. We used two grid-based codes: the fargOCA code (Lega et al. 2014), which integrates a full energy equation taking into account viscous heating and radiative cooling and the FARGO3D code (Benítez-Llambay & Masset 2016), which includes non-ideal MHD effects. The two codes are very similar in the sense that they use the same hydrodynamical solvers based on the operator splitting procedure (Stone et al. 1992) and they both benefit from the FARGO algorithm (Masset 2000) for a computation of the time step suited to gas rotating with quasi-Keplerian motion. The reason for the use of the two codes solely depends on the different physical modules that they implement.

In both cases, we used spherical coordinates (r, θ, φ), where r is the radial distance from the star (which is at the origin of the coordinate system), θ the polar angle measured from the z axis (the colatitude), and φ is the azimuthal coordinate measured from the x axis. The midplane of the disc is located at the equator, $\theta =\frac{\pi }{2}$. We are working in a coordinate system that rotates with the angular velocity of a planet of mass, m_p, $${\mathrm{\Omega }}_k=\sqrt{\frac{G(M_{\star }+m_p)}{r_p{}^3}},$$

where M_⋆ is the mass of the central star, G is the gravitational constant, and r_p is the semi-major axis of a planet assumed to be on a circular orbit. We consider a planet orbiting a Solar mass star and, therefore, in the following, we refer to the mass of the star as M_⊙.

The dynamics of the gas is provided by the integration of the Navier-Stokes equations composed by a continuity equation and a set of three equations for the momenta (Lega et al. 2014; Benítez-Llambay & Masset 2016). Additionally, we consider the gas internal energy density as e = ρc_vT, where ρ and T are the gas volume density and the gas temperature, while c_v is the specific heat at constant volume. The equation of state is that of an ideal gas of pressure of P = (γ − 1)e with adiabatic index of γ = 1.4. In the following, we use either of two prescriptions for the integration of the energy equation:

1. Changes in the internal energy density e are due to adiabatic (de)compression, viscous heating, (Mihalas & Mihalas 1984) and a relaxation to the initial temperature, T₀, on a cooling timescale of ${\tau }_{\mathrm{c}}={\beta }_c{\mathrm{\Omega }}_k^{-1}$, $${\mathrm{\partial }}_{t}e+\mathrm{\nabla }\cdot (e\mathbf{v})=-P\mathrm{\nabla }\cdot \mathbf{v}+Q^{+}-\rho c_v\frac{(T-T_0)}{{\tau }_{\mathrm{c}}},$$(1)

   where v is the gas velocity, while the terms −P∇ · v and $Q^{+}\equiv \overline{\overline{\tau }}\mathrm{\nabla }v$ are the compressional and viscous heating ($\overline{\overline{\tau }}$ is the viscous stress tensor defined in Mihalas & Mihalas 1984), respectively. The dimensionless cooling time β_c is either a constant parameter or a function depending on disc properties, as described below.

2. We consider the evolution of the internal energy and of the thermal radiation energy, E_r (code fargOCA), in the flux-limited diffusion approximation (FLD, Levermore & Pomraning 1981) using the so-called two temperature approach (Commerçon et al. 2011), expressed as $$\{\begin{array}{l}{\mathrm{\partial }}_t{E}_{\mathrm{r}}-\mathrm{\nabla }\cdot \mathit{F}=\rho {\kappa }_{p}c(a_r{T}^4-E_{\mathrm{r}}),\\ {\mathrm{\partial }}_{t}e+\mathrm{\nabla }\cdot (e\mathbf{v})=-P\mathrm{\nabla }\cdot \mathit{v}-\rho {\kappa }_{p}c(a_r{T}^4-E_{\mathrm{r}})+Q^{+},\end{array}$$(2)

   where $\mathit{F}=\frac{c\lambda }{\rho {\kappa }_r}\mathrm{\nabla }{E}_{\mathrm{r}}$ is the radiation flux vector and λ the flux limiter (Kley et al. 2009). We use κ_p and κ_r to indicate the Planck and the Rosseland mean opacity, respectively, with a_r as the radiation constant and with c as the speed of light. In this paper, we consider κ = κ_p ≡ κ_r (see Bitsch et al. 2013) and use the opacity law of Bell & Lin (1994), where a dust to gas ratio of 0.01 is assumed.

   We did not include the irradiation heating from the central star, so that we could reduce the computational cost with no significant impact on the gas dynamics in the planet vicinity (see Lega et al. 2015). However, we did consider low-viscosity disc with stellar irradiation as the main source of heating. Therefore, in our 3D simulations, we mimicked the thermal equilibrium structure obtained in discs heated by the star, by suitably fixing the disc surface temperature (see Section 4 and Lega et al. 2024). When considering the energy evolving according to Eq. (2), the disc cooling time, β_c, fully depends on the disc properties and can be estimated a posteriori as explained below. In the following we refer to simulations as “fully radiative” with the above prescription for the energy computation.

2.2 Cooling time approximation

When computing an energy equation with temperature relaxing towards the initial values (Eq. (1)), it is customary to consider the cooling time, β_c, as a constant parameter (see e.g. Miranda & Rafikov 2020b; Ziampras et al. 2023). However, when studying specific discs (or aiming to make comparisons to the observed disc structures), it is important to compute the cooling time from disc properties (Zhang et al. 2018; Miranda & Rafikov 2020a; Ziampras et al. 2023).

2.2.1 Cooling through disc surfaces

The cooling process occurs through blackbody emission of dust grains. In vertically integrated discs the radiative cooling power per unit surface is given by Menou & Goodman (2004), expressed as $$Q_{cool}=-\frac{\sigma T^4}{{\tau }_{eff}},$$

where σ = a_rc/4 and τ_eff is the effective optical depth (Hubeny 1990), expressed as $${\tau }_{eff}=\frac{3}{8}\tau +\frac{\sqrt{3}}{4}+\frac{1}{4\tau },$$

to take into account the optically thick and thin limits, while τ = 1/2κΣ is the optical depth at disc midplane.

The discs cools through surfaces and, following Ziampras et al. (2020) the surface cooling time can be written as $${\beta }_c^{surf}={\mathrm{\Omega }}_K\frac{\mathrm{\Sigma }{c}_{v}T}{|2Q_{cool|}}\equiv {\mathrm{\Omega }}_K\frac{\mathrm{\Sigma }{c}_v{\tau }_{eff}}{2\sigma T^3}.$$(3)

The formula has been extended to a three dimensional (3D) disc by Lyra et al. (2016) by considering the radiative timescale $t_r=\frac{e}{\mathrm{\nabla }\cdot \mathit{F}}$ (with e and F defined in Eq. (2)) and integrating over a spherical volume of radius, H¹, $${\beta }_c^H={\mathrm{\Omega }}_K\frac{\rho c_{v}H{\tau }_{eff}}{3\sigma T^3}.$$(4)

In this case, the optical depth τ_eff is computed for grid cells at height z above (below) the midplane as $\tau (z)={\int }_z^{z_{max}}\rho (z^{\mathrm{\prime }})\kappa (z^{\mathrm{\prime }})\mathrm{d}{z}^{\mathrm{\prime }}(\tau (z)={\int }_{z_{min}}^z\rho (z^{\mathrm{\prime }})\kappa (z^{\mathrm{\prime }})\mathrm{d}{z}^{\mathrm{\prime }})$.

2.2.2 In-plane cooling

Dust grains emission actually occurs in all directions, namely, the radiative flux has also a radial component in addition to the previously described surface cooling. A detailed study on the role of the radial radiative flux or in-plane cooling on the evaluation of the cooling timescale was performed by Miranda & Rafikov (2020a) and Ziampras et al. (2023) based on two-dimensional (2D) vertically integrated disc models. In those papers, an effective cooling time scale was obtained by combining the surface cooling of Eq. (3) and the in-plane cooling obtained considering heat diffusion through the mid-plane by Miranda & Rafikov (2020a). The results obtained underline the importance of considering the in-plane cooling component, suggesting that a study of 3D discs with heat diffusion from the FLD approximation is required to evaluate thermal processes and correctly estimate the cooling time, β_c.

2.3 Cooling time in 3D discs with FLD approximation

When dealing with Eq. (2), heat diffusion is fully taken into account by computing the radiation flux on every grid-cell at any time, t. The cooling time can therefore be derived a posteriori from the timescale over which radiation is diffused over a typical length scale, l_r, $${\beta }_c={\mathrm{\Omega }}_K\frac{l_r^2}{D_r},$$(5)

where D_r is the radiative diffusion coefficient of a grid cell at time, t. In the optically thick limit, the diffusion coefficient associated with the radiation flux vector in Eq. (2) is $$D_r=\frac{16\sigma T^3}{3c_v{\rho }^2\kappa },$$(6)

where we consider E_r = a_rT⁴ and λ = 1/3. Although Eq. (5) overestimates the cooling time in the optically thin regions, we chose to keep this expression since typical discs are optically thick in the planet’s formation region near the midplane.

By considering the disc scale height, H, as the characteristic diffusion length scale, we obtain (see also Miranda & Rafikov 2020a; Ziampras et al. 2023) the following expression, $${\beta }_c={\mathrm{\Omega }}_K\frac{3c_v{\rho }^2\kappa H^2}{16\sigma T^3}.$$(7)

We computed this quantity to compare the results obtained with the energy evolving with FLD approximation to those obtained with Eq. (1). We remark that the cooling time obtained from diffusion Eq. (7) in an optically thick disc is similar to the cooling from the disc’s photosphere, as given in Eqs. (3) and (4).

2.4 Non-ideal magnetic effects

To include MHD effects (with the code FARGO3D), we also considered the induction equation for the evolution of the magnetic field, B, as $${\mathrm{\partial }}_t\mathit{B}=\mathrm{\nabla }\times (\mathit{v}\times \mathit{B}-{\eta }_O\mathit{J}+{\eta }_A\mathit{J}\times {\mathit{e}}_b\times {\mathit{e}}_b).$$(8)

With respect to the pure hydrodynamical case, the equation for the momenta have an additional source term, J × B (Benítez-Llambay & Masset 2016), where J ≡ ∇ × B is the electric current.

The terms η_O and η_A are the Ohmic and ambipolar diffusivities and e_b is the unit vector parallel to the magnetic field line. We chose not to consider the Hall effect for the following reasons: (i) Ohmic resistivity is considered dominant in the region of interest near the midplane where the density is high; in addition, it was recently supposed that it is dominant even at low density values (Hopkins et al. 2024); (ii) to our knowledge, spiral arm propagation in MHD discs has not yet been investigated, so that we chose a relatively simple setting; (iii) finally, from the computational point of view, 3D MHD simulations with embedded planets are extremely expensive and Hall effect can contribute to further decreasing the time step that would make the integration of the system for about 100 planetary orbits prohibitive.

Following Lesur (2021b) and Wafflard-Fernandez & Lesur (2023), we prescribed the Ohmic and ambipolar profiles instead of considering a model for ionisation that would also make the 3D computations extremely expensive. We considered the dimensionless Ohmic Reynolds number, ℛ_O, and Elsasser number, Λ_A, defined at disc midplane via $${\mathrm{\Lambda }}_A(0)\equiv \frac{v_A^2}{{\eta }_A{\mathrm{\Omega }}_k}{\mathcal{R}}_O(0)\equiv \frac{{\mathrm{\Omega }}_K{H}^2}{{\eta }_O},$$

where $v_A=B/\sqrt{\rho {\mu }_0}$ is the Alfvén velocity. According to Lesur (2021b) and Wafflard-Fernandez & Lesur (2023), we can consider the gas to be fully ionised in the corona. This is set via $$\begin{array}{rl}& {\mathrm{\Lambda }}_A(z)=max\{{\mathrm{\Lambda }}_A(0)\exp \left[\frac{z^4}{(\xi H{)}^4}\right],\frac{1}{10}\frac{v_A^2}{c_s^2}\},\\ & {\mathcal{R}}_O(z)={\mathcal{R}}_O(0)\exp \left[\frac{z^4}{(\xi H{)}^4}\right]\left(\frac{\rho (0)}{\rho (z)}\right),\end{array}$$

where ξ quantifies the thickness of the non-ideal layer in units of the pressure scale height H (ξ = 3 in our simulations).

We noticed that the models used to compute the Ohmic and ambipolar coefficients, while taking into account the ionisation of gas and chemical models (Gressel et al. 2015; Béthune et al. 2017; Bai 2017; Lesur 2021b,a), do not provide a clear consensus on the values of the strength of non-ideal effects. This is mainly due to a strong uncertainty about the ionisation rate from cosmic rays (the main source of ionisation below two scale heights) and on the size and abundance of grains in discs that tend to reduce the ionisation fraction.

In the simulations presented in Section 5, we considered ℛ_O = 1 and different values of Λ_A ∈ [0.01 : 10]. This choice is motivated by the fact that values of R_O ~ O(1) and strong (Λ_A ~ 0.01) to weak ambipolar diffusion (Λ_A ~ 10) are plausible (Gressel et al. 2015; Béthune et al. 2017; Bai 2017) in the region of the protoplanetary disc considered in this paper.

2.5 Disc models

The code units are such that G = M_* ≡ M_⊙ = 1 and the unit of distance r₀ is arbitrary. The unit of time is therefore (r₀/au)^3/2/(2π) yr. For the simulations we adopt as the unit of length r₀ = 5.2 au except for the MHD case for which r₀ = 1 au. Planets are considered, in all cases, at a position of r_p = 5.2 au. We call T₀ the orbital period at r = r₀. The simulations’ domain extend radially from r_min to r_max and vertically from the colatitude θ₀ to the midplane at θ = π/2. We consider discs with radial profile of surface density Σ(r) = Σ₀(r/r₀)^−α_Σ. The values of these parameters as well as simulation names are reported in Table 1. The fully radiative simulations are named HD_rad followed by J, S, and 30, for planets of respectively the mass of Jupiter, of Saturn, and of a super-Earth of 30 Earth masses. We considered non-zero α viscosity described by the parameter α_ν for the whole set of hydro-dynamical simulations – except HD_Bae. In the MHD simulations, angular momentum is transported by the magnetic field without the need to prescribe a viscosity (e.g. of turbulent origin).

Other parameters such as the boundary conditions and the numerical resolution do also depend on the simulation (specified in Table 2). In the azimuthal direction we use periodic boundary conditions. The label “damping” refers to the implementation of a wave damping region according to de Val-Borro et al. (2006). The radial grid spacing is logarithmic, the azimuthal one is constant, while the vertical spacing is constant except for the MHD simulations for which we use a nonuniform grid. The reason for this is attributed to magnetic effects, such as winds, developing on a large vertical domain and a non-uniform grid meeting the requirements of having a moderate total number of grid cells together with a suited resolution (~8 grid cells) over a disc scale height above the midplane. Precisely, our grid has constant resolution from the midplane up to one scale height and then it varies quadratically with a maximum resolution ratio of 4. Moreover, in the case of MHD simulations additional boundary conditions are required for the magnetic field and for the electromotive field (EMF). We report the details in Appendix C.

2.6 Discs with embedded planets: Numerical procedure

In all the simulations, we considered a planet on a fixed circular orbit located on the midplane at r_p = 5.2 au (i.e. r_p = r₀ for all the simulations – except the MHD ones, for which r_p = 5.2r₀ with r₀ = 1 au) and azimuth of φ = 0; or, equivalently, at (x_p, y_p, z_p) = (5.2, 0, 0) au. We embedded in each disc model a planet of 20 Earth masses at t = 0 and smoothly increased this value over a time interval of 10 or 20 orbits until the mass until the final mass was reached. The fully radiative and MHD cases require preliminary 2D (r, θ) simulations to reach an equilibrium configuration prior to the planet insertion. This phase is described in the corresponding sections.

The gravitational potential of the planet acting on the disc (Φ_p) is modelled as in Kley et al. (2009) in the fargOCA code: the full gravitational potential is computed for disc elements having distance, d, from the planet larger than a fraction, ϵ, of the Hill radius, while it is smoothed for disc elements with d < r_sm ≡ ϵR_H according to $${\mathrm{\Phi }}_p=\{\begin{array}{ll}-\frac{m_{p}G}{d}& d>r_{\mathrm{s}\mathrm{m}}\\ -\frac{m_{p}G}{d}f\left(\frac{d}{r_{\mathrm{s}\mathrm{m}}}\right)& d\le r_{\mathrm{s}\mathrm{m}}\end{array},$$(9)

with $f\left(\frac{d}{r_{\mathrm{s}\mathrm{m}}}\right)=[{\left(\frac{d}{r_{\mathrm{s}\mathrm{m}}}\right)}^4-2{\left(\frac{d}{r_{\mathrm{s}\mathrm{m}}}\right)}^3+2\frac{d}{r_{\mathrm{s}\mathrm{m}}}]$.

In FARGO3D the potential is smoothed according to $${\mathrm{\Phi }}_p=-\frac{m_{p}G}{\sqrt{(d^2+r_{sm}^2)}}.$$(10)

We recall that the Hill radius is defined as: R_H = r_p(m_p/3M_⊙)^1/3. The values of ϵ are respectively (0.13, 0.2, 0.1) in simulation sets (HD_Bae, HD,MHD).

3 Hydrodynamical case with prescribed cooling

The only 3D simulation, to our knowledge, which shows the formation of multiple rings and gaps inside the orbit of a Jupiter-mass planet is the one in Bae et al. (2016). This paper was mainly devoted to the study of spiral waves induced by a massive planet; however, it offers an important starting point for our study. The aim of this section is to produce a simulation similar to the one presented in Section 4.2 of Bae et al. (2016) and to check whether the density maxima that form inside the Jupiter’s orbit correspond to pressure bumps with the ability of stopping dust from drifting inwards.

Although we did not consider dust dynamics in our simulations, we noticed that dust velocity would be affected by the gas drag. Thus, by introducing a radial dust velocity, u_r, we have (Nakagawa et al. 1986; Takeuchi & Lin 2002): $$u_r=\frac{v_r}{1+S_t^2}-2S_t\eta v_K,$$(11)

where S_t is the Stokes number, namely, the value of the friction timescale in units of the orbital frequency, and η is the fraction of the Keplerian velocity, v_K, corresponding to the headwind experienced by dust particles, $$\eta =-\frac{1}{2}{\left(\frac{H}{r}\right)}^2\frac{\mathrm{\partial }\log P}{\mathrm{\partial }\log r}.$$(12)

Neglecting the gas radial velocity, v_r, which is expected to be low, a negative derivative of the pressure with respect to r corresponds to η > 0 so that the dust moves toward the star; whereas for a positive pressure gradient (η < 0) the situation is reversed. A dust trap is located at the pressure maximum where η = 0.

In locally isothermal discs, $P=c_s^2\rho$; whereas for more generic equations of state (EoS), we have P = (γ − 1)e. The disc we considered in the HD_Bae simulation was initialised with a locally isothermal pressure, giving a value of η in the unperturbed disc of 0.003 (considering the parameters in Table 1).

In our HD_Bae simulation, we considered a non-viscous disc, along with Eq. (1) for the time evolution of the internal energy density as in Bae et al. (2016). However, we did not implement the computation of ${\beta }_c^{surf}$ from Eq. (4) and we considered a constant prescribed β_c corresponding to ${\beta }_c\sim 100\sim {\beta }_c^H$ (see Eq. (4)) obtained by Bae et al. (2016) at the midplane at the planet location.

In Fig. 1 we show the midplane density in the (r, φ) plane when the planet has just reached its final mass. Multiple spiral arms are clearly shown (cf. Fig. 6 of Bae et al. 2016).

We considered the same snapshots of Fig. 5 of Bae et al. (2016). We show in Fig. 2, the azimuthally averaged surface density profile (top panel) and the η profile (bottom panel). In addition to the first bump at the gap’s inner edge (at about r/r_p = 0.7), another density maximum associated to a pressure bump is located at r/r_p = 0.54 and a third maximum at about 0.4, where η gets close to zero. Two additional pressure maxima appear inside r/r_p = 0.37, but they do not correspond to rings. Specifically, in Fig. 3, we see the perturbed surface density plotted in the same snapshots of Fig. 2, indicating that the averaged maxima correspond to rings extending over the full disc in azimuth (excluding the two innermost maxima, which correspond to vortices).

In conclusion, a Jupiter-mass planet can form multiple rings and gaps inside its orbit in a non-viscous disc. Such rings are effective pressure bumps and can potentially stop dust from drifting towards the star. We stress the term ‘potentially’ because η is computed from the azimuthally averaged pressure and the actual response of dust to non-axisymmetric features (e.g. spiral waves) ought to be investigated in a future work.

Fig. 1 HDBae simulation. 2D distribution in the (r, φ) plane of the mid-plane density, ρ, normalised over the azimuthally averaged value < ρ > at t = 10 orbits when the planet has just reached a Jovian mass. For the purposes of comparison with Fig. 6 of Bae et al. (2016), the y axis is plotted using a logarithmic scale.

Fig. 2 Top panel: azimuthally averaged surface density profile for simulation HDBae. Bottom panel: azimuthally averaged η parameter. The planet is kept on a fixed circular orbit at r/rp = 1 and the ticks on the x axis (as well as the vertical black lines) indicate the values of η ~ 0 with a positive slope. The horizontal black line corresponds to η = 0. The damping region extends radially from rmin/r0 to the dot dashed line.

4 Hydrodynamical case with complete treatment of radiative transport

The aim of this section is to study the propagation of density waves with a self-consistent treatment of radiation transport within the FLD approximation. With respect to the study presented in Bae et al. (2016) (and revisited in Section 3) we also considered a non-zero viscosity. In fact, discs have a low (but non-zero) bulk viscosity, as evidenced by observations (Pinte et al. 2016; Villenave et al. 2022) and supported by mechanisms such as the vertical shear instability (Nelson et al. 2013) leading to low viscosity in dead zones of the discs (where viscosity is not sustained by magnetorotational instability). The value of the alpha parameter α_ν is expected to be lower than 10⁻³ and possibly take values of the order of 10⁻⁴. In the following, we consider α_ν = 10⁻⁴.

Before inserting the planet, we ran a 2D axisymmetric (r, θ) simulation until the disc reaches thermal equilibrium. We notice that the main source of heating in a low-viscosity disc is the stellar heating. To recover the temperature profile of a stellar irradiated disc, we ran a 2D (r, θ) simulation, including the heating from a Solar-type star and using the temperature values of the resulting equilibrium disc as boundary for our 3D disc at the disc’s surface (as described in Lega et al. 2024). The equilibrium disc is flared (with flaring index 2/7) and has an aspect ratio of 0.03 at r_p. The density slope, α_Σ, (see Table 1) was chosen on the basis of a disc with constant star accretion rate (Lega et al. 2015). The mass flow carried by the disc was set to be $\dot{M}\sim {10}^{-9}{M}_{\odot }/\mathrm{y}\mathrm{r}$, which is about one order of magnitude smaller than the typical accretion rate observed in young stars (Hartmann et al. 1998; Manara et al. 2016). We recovered typical star accretion rates by including non-ideal MHD effects (described in Section 5).

Figure 4 shows the values of the equilibrium disc cooling time, density, and opacity in the (R, Z) plane. When looking at values of β_c in Fig. 4 (left panel), the disc can be divided in two main regions based on the iso-contour line, logβ_c = 1. Specifically, inside the contour of logβ_c = 1 (Fig. 4, left panel), the cooling time values correspond to the regime of density wave propagation referred to as ‘adiabatic’ in Miranda & Rafikov (2020a) (see their Fig. 4). When inserting a planet in this region, we expect to find a gap at the planetary orbit and multiple rings and gaps in the inner region. Near the logβ_c = 1 contour, density waves are expected to be radiatively damped (see Ziampras et al. 2023). Finally, for logβ_c ≪ 1, the disc is in the locally isothermal regime. This regime is found near the disc surface, far from the planet’s forming region. We remark that the β_c values are expected to change with time according to the perturbation introduced by the planet; thus, wave propagation regimes can be modified accordingly.

Fig. 3 Perturbed surface density of the HDBae simulation at the same snapshots represented in Fig. 2. The dashed lines indicate the radial position of the pressure bumps, i.e. the values of r for which η = 0 in Fig. 2. The dot-dashed line indicates the boundary of the damping region.

Fig. 4 Disc cooling time, density, and opacity in the (R, Z) ≡ (r sin(θ), r cos(θ)) plane of the disc in thermal equilibrium before the insertion of the planet.

4.1 Jupiter-mass planet

Starting from the (R, Z) disc in thermal equilibrium, we expanded it in azimuth over the interval [−π, π] and restarted our simulation with an embedded planet. We ran the HD_radJ simulation over 150 orbits at the planet location and obtained the azimuthally averaged values of β_c, as well as the volume density and the optical thickness at the end of the integration (see Fig. 5). At the point, the planet has carved a gap and deeply modified the disc structure. However, the values of β_c inside the planetary orbits are still compatible with the ‘adiabatic’ regime of wave propagation; thus, we would expect to find multiple gap and rings in the inner disc region (r < 4 au).

In Fig. 6, we show the perturbed surface density of the HD_radJ simulation at different snapshots indicated in orbital period at the planet location: in the left panel, at t = 20 orbits, we clearly see the gap at the planet’s orbit and a primary ring of gas at the inner and outer gap edges (with the density maxima, respectively, at r ~ 4 and ~ 6.7 au) as well as the primary spiral arm. The secondary and a tertiary arms are clearly present inside the planetary orbit launched, respectively, at (r, φ) ~ (3.8, −π) and (r, φ) ~ (3.8, −1.8), propagating in the inner disc. The presence of multiple spiral arms is a necessary condition to have multiple gaps (Miranda & Rafikov 2019, 2020b). As for the primary gap at the planet’s orbit, multiple gaps form when density waves steepen into shocks depositing angular momentum into the gas. A second and a third ring at r ~ 3.3 and 2 au, respectively, are clearly formed at t = 100 (Fig. 6, second panel, and Fig. 7, top panel). Over the total integration time of 150 orbits, the primary gap is still evolving and rings of density maxima slightly move inwards and outwards as indicated by the vertical dotted lines in Fig. 6. The density contrast of the innermost ring increases with time while the second ring at r ~ 3.3 partially merges with the primary inner ring as the former eventually moves inwards.

From the computation of η (Fig. 7, bottom panel), we can see there is a pressure bump at the innermost ring at about 2 au, so that a reservoir of dust would be assumed to be trapped there, while the primary and secondary rings are marginally separated at t = 150 orbits. The pressure bump at about 3.3 au should stop dust from drifting inwards, thereby making a distinct reservoir with respect to the one centered at 2 au. The 3D fully radiative simulations require about a longer computational time by a factor of 10 with respect to isothermal or adiabatic simulations with controlled cooling². The important result here is that distinct dust reservoirs can be formed in a realistic scenario, where thermal processes are computed in a self-consistent way.

Fig. 5 Same as Fig. 4 for azimuthally averaged quantities for the 3D simulation HDradJ with an embedded Jupiter-mass planet. The snapshot corresponds to the end of the integration at 150 orbital periods at the planet location (5.2 au).

Fig. 6 Perturbed surface density of the HDradJ simulation at different snapshot indicated in orbital period at the planet location. The letters in the left and central panels indicate the primary (P), secondary (S), and tertiary (T) spiral arms. The dashed lines indicate the position the radial position of the pressure bumps, i.e. the values of r for which η = 0, with a positive slope, shown in Fig. 7.

Fig. 7 Top panel: azimuthally averaged surface density profile for simulation HDradJ at the same snapshot represented in Fig. 6. Bottom panel: azimuthally averaged η values. The dashed line indicates the boundary of the damping region.

4.2 Comparison with various EoSs

In this subsection, we consider simulations with an embedded Jupiter-mass planet (as in the HD_radJ case, but using a different EoS) to compare the results obtained in HD_radJ with different wave-damping regimes (see Miranda & Rafikov 2020a; Ziampras et al. 2023). In the following, we consider discs with respectively an isothermal EoS (HD_iso), an adiabatic EoS with temperature relaxation to the initial value as in Eq. (1) (HD_βc with β_c = 1 and β_c = 100) and an adiabatic case without thermal relaxation (HD_adia).

We integrated the new simulations on longer timescales with respect to HD_radJ taking advantage of shorter running times³ to check the persistence of rings and gap structures. Moreover, the disc domain is extended further in to investigate the possible formation of rings inside the one at 2 au obtained in the case of HD_radJ; in addition, we sought to confirm that this inner-most ring at 2 au is not affected by the close inner disc boundary (which is at 1.5 au in HD_radJ).

We show in Fig. 8 the perturbed surface density at t = 150 orbits at the planet location (top panels) and at t = 500 orbits (bottom panels) for all these cases; for comparison, we also inserted the results obtained before for the HD_radJ case at t = 150 orbits between the two HD_βc cases. As for the 2D case (Ziampras et al. 2020), the locally isothermal disc is the most efficient in forming multiple rings and gaps: three in our disc, which are also very contrasted (see Fig. 9, top panel), and are effective dust barriers based on the radial profile of η (Fig. 9, bottom panel). We confirm that we observed spiral arms extending further in, but no additional ring was observed inside r = 2 au.

At t = 500 orbits, the rings appear more contrasted with respect to t = 150 orbits, which is a good indication for their persistence. In all the cases, an innermost ring is formed close to r = 2 au, except in the case with β_c = 1. According to the results obtained for 2D discs, this is the case where wave damping appears to be the most efficient (when no radiative effects are taken into account; see Ziampras et al. 2023).

We note the similarity between the fully radiative case and simulation HD_βc with β_c = 100 (i.e with a cooling time consistent with typical cooling times inside the Jupiter orbit; as shown in Fig. 5). We can also state that the ring at r = 2 au is not affected by the close radial boundary, since a very similar structure has been obtained in discs extending further in.

5 MHD simulation with prescribed cooling

In the previous section, we describe how the opening of secondary gaps by a Jupiter-mass planet is observed in low-viscosity 3D discs with a self-consistent treatment of the diffusion of heat. We consider the formation history of Jupiter by analyzing smaller masses in Section 6. In this section, we describe how low-viscosity discs transport gas towards the star at very low rates with respect to typical accretion rates observed in young stars (Hartmann et al. 1998; Manara et al. 2016). Instead, in magnetised discs embedded in warm fully ionised corona, gas transport at disc surface layers (at about two to three disc scale heights) can provide the required accretion flow onto the star, which we did not include in our simulation HD_radJ. In fact, gas can flow inwards because of angular momentum removal by disc winds (Béthune et al. 2017; Lesur 2021b).

In the simulations presented here, we did not solve the radiation field; instead, we prescribed a warm corona via a steep vertical increase of the disc midplane temperature from three to five pressure scale heights. Precisely, similarly to Béthune et al. (2017), the disc is initialised with a temperature ratio between the corona and the disc T_cd = 6, as follows, $$T(\tilde{\theta })=1+(T_{cd}-1){\left(\frac{1}{1+\exp [-6(|\tilde{\theta }|-4h)/h]}\right)}^2,$$(13)

where $\tilde{\theta }\equiv \pi /2-\theta$ is the latitude, h is the disc aspect ratio (constant as our MHD discs are considered non-flared; see Appendix A). To this aim, we assumed the following function for the sound speed, $$c_s(r,\theta )=\sqrt{\frac{G{M}_{\ast }}{r}{h}^{2}T(\tilde{\theta })}.$$(14)

Considering the equations for vertical and rotational equilibrium, we provide in Appendix A the initial volume density and the initial azimuthal velocity we used in our MHD simulations. After the disc initialisation, we increased T_cd to 10 in the time interval [0 : 30] (code units).

The initial magnetic field is vertical B = B_z. Although there is a high level of uncertainty about the topology of the magnetic field of protoplanetary discs, their intensity is known to be weak (Lesur 2021b,a) in the sense that the midplane ratio of thermal over magnetic pressure: $\beta \equiv 2{\mu }_{0}P/B_z^2\gg 1$. Therefore, we considered discs with initial midplane values of β = 10⁶, 10⁵, 5 × 10³ and computed the radial profile of the magnetic field at midplane from β and from the midplane pressure, P = ρ(r, 0)c_s(0)².

We did not solve the energy equation with the FLD approach (Eq. (2)) as in simulation HD_radJ since this would require extremely long integration times for 3D MHD simulations. However, as we show in Fig. 8, when relaxing the internal energy to the initial temperature profile according to Eq. (1) on a cooling time β_c = 100, we obtain results that are very similar to the fully radiative case. Therefore, we considered Eq. (1) in our MHD simulations with cooling time: β_c = 100.

Before inserting the planet, we ran 2D (r, z) simulations for a few hundreds orbits at r = 1 au. At t = 300T₀ (T₀ = 2π yr) we show (Fig. 10, top panel) the density distribution on the (r, Z) plane for the case β = 5 × 10³ and Λ_A = 1 (fiducial simulation). In Fig. 10 (top panel) we see that in the regime of Ohmic diffusion dominating over the Ambipolar one, there is no evidence of spontaneous generation. In other words, in the absence of the planet, the density distribution does not show rings and gaps (see e.g. Béthune et al. 2017). Therefore, when inserting a planet, any gap and ring formation will be a consequence of the perturbation induced by the planet itself. We also computed the mass flux as a function of r as $$\dot{M}=2({\dot{M}}_d+{\dot{M}}_w),$$(15)

where the factor of 2 takes into account the full disc. The quantities ${\dot{M}}_d$ and ${\dot{M}}_w$ are the mass flux of the disc region up to 2H and the corona from 2H to the disc surface, respectively. These quantities are computed from $${\dot{M}}_d={\int }_0^{2\pi }{\int }_{{\theta }_{2H}}^{\pi /2}\rho v_r{r}^{2}sin(\theta )\mathrm{d}\theta \mathrm{d}\phi .$$(16)

Here, θ_2H is the colatitude at two disc scale heights and then we have $${\dot{M}}_w={\int }_0^{2\pi }{\int }_{{\theta }_0}^{{\theta }_{2H}}\rho v_r{r}^{2}sin(\theta )\mathrm{d}\theta \mathrm{d}\phi ,$$(17)

with θ(0) the value of the colatitude at the disc surface. In the bottom panel of Fig. 10, we plot this quantity averaged in time from 100 to 300T₀ for our fiducial simulation. Negative values correspond to flux directed towards the star, while positive ones correspond to gas outflow. The corona appears to be dominated by disc out flowing wind while the disc accretes with values of $\dot{M}\sim {10}^{-7}{M}_{\odot }/\mathrm{y}$, similar to typical accretion rate observed in young stars (Hartmann et al. 1998; Manara et al. 2016). We provide in Appendix B an analysis of the mass flux from angular momentum conservation arguments.

We let the planet growing for t = 10 orbital periods at the planet location and continued our simulation over 100 orbits. In Fig. 11, we show the perturbed surface density of our MHD runs at t = 100 orbits at r_p = 5.2 au.

From left to right, we show the increased strength of the magnetic field and for a given value of the magnetic field, we can also increase the ambipolar diffusivity depending on the results obtained for multiple ring formation. For the lowest value of the magnetic field (β = 10⁶, left panel), we see the formation of a secondary ring t about r = 2 au. In this case, we considered Λ_A = 10, which corresponds to a very low effect of the non-ideal ambipolar diffusion effect. Keeping this same value of Λ_A, and increasing the strength of the magnetic field (β = 10⁵, second panel), we can scarcely identify the formation of a secondary ring. Instead, by increasing the ambipolar coefficient (Λ_A = 1, third panel), we observe a significant enhancement of the secondary ring. By further increasing the magnetic field (β = 5 × 10³, right panels), we observe that the inner planet-induced spiral arms can be damped; in this case, no secondary ring would ever form at Λ_A = 10 (not shown) or Λ_A = 1 (fourth panel), while a weakly contrasted secondary ring appears for Λ_A = 0.1 (fifth panel). Finally, a well-contrasted secondary ring appears in the right-most panel of Fig. 11 for Λ_A = 0.01.

For the largest value of β, the magnetic field is immaterial and the formation of secondary rings closely resembles that occurring in purely hydrodynamical discs, regardless of the Ohmic and ambipolar diffusion coefficients. As the magnetic field increases (for smaller values of β), MHD effects come into play and the magnitude of ambipolar diffusion is found to have a strong impact on the existence of secondary rings.

This effect is similar to the one observed in the formation of zonal flow (Béthune et al. 2017; Riols et al. 2020) in disc regions where Ambipolar diffusion is the dominant non-ideal effect. In particular, spontaneous accumulations of gas are observed with anti-correlated concentrations of the vertical component of the magnetic field B_z. It has been shown (Béthune et al. 2017) that ambipolar diffusion is responsible for the accumulation of B_z. Therefore, if one or multiple gaps forms in a non-spontaneous way (e.g. because a massive planet perturbs the gas), then B_z will accumulate in the gaps (Wafflard-Fernandez & Lesur 2023) and lead to gas accumulation in rings, provided that the ambipolar diffusion is strong enough.

In Fig. 12, we show the azimuthally averaged surface density (top panel) and the η values (bottom panel) as a function of r at the same snapshot of Fig. 11. The plots correspond to the simulations presented in Fig. 11 for which a clearly contrasted secondary ring is formed (i.e. for β = 10⁶ with Λ_A = 10, for β = 10⁵ with Λ_A = 1, and for β = 5 × 10³ with Λ_A = 0.01).

The values of η for both the primary and secondary rings at r 2 au reach values close to zero, although we do not observe a clear transition from negative to positive η, which is necessary for dust-trapping. A longer integration time would be required to verify whether dust-trapping is possible in these rings. Nevertheless, we consider very interesting that rings form for values of non-ideal parameters that are plausible in the disc region inside the Jupiter orbit.

Fig. 8 Perturbed surface density at t = 150 orbits at the planet location (top panels) and at t = 500 (bottom panels) orbits. From left to right: simulation sets HDiso, HDβc with βc = 1, HDrad, and HDβc with βc = 100, HDadia. We observe multiple rings formation that depend on the EoS. Rings are persistent on longer integration times. For the radiative case, we report the last snapshot, i.e. the one at 150 orbits in the bottom panel.

Fig. 9 Azimuthally averaged surface density profile (top) and the η values (bottom) for simulation sets HDiso, HDβc, and HDadia at t = 500 orbits at the planet’s location.

Fig. 10 Top panel: density distribution in the (r, Z) plane obtained for our simulation MHDβ5e3 with ΛA 1 at t = 300T0. Bottom panel: computation of $\dot{M}(r,\theta )$ (Eq. (15)), averaged over the time interval [100 : 300]T0.

Fig. 11 Perturbed surface density of the MHD simulations. From left to right: increase in the strength of the magnetic field. In addition, for a given value of the magnetic field, we can also increase the Ambipolar coefficient value. For the lowest value of the magnetic field (left panel), we see the formation of a secondary ring at r = 2 au with our lower value of ambipolar diffusion (ΛA = 10); thus, we do not run additional simulations in this case. In the case of βm = 105 (second panel and third panels) we observe a secondary ring for ΛA = 10; however, the formation of a secondary well contrasted ring requires a higher ambipolar diffusion strength (ΛA = 1, third panel). Finally, for βm = 5 × 103, there is no secondary ring formed for ΛA = 10 (not shown) nor for ΛA = 1 (fourth panel). We observe a poorly contrasted secondary ring for ΛA = 0.1 (fifth panel). Finally, we observe a well contrasted secondary ring for ΛA = 0.01 (right-most panel).

Fig. 12 Azimuthally averaged surface density profile (top) and azimuthally averaged η values (bottom) for MHD simulations presented in Fig. 11. Here, a well-contrasted secondary ring is formed. The planet is kept on a fixed circular orbit at rp = 5.2 au and the ticks on the x axis (as well as the vertical black lines) indicate the values of η ~ 0 with positive slope. The horizontal black line corresponds to η = 0. The damping region extends radially from rmin to the dotted line.

6 Minimum planet mass for forming dust-trapping rings

It is well known (Lambrechts et al. 2014; Bitsch et al. 2018) that when a giant planet core reaches few tens of Earth masses, the so-called pebble isolation mass is reached and dust (at suited Stokes numbers) is blocked at the outer gap’s edge opened by the core. Moreover, in low-viscosity discs even cores of less than 20 Earth masses can open small gaps lowering the limit of the pebble isolation mass (see Bitsch et al. 2018). On the other hand, the dust inside the planet’s orbit drift towards the star and is eventually stopped at the inner (close to 1 au) pressure bump due do gas removal by magnetised winds (Suzuki et al. 2010; Ogihara et al. 2018).

We show in this paper that a Jupiter-mass planet has the ability to form multiple rings, according to disc thermal properties and to non-ideal magnetic diffusivities. However, such rings can be effective dust’s reservoirs if they form early on in the Jupiter formation history, possibly when Jupiter core has reached or slightly overcome the pebble isolation mass. Therefore, in this section, we describe the simulations we performed for a planet of a Saturn mass (simulation HD_radS) and for a super-Earth of 30 Earth masses (simulation HD_rad30) to check their ability to form multiple rings.

The super-Earth opens a gap and forms a secondary ring (at ~3.8 au). which increases in density contrast on long integration times (see Fig. 13). A third density maximum appears which, however, is not associated to a local pressure maximum since it has η > 0. Instead, three rings appear for Saturn at 150 orbits and two for Jupiter, as previously discussed. The values of r for which η = 0 move inwards when we are increasing the planet’s mass.

From Fig. 13, we can propose the following evolutionary picture. When a planet reaches the pebble isolation mass, it starts trapping dust at two locations near the inner edge of its gap. As the planet grows in mass, the locations of these two rings shift inwards, together with their dust. Near the mass of Saturn, a third ring appears near 3 au, trapping the dust that had not been trapped in the two other rings (unless the planet grows so slowly that all the untrapped dust has already drifted into the star). As the planet grows further, the two initial rings merge and the last one moves from ~3 to ~2 au. Obviously, this picture does not take into account that the planet is migrating at the same time. However, in the case of a proto-Jupiter, its migration may have been stalled or even reversed by the appearance of a proto-Saturn (Walsh et al. 2011).

Fig. 13 Top panel: azimuthally averaged surface density profile for simulations Hrad30, HradS and HradJ at 150 orbital periods at the planet location. Bottom panel: azimuthally averaged η values. Colour markers at η = 0 with a positive slope identify the position of the pressure bumps on the top panel.

7 Conclusion

Cosmochemical observations provide evidence that millimetresize dust remained trapped in the inner part of the disc for a timescale of at least a million years. Multiple trapping sites would have had to exist to explain the formation of four families of chondrites, which are chemically and isotopically distinct from each other.

In this work, we have explored the possibility that Jupiter formed local maxima in the azimuthally averaged pressure in the inner part of the protoplanetary disc (inwards of its orbit). These averaged pressure bumps are potential sites of dust-trapping, but the actual response of dust to the non-axisymmetric features that characterise these bumps remains to be investigated.

We have shown that in a low-viscosity disc, the opening of secondary gaps and density bumps is robust. This was also observed in our 3D simulations when we self-consistently simulated the diffusion of heat, giving a cooling timescale of about 100 orbital times in the inner disc.

When considering a strong magnetic field, we observe that the inner planet-induced spiral arm can be damped and a secondary ring forms only if ambipolar diffusion is strong enough. Further work is warranted to assess the origin of this damping. We note that it could arise from dissipation, when the diffusion coefficients are at adequate levels, or from a removal of the angular momentum flux by the wind. We remark that self-consistent disc ionisation models show that the appropriate parameters for non-ideal MHD are consistent with those for which we observed the effective formation of secondary gaps and rings in our simulations.

A planet seems to be able to open secondary gaps and trap dust in rings when it reaches a mass close to the pebble isolation mass. Thus, as it cuts the flow of dust from the outer disc across its orbit, it also traps part of the dust in rings located inwards of its orbit. As the planet grows, these rings move away from the planet and can merge with each other, while a new ring can be formed farther away. Coupled with planet migration, these results provide a new view of the non-trivial evolution of dust in the inner Solar System, which may be consistent with cosmochemical constraints on dust preservation and confinement on extended timescales.

Acknowledgements

A.M. is grateful for support from the ERC advanced grant HolyEarth N. 101019380. This work was granted access to the HPC resources of IDRIS and CINES under the allocation A0160407233 made by GENCI. F.M. acknowledges support from UNAM’s grant PAPIIT 107723, UNAM’s DGAPA PASPA program and the Laboratoire Lagrange at Observatoire de la Côte d’Azur for hospitality during a one-year sabbatical stay. E.L. wish to thank Alain Miniussi for maintenance and re-factorisation of the code fargOCA. This work was supported by the French government through the France 2030 investment plan managed by the National Research Agency (ANR), as part of the Initiative of Excellence Université Côte d’Azur under reference number ANR-15-IDEX-01.

Appendix A Setting initial conditions for MHD simulations

We describe the method we used to determine the initial conditions for our MHD simulations described in Section 5 in order to model discs with a warm corona. To this aim we start with a corona mildly hotter than the disc’s midplane, and determine the vertical density and rotational velocity profiles that correspond to strict hydrostatic and centrifugal equilibria. We then progressively increase the corona’s temperature until we reach the desired corona to disc temperature ratio.

We assume that the sound speed is inversely proportional to the spherical radius and has a separable dependence on r and θ: $$c_s^2(r)={\left(c_{s0}^0\right)}^2{\left(\frac{r}{r_0}\right)}^{-1}g(\theta ),$$(A.1)

where g(θ) is an even function with respect to π/2 and g(π/2) = 1 and r₀ is an arbitrary reference radius where the sound speed, at the midplane, is $c_{s0}^0$. We have considered for g(θ) the function $T(\tilde{\theta })$ defined in Eq. 13. We consider a constant disc aspect ratio: $$h(r)=\frac{c_{s0}(r)}{v_K(r)}=\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t},$$(A.2)

where $v_K(r)=\sqrt{\frac{G{M}_{\ast }}{r}}$ is the circular Keplerian velocity at mid-plane, which implies that $c_{s0}^0\equiv h{v}_k\left(r_0\right)$.

The equations for rotational and vertical equilibrium in spherical coordinates are $$-\frac{{\mathrm{\partial }}_r\left({\rho }_0{c}_s^2\right)}{{\rho }_0}+\frac{v_{\varphi }^2}{r}-\frac{G{M}_{\ast }}{r^2}=0,$$(A.3) $$-\frac{1}{r}\frac{{\mathrm{\partial }}_{\theta }\left({\rho }_0{c}_s^2\right)}{{\rho }_0}+\frac{v_{\varphi }^2}{r}\cot \theta =0.$$(A.4)

Here, we have ρ₀ and $c_{s0}(r)=c_{s0}^0\times {\left(\frac{r_0}{r}\right)}^{1/2}$ respectively the density and the sound speed at the midplane at radius r. We introduce notations similar to those of Masset & Benítez-Llambay (2016) and define $$L=\log {\rho }_0,v=\log \left(\frac{r}{r_0}\right),u=-\log (\sin \theta ).$$

Equations (A.3) and (A.4) become $$-c_s^2{\mathrm{\partial }}_{v}L+c_s^2+v_{\varphi }^2-\frac{G{M}_{\ast }}{r}=0,$$(A.5) $$c_{s0}^2(r)\frac{{\mathrm{\partial }}_u[{\rho }_{0}g(\theta )]}{{\rho }_0}+v_{\varphi }^2=0.$$(A.6)

Introducing $m\equiv \frac{v_{\varphi }^2}{c_s^2}$, we get: $$-{\mathrm{\partial }}_{v}L+1+m-\frac{G{M}_{\ast }}{r{c}_s^2}=0,$$(A.7) $$c_{s0}^{2}g(u){\mathrm{\partial }}_{u}L+c_{s0}^2{\mathrm{\partial }}_{u}g(u)+v_{\varphi }^2=0.$$(A.8)

As per our assumption of a constant aspect ratio, $K\equiv \frac{G{M}_{\ast }}{r{c}_{s0}^2}$ is a constant and the equations above can be recast as: $$-{\mathrm{\partial }}_{v}L+1+m-\frac{K}{g(u)}=0,$$(A.9) $${\mathrm{\partial }}_{u}L+\frac{{\mathrm{\partial }}_{u}g(u)}{g(u)}+m=0.$$(A.10)

Define $G(u)=\frac{1}{g(u)}$, and L′ = L + log g(u), giving: $$-{\mathrm{\partial }}_v{L}^{\mathrm{\prime }}+1+m-KG(u)=0,$$(A.11) $${\mathrm{\partial }}_u{L}^{\mathrm{\prime }}+m=0.$$(A.12)

Let m′ = m − KG(u), thus, we have $$-{\mathrm{\partial }}_v{L}^{\mathrm{\prime }}+1+m^{\mathrm{\prime }}=0,$$(A.13) $${\mathrm{\partial }}_u{L}^{\mathrm{\prime }}+m^{\mathrm{\prime }}+KG(u)=0.$$(A.14)

The values of m′ and L′ are determined using a method of characteristics. We introduce the following combinations of the coordinates u and v: $$s=u+v,$$(A.15) $$s^{\mathrm{\prime }}=u-v,$$(A.16)

so $$u=\frac{1}{2}(s+s^{\mathrm{\prime }}),$$(A.17) $$v=\frac{1}{2}(s-s^{\mathrm{\prime }}).$$(A.18)

Using Eqs. A.13 and A.14, we derive $${\mathrm{\partial }}_s{m}^{\mathrm{\prime }}=\frac{1}{2}({\mathrm{\partial }}_u{m}^{\mathrm{\prime }}+{\mathrm{\partial }}_v{m}^{\mathrm{\prime }})=0$$

Thus, m′ depends only on s′, and its dependence on s′ can be determined at the midplane. Considering a constant h and a midplane density, ${\rho }_{mid}=\frac{{\mathrm{\Sigma }}_0{\left(r/r_0\right)}^{-{\alpha }_{\mathrm{\Sigma }}-1}}{\sqrt{2\pi }h}$, we obtain from Eq. A.13: $$\xi +1+m^{\mathrm{\prime }}=0\Rightarrow m^{\mathrm{\prime }}=-\xi -1,$$

with ξ = α_Σ + 1 Hence, $$m=h^{-2}G(u)-1-\xi .$$(A.19)

From Eq. (A.14), we get $${\mathrm{\partial }}_u{L}^{\mathrm{\prime }}+h^{-2}G(u)-1-\xi =0.$$(A.20)

Letting $W(u)={\int }_0^{u}G\left(u^{\mathrm{\prime }}\right)d{u}^{\mathrm{\prime }}$, we have $$L^{\mathrm{\prime }}(u)=L^{\mathrm{\prime }}(0)-h^{-2}W(u)+(1+\xi )u,$$

hence $$L(u)=-\log g(u)+L(0)-h^{-2}W(u)+(1+\xi )u.$$

Thus, the density profile is $$\rho (u)=\frac{1}{g(u)}{\rho }_{\text{mid\hspace{0.17em}}}\exp [-h^{-2}W(u)+(1+\xi )u],$$(A.21) $$={\rho }_{00}{\left(\frac{r}{r_0}\right)}^{-\xi }(\sin \theta {)}^{-1-\xi }G(u)\exp [-h^{-2}W(u)],$$(A.22)

with ${\rho }_{00}=\frac{{\mathrm{\Sigma }}_0}{\sqrt{2\pi }h}$. Finally, from the expression for m (Eq. A.19), the azimuthal velocity is $$v_{\varphi }^2={\left(c_{s0}\right)}^2{h}^{-2}-(1+\xi )c_s^2\text{,\hspace{0.17em}}$$(A.23) $$={\left(c_{s0}^0\right)}^2\left(\frac{r_0}{r}\right)\frac{G{M}_{\ast }}{r_0{\left(c_{s0}^0\right)}^2}-(1+\xi )c_s^2,$$(A.24) $$=\frac{G{M}_{\ast }}{r}-(1+\xi )c_s^2.$$(A.25)

Appendix B MHD simulations diagnostics

In the following, we provide some diagnostics on our MHD simulations based on the study of mass transport through the disc in terms or mass flux or mass accretion (or mass loss through the wind) rates.

The radial mass flux is generated by torques that determine a loss of angular momentum. We follow Béthune et al. (2017) and compute quantities defined in their Eqs. 8 and 9 and reported here:

• The azimuthal average of a quantity X is denoted as $$<X{>}_{\phi }=\frac{1}{\mathrm{\Delta }\phi }{\int }_0^{\mathrm{\Delta }\phi }X(r,\theta ,\varphi ,t)d\varphi$$

• Radial profiles are obtained by vertical integration in the disc domain, $$<X{>}_{\phi ,\theta }=\frac{1}{2H(r)}{\int }_{-H(r)}^{H(r)}<X{>}_{\phi }(r,\theta ,t)d\theta$$

• A density weighted quantity is denoted as $$<X{>}_{\rho }=\frac{<\rho X{>}_{\phi ,\theta }}{<\rho {>}_{\phi ,\theta }}$$

• The fluctuating Reynold stress tensor is defined by: $$\mathcal{R}\equiv \rho \tilde{\mathbf{v}}\otimes \tilde{\mathbf{v}}$$

  where $\tilde{\mathbf{v}}=\mathbf{v}-<v{>}_{\rho }$.

• The Maxwell stress tensor is: $$\mathcal{M}\equiv -\mathbf{B}\otimes \mathbf{B}$$

  and we call 𝒯 ≡ ℛ + ℳ

• We use the same notation as Shakura & Sunyaev (1973) to define α values by normalising with vertically integrated pressure: $${\alpha }^{\mathcal{R}}(r,\theta ,t)=\frac{<\mathcal{R}{>}_{\phi }}{<P{>}_{\phi ,\theta }}$$(B.1)

  and similarly for values associated to the magnetic field, denoted α^ℳ.

Fig. B.1 Computation of the α values from Eq.B.1 averaged over time for the fiducial MHD simulation with β = 5 103 and ΛA = 1.

The radial mass flux can be deduced from the conservation of angular momentum, provided Σ(r) = 2H(r) < ρ >_φ,θ: $$\{\begin{array}{ll}\mathrm{\Sigma }(r)<v_r{>}_{\rho }& \equiv {\tau }_r+{\tau }_z\\ {\tau }_r& \simeq -\frac{1}{r({\mathrm{\partial }}_r<r{v}_{\phi }{>}_{\phi ,\theta })}{\mathrm{\partial }}_r(2r^{2}H(r)<{\mathcal{T}}_{r,\phi }{>}_{\phi ,\theta })\\ {\tau }_z& \simeq -\frac{r}{[{\mathrm{\partial }}_r<r{v}_{\phi }{>}_{\phi ,\theta }]}{(<{\mathcal{T}}_{\phi ,\theta }{>}_{\phi })}_{-H}^{+H}\end{array}$$(B.2)

where τ_r gives the mass flux due to the radial component of the angular momentum and τ_z gives the mass transport due to angular momentum extracted at disc surfaces z = ±H

In Figure B.1, we clearly see that the α values associated to the magnetic field dominate with respect to the Reynold α values. From Fig. B.2, we can appreciate that the 1.h.s. and the sum of the τ component of the r.h.s. of the equation follow nicely the same trend both when averaged on time and on space in the interval r ∈ [2 : 8]r₀. Moreover, in this radial interval there is no contribution to the mass flux from the radial transport of the angular momentum. Precisely, let’s notice (Fig.B.1), that $<{\alpha }_{r\phi }^{\mathcal{M}}{>}_{\phi ,\theta }$ is practically constant for r ∈ [2 : 8]r₀, and it is the divergence of this term (i.e of 𝒯_rφ in Eq.B.2) that contributes to the radial mass transport.

Therefore, the mass flux is completely due to the angular momentum extracted from the disc in the vertical direction (in this case it is the value of 𝒯_φ,θ in Eq.B.2 at z = ±H that contributes to τ_θ). This result is in agreement with the one described in the fiducial simulation of Béthune et al. (2017) (their Fig.7).

Appendix C Boundary conditions for MHD simulations

In our MHD simulations, we considered the following boundary conditions:

• Hydrodynamical quantities:

  1. Inner radial boundary: We extrapolated the density and the energy from the active cells into the ghost cells according to the initial power law profile. We extrapolated the azimuthal velocity component of the active cells into the ghost cells according to the Keplerian profile. The polar component of the velocity was copied from the active cell into the ghosts cells. We closed the inner boundary by setting to zero the radial component of the velocity.

     Fig. B.2 Computation of the radial and vertical contribution to the mass flux (Eq.B.2) compared to the measured mass flux (1.h.s. of Eq.B.2). Top panel: Temporal average. bottom panel: Spatial average.

  2. Outer radial boundary: We extrapolated the density and the energy from the active cells into the ghost cells according to the initial power law profile. The azimuthal component of the velocity was extrapolated according to the Keplerian profile. The polar component of the velocity was copied from the active cell into the ghosts cells. For the radial velocity component, we compared the value of the active cell to a threshold value ${\tilde{v}}_r$ proportional to a fraction of the sound speed. We set the velocity in the ghost cells equal to ${\tilde{v}}_r$ if the velocity of the active cell is smaller than the threshold value; otherwise, we simply copied the value of the active cell into the ghosts.

  3. Polar boundaries We extrapolated the density and the energy fields according to the hydrostatic vertical equilibrium defined in Appendix A. We extrapolated the azimuthal velocity according to rotational equilibrium (Appendix A). We copied the radial and the vertical components of the velocity from the active cells into the ghost cells.

• Magnetic field components and EMF

  1. Inner and outer radial boundary: We linearly extrapolated the values of the electromotive forces from the active cells into the ghost cells. We copied the azimuthal and the vertical components of the magnetic field from the active cells into the ghosts. We derived the radial component of the magnetic field by requiring that the divergence of the magnetic field in the ghost cells is zero.

  2. Polar boundaries We linearly extrapolated the values of the electromotive forces from the active cells into the ghost cells. We copied the azimuthal and the radial components of the magnetic field from the active cells into the ghosts. We derived the polar component of the magnetic field by requiring that the divergence of the magnetic field in the ghost cells is zero.

Finally, we imposed a density floor by requiring that the minimum value of the density is given by: $${\rho }_{min}=\frac{B_r^2+B_{\theta }^2+B_{\phi }^2}{4v_{kin}{\mu }_0}$$(C.1)

where v_kin is the Keplerian velocity at the inner radial boundary. In other words, we forced the Alfvén speed to remain lower than few times the Keplerian velocity at the inner boundary.

References

All Tables

All Figures

	Fig. 1 HDBae simulation. 2D distribution in the (r, φ) plane of the mid-plane density, ρ, normalised over the azimuthally averaged value < ρ > at t = 10 orbits when the planet has just reached a Jovian mass. For the purposes of comparison with Fig. 6 of Bae et al. (2016), the y axis is plotted using a logarithmic scale.
In the text

Fig. 2 Top panel: azimuthally averaged surface density profile for simulation HDBae. Bottom panel: azimuthally averaged η parameter. The planet is kept on a fixed circular orbit at r/rp = 1 and the ticks on the x axis (as well as the vertical black lines) indicate the values of η ~ 0 with a positive slope. The horizontal black line corresponds to η = 0. The damping region extends radially from rmin/r0 to the dot dashed line.

In the text

	Fig. 3 Perturbed surface density of the HDBae simulation at the same snapshots represented in Fig. 2. The dashed lines indicate the radial position of the pressure bumps, i.e. the values of r for which η = 0 in Fig. 2. The dot-dashed line indicates the boundary of the damping region.
In the text

	Fig. 4 Disc cooling time, density, and opacity in the (R, Z) ≡ (r sin(θ), r cos(θ)) plane of the disc in thermal equilibrium before the insertion of the planet.
In the text

	Fig. 5 Same as Fig. 4 for azimuthally averaged quantities for the 3D simulation HDradJ with an embedded Jupiter-mass planet. The snapshot corresponds to the end of the integration at 150 orbital periods at the planet location (5.2 au).
In the text

Fig. 6 Perturbed surface density of the HDradJ simulation at different snapshot indicated in orbital period at the planet location. The letters in the left and central panels indicate the primary (P), secondary (S), and tertiary (T) spiral arms. The dashed lines indicate the position the radial position of the pressure bumps, i.e. the values of r for which η = 0, with a positive slope, shown in Fig. 7.

In the text

	Fig. 7 Top panel: azimuthally averaged surface density profile for simulation HDradJ at the same snapshot represented in Fig. 6. Bottom panel: azimuthally averaged η values. The dashed line indicates the boundary of the damping region.
In the text

Fig. 8 Perturbed surface density at t = 150 orbits at the planet location (top panels) and at t = 500 (bottom panels) orbits. From left to right: simulation sets HDiso, HDβc with βc = 1, HDrad, and HDβc with βc = 100, HDadia. We observe multiple rings formation that depend on the EoS. Rings are persistent on longer integration times. For the radiative case, we report the last snapshot, i.e. the one at 150 orbits in the bottom panel.

In the text

	Fig. 9 Azimuthally averaged surface density profile (top) and the η values (bottom) for simulation sets HDiso, HDβc, and HDadia at t = 500 orbits at the planet’s location.
In the text

	Fig. 10 Top panel: density distribution in the (r, Z) plane obtained for our simulation MHDβ5e3 with ΛA 1 at t = 300T0. Bottom panel: computation of $\dot{M}(r,\theta )$ (Eq. (15)), averaged over the time interval [100 : 300]T0.
In the text

Fig. 11 Perturbed surface density of the MHD simulations. From left to right: increase in the strength of the magnetic field. In addition, for a given value of the magnetic field, we can also increase the Ambipolar coefficient value. For the lowest value of the magnetic field (left panel), we see the formation of a secondary ring at r = 2 au with our lower value of ambipolar diffusion (ΛA = 10); thus, we do not run additional simulations in this case. In the case of βm = 105 (second panel and third panels) we observe a secondary ring for ΛA = 10; however, the formation of a secondary well contrasted ring requires a higher ambipolar diffusion strength (ΛA = 1, third panel). Finally, for βm = 5 × 103, there is no secondary ring formed for ΛA = 10 (not shown) nor for ΛA = 1 (fourth panel). We observe a poorly contrasted secondary ring for ΛA = 0.1 (fifth panel). Finally, we observe a well contrasted secondary ring for ΛA = 0.01 (right-most panel).

In the text

Fig. 12 Azimuthally averaged surface density profile (top) and azimuthally averaged η values (bottom) for MHD simulations presented in Fig. 11. Here, a well-contrasted secondary ring is formed. The planet is kept on a fixed circular orbit at rp = 5.2 au and the ticks on the x axis (as well as the vertical black lines) indicate the values of η ~ 0 with positive slope. The horizontal black line corresponds to η = 0. The damping region extends radially from rmin to the dotted line.

In the text

	Fig. 13 Top panel: azimuthally averaged surface density profile for simulations Hrad30, HradS and HradJ at 150 orbital periods at the planet location. Bottom panel: azimuthally averaged η values. Colour markers at η = 0 with a positive slope identify the position of the pressure bumps on the top panel.
In the text

	Fig. B.1 Computation of the α values from Eq.B.1 averaged over time for the fiducial MHD simulation with β = 5 103 and ΛA = 1.
In the text

	Fig. B.2 Computation of the radial and vertical contribution to the mass flux (Eq.B.2) compared to the measured mass flux (1.h.s. of Eq.B.2). Top panel: Temporal average. bottom panel: Spatial average.
In the text