© ESO 2020

1 Introduction

In 1772 Joseph-Louis Lagrange identified five equilibrium points (L₁, L₂, L₃, L₄, and L₅) derived from the restricted three-body problem (Lagrange 1772), in which a particle of negligible mass orbits under the action of two massive bodies (e.g., a star–planet system). Two of these points, L₄ and L₅, lie in the orbit of the smaller body (planet), each one at the vertex of an equilateral triangle with the opposing joint base formed by the line of the two massive bodies. L₄ is located ≃+π∕3 rad at the leading position of the planet, while L₅ is located at the trailing co-orbital region at ≃−π∕3 rad (with respect to the planet).

The geometry of L₄ and L₅ implies thatthe ratio of their distances to the barycenter is equal to the ratio of the two masses. Consequently, the net gravitational force from the planet–star system is zero at these locations. Hence, L₄ and L₅ should be linearly stable under small perturbations. Gascheau (1843) determined that for sufficiently small ratios (< 1∕27) of the star–planet mass system, it should accumulate nonmassive objects called Trojans (see Brouwer & Clemence 1961; Szebehely 1967; Sosnitskii 1996 for more details on the stability). The first Jovian Trojan (588 Achilles) was discovered by Max Wolf in 1906 (Wolf 1907). Since then, thousands have been reported¹. More recently, the first Neptune Trojan asteroids 2001 QR322 and 2008 LC15 (Sheppard & Trujillo 2010) were discovered.

Several efforts have been made to understand the dynamics of Trojans and their origins. For instance, Morbidelli et al. (2005) suggest that Trojan asteroids could have been formed in distant regions to be later captured into co-orbital motion during the migration of the giant planet in the context of the Nice model, where the evolution of the planets is followed after the gas disk hasdissipated (Tsiganis et al. 2005).

The origin of Trojans could also be connected to an even earlier stage of the Solar System, where there is a gas-rich environment. Several hydrodynamical simulations indeed show that asymmetric overdensities in the gas are produced in co-rotation with the planet, favoring L₅ over L₄ (e.g., de Val-Borro et al. 2006; Lyra et al. 2009). Therefore, the accumulation of Trojan dust should be larger in L₅ compared to L₄. However, this is contrary to observations of Jovian Trojans, which show a ratio of the number of asteroids N(L₄)∕N(L₅) ≈ 1.8 for Trojans with diameters D > 2 km (Yoshida & Nakamura 2005).

One possible explanation for this discrepancy is that we now observe the end result of multiple physical processes (e.g., drag forces, collisions, grain growth) that have modified the Trojan population since formation of the Solar System (e.g., Milani et al. 2017). For instance, Di Sisto et al. (2019) analyzed the observed Trojan population through numerical simulations and concluded that the Trojan escape rate from L₅ in the lifetime of the Solar System is ~1.1 times greater than that from L₄, and is mainly produced by gravitational interactions with the other planets (from Venus to Neptune). Pirani et al. (2019) studied the consequences of planetary migration on the minor bodies of the Solar System through N-body simulations. They find that inward migration produces a more populated leading swarm (L₄) than the trailing one (L₅), in agreement with observations, while a nonmigrating planet results in symmetric swarms. However, the study of these latter authors is limited by a simplified treatment of the drag force that mimics the effect of the gas-phase of the protoplanetary disk. A complementary explanation could be the lack of systematic observations, or bias in them.

The bias problem concerns two critical aspects. The first is related to the limited amount of time for the observation, which translates to a limitation in the absolute magnitude H that can be reached. For instance, Lagerkvist et al. (2002) show that observing Trojans to a limit of H = 11 mag indicates that the L₅ swarm is 75% of L₄, while down to H = 13 mag shows that L₅ swarm is 76% of L₄. The other aspect concerns the covered area to observe a Lagrangian swarm. Unfitted estimations of the density area lead to inaccurate population estimations. For instance, Jewitt et al. (2000) estimate the L₄ Trojan population by analyzing an area between L₄ and L₃, where the distribution is more spread than between L₄ and Jupiter. Follow-up observations indicate that their results were overestimated by 40% (Lagerkvist et al. 2002). Different inclination distributions of Trojans and low albedo (e.g., Yoshida & Nakamura 2005) also contribute to bias these detections, notably in the case of small Jovian objects with D ≲ 1 m (e.g., Jedicke et al. 2002). Large correction factors are therefore required to overcome these observational biases (Karlsson 2010).

An interesting scenario is the possibility of a co-orbital planet located at a Lagrangian point. In that context, Lyra et al. (2009) studied the effect of high-pressure regions around the Lagrangian points L₄ and L₅ of a giant planet. They show that large bodies (m–km) accumulate in tadpole orbits, suggesting the formation of a Trojan Earth-mass planet in situ. This model requires a considerably massive and dense (self-gravitating) disk. In a closely related context, Cresswell & Nelson (2009) investigated the evolution of Trojan planets embedded in a gaseousdisk, from which they can grow to become gas giant planets. The 2D numerical simulations of these latter authors show than once a Trojan planet is placed to grow, the system is stable for about ~ 10⁹ yr.

Despite the strong assumptions needed to produce high-mass Trojan planets, active searches for these objects in extra-solar systems are currently taking place. For instance, Giuppone et al. (2012) analyzed the possible detection of exoplanets in co-orbital motion, namely the TROY project², an observational and theoretical effort to understand the evolution of planetary systems from the characterization of (still not detected) exo-Trojan planets (Lillo-Box et al. 2018a,b). This characterization raises several open questions, such as: are large Trojan asteroids formed in situ (by aggregates of micron-cm particles to km rocks) or captured, for example during migration? What could be the maximum mass reservoir of dust that can accumulate in Lagrangian points? What is the dynamical evolution of the early dust grains initially present when a protoplanetary object appears?

In this work, we study the early evolution of Trojan dust with the aim being to understand the current configuration of Trojans in the Solar System and the possibility of finding Trojan exoplanets. We model the dynamical evolution of the primordial dust present at the early stage of a gaseous protoplanetary disk, with a particular focus on the dust concentration around the Lagrangian points L₄ and L₅, and its stability on short timescales (10⁴ yr).

The paper is structured as follows: In Sect. 2, we provide a description of the numerical setup. In Sect. 3, we discuss the physical considerations regarding vortices and Lagrangian points. In Sect. 4, we present the results of the gas and dust simulations. A discussion is provided in Sect. 5, and our main conclusions are drawn in Sect. 6.

2 Numerical model

We follow the evolution of a dusty, gaseous, viscous self-gravitating protoplanetary disk with an embedded planet. For that purpose, we divide our simulations into two stages; the first computes the gas hydrodynamics alone by solving the Navier–Stokes equation and a nonstationary energy equation using a revised version of the 2D FARGO-ADSG code (Masset 2000; Baruteau & Masset 2008). In our version, we model a passively heated disk irradiated by the central star. The disk includes a radiative cooling mechanism assuming black-body radiation, with a nonstationary energy equation. For more details, see Montesinos et al. (2016). The second stage follows the dust dynamics computed on top of the gaseous stage, where the outputs of the hydro-simulation (first stage) are used as inputs for the dust code (second stage), which computes the dynamics of the Lagrangian particles. A complete description of the dust code and our methodology can be found in Cuello et al. (2019).

The evolution of the disk is followed for about 10 000 yr (or 460 planetary orbits). To check that our results correspond to a stationary regime, we also run a supplementary model running up to 1000 orbits for our fiducial model, obtaining almost the same results. We do not present such a model in this work. A description of the dust and gas stage is outlined in the following section.

2.1 First phase: gas dynamic setup

Our fiducial model includes a Jupiter-mass planet located at 7.8 au. The choice of the planet location was motivated by its generated gap, which may be scalable to some observations such as the cavity reported in HD 100546, which could be carved by a Jupiter-mass planet (e.g., Bouwman et al. 2003; Tatulli et al. 2011).

The orbital period of the planet is about T_p = 21.8 yr, and its gravitational effect is introduced smoothly by a taper function such that its final mass is reached after N_taper = 10 orbits; $$M_{\text{p}}(t)=\{\begin{array}{ll}{M}_{\text{p}}\cdot {\sin }^2\left(\frac{\pi }{2}\frac{t}{N_{\text{taper}}{T}_{\text{p}}}\right)\hfill & t<10T_{\text{p}}\hfill \\ M_{\text{p}}\hfill & \text{otherwise}.\hfill \end{array}$$(1)

The initial density profile of the disk is assumed to be $$\Sigma (r)=8.9\left(\frac{\text{au}}{r}\right)\text{gr} {\text{cm}}^{\text{-2}},$$(2)

distributed in a radially logarithmic grid of 512 (radial) × 1024 (azimuth) sectors with cylindrical coordinates r, ϕ. The inner radius r_in of the disk is located at 2.5 au from the central star, and the outer radius r_out at 15 au. From the density profile (2) and the grid limits, the initial disk mass gives $M_{\text{disk}}^{\text{gas}}={\displaystyle {\int }_{r_{\text{in}}}^{r_{\text{out}}}2}\pi \Sigma (r)r\text{d}r=2.4\times {10}^{-3}{M}_{\odot }$, which is compatible with typical circumstellar disk masses (e.g., Andrews & Williams 2005).

Despite the fact that this mass corresponds to a low disk mass, we include self-gravitating effects in our calculations. These effects could impact the dust dynamics, even for a large Toomre parameter such as that obtained from our disk parameters (Baruteau & Zhu 2016).

We solve a nonstationary energy equation, in which the source term is given by stellar irradiation $$Q_{\star }^{+}=(1-\beta )\frac{L_{\star }}{4\pi r^2}\frac{H}{r}\left(\frac{\text{d}\ln H}{\text{d}\ln r}-1\right),$$(3)

where β is the albedo (set to zero), L_⋆ the stellar luminosity, and H the scale-height computed assuming local hydrostatic equilibrium $H=\frac{c_{\text{s}}}{v_{\varphi }}r$, with c_s and v_ϕ being the sound speed and azimuthal velocity, respectively.

The radiative cooling mechanism is given by black-body emission, $$Q^{-}=\frac{2{\sigma }_{\text{SB}}{T}^4}{\tau },$$(4)

where σ_SB is the Stefan-Boltzmann constant, $\tau =\frac{1}{2}\kappa \Sigma$ the optical depth, and T the midplane temperature. We simplify our calculations by using a constant disk opacity κ = 1 cm² g⁻¹. This choice is justified by noting that the Rosseland mean opacity in the temperature and density range of this work gives a mean value κ ~ 1 cm² g⁻¹ (Semenov et al. 2003).

The disk temperature is initialized by $$T(t=0,r)=h^2\frac{G{M}_{\star }}{r}\frac{\mu }{\mathcal{R}},$$(5)

where h = h₀r^f is the aspect ratio of the disk. The factor f is the flare of the disk, which is assumed to be f = 1∕7, obtained from equating $Q_{\star }^{+}=Q^{-}$, corresponding to a quasi-steady regime. This choice helps to reach a quasi-steady situation as fast as possible. The initial aspect ratio is set to h₀ = 0.05. A complete description of the energy equation and its initialization can be found in Montesinos et al. (2016) and Cuello et al. (2019).

To explore other disk configurations we compute different models taking the same density distribution (Eq. (2)), but changing the mass of the planet M_p = {1, 5}M_J, the stellar luminosity L_⋆ = {1, 5} L_⊙, and the disk turbulent viscosity α = {10⁻⁴, 10⁻³, 10⁻²}³, where α is the turbulent viscosity prescription from Shakura & Sunyaev (1973). We also run an extra migrating model, where the planet was initially located at 5.2 au. We summarize the explored parameter space in Table 1.

2.2 Second phase: dust dynamics setup

The dust phase is a post-processing calculation where an ensemble of N-independent dust particles react to drag forces produced by the gas and gravitational forces from the star, the gaseous disk, and the embedded planet. The particles are described as Lagrangian test particles that do not interact with one another. The dust code is based on previous work by Paardekooper (2007); Zsom et al. (2011).

The drag force on each particle is computed as $$F_{\text{drag}}=-\frac{\Omega (r)}{\text{St}}\Delta v,$$(6)

where Ω(r) is the angular velocity of the gas, Δv the relative velocity between the gas and the dust, and St the Stokes number, which is computed as an interpolation between Epstein and Stokes regimes, which is valid for small and large particles, respectively (see Appendix C from Cuello et al. 2019, and Stammler 2017).

Although we use a more convenient way to compute the Stokes number, most of the particles in the simulation follow the Epstein regime, which is valid for particles with sizes smaller than 9/4 of their mean free path (Weidenschilling 1977). In the Epstein regime, the Stokes number for each particle i is given by $${\text{St}}_i=\frac{{\rho }_{\text{d}}{s}_i}{{\rho }_{\text{g}}{c}_{\text{s}}}\Omega (r),$$(7)

where ρ_d is the bulk density of the dust, s_i the size of the particle i, ρ_g the local gas density, and c_s the local sound speed. The Stokes number for each particle varies depending on the position r of that particle.

In our calculations, we also include turbulent diffusion of the dust, which is introduced as a random walk of the dust particles. At each time-step d t, particles are displaced in a random direction by a length $l_{\text{turb}}=\sqrt{(D_{\text{d}}\text{d}t)}$, where D_d = αc_sH∕(1 + St²) corresponds to the diffusivity coefficient for dust (Youdin & Lithwick 2007).

Our dust simulation uses 150 000 particles. The dust sizes s range from s_min = 1 μm to s_max = 1 cm, covering the whole dust-size spectrum logarithmically. We set an initial power-law size distribution n(s) ∝ s⁻¹, which gives a power-law dust-density distribution equal to the gas surface distribution. This choice implies approximately the same number of particles perbin size, which allows us to better follow the dust kinematics by directly comparing the assumed dust species. The dust content is introduced at the beginning of the simulation, and is assumed to be a mixture of pyroxene and ices, with an intrinsic bulk density of ρ_d = 2.0 g cm⁻³. The initial gas-to-dust ratio is fixed to 100:1.

To calculate a physical quantity such as the effective mass accumulated in a region, or for instance the dust continuum emission, we are required to specify the size distribution n(s) which can be treated as a free parameter, as well as the total dust mass. This free choice does not affect the initial dust size scaling n(s) ∝ s⁻¹, which was chosen for computational reasons. In other words, the simulations give us the spatial distribution of particles, but we need to use a realistic size distribution to compute realistic physical quantities.

Consequently, once the simulation is complete, we re-scale the size distribution by n(s) ∝ s^−p, using p = 3.5 (Dohnanyi 1969). Therefore, the dust mass per bin size M_i can be computed as $$M_i=\varepsilon M_{\text{gas}}\frac{s_i^{4-p}}{{\sum }_i{s}_i^{4-p}},$$(8)

where ɛ is the initial dust-to-gas mass ratio (i.e., 1:100), and M_gas the total mass gas. The surface density σ_i for each i-bin can then be computed as $${\sigma }_i=\frac{N_i(r,\varphi )}{A(r,\varphi )}\frac{M_i}{{\sum }_i{N}_i(r,\varphi )},$$(9)

where N_i(r, ϕ) is the number of dust particles per bin size, A(r, ϕ) the surface area of each cell grid, M_i the dust mass per bin size computed above.

3 Physical considerations

3.1 Vortices

The vortensity gives a measure of the local rotation of the fluid, and its minimum indicates the presence of a local pressure bump. Particles are attracted to higher pressure regions, even for a slight enhancement from the background (Klahr & Lin 2001). Therefore, dust accumulation is expected at the center of vortices (e.g., Barge & Sommeria 1995; Chavanis 2000; Pinilla et al. 2012; Baruteau et al. 2019).

Vortices are expected to appear in the disk as a consequence of Rossby wave instabilities RWI (Lovelace et al. 1999; Li et al. 2000)or baroclinic instabilities Klahr & Bodenheimer (2003). Inside them, pressure maxima are present, with a minimum of the gas vortensity. Typical regions where these instabilities appear, excited by density gradient, are at the edges of the gap carved by a planet (e.g., Lin 2012).

To visualize the vortices generated in the disk, it is useful to compute the amplitude of the vortensity perturbations relative to its initial profile, that is, (ω − ω₀)∕ω₀. The vortensity field is defined as $$\omega =\frac{\nabla \times v}{\rho },$$(10)

where v is the local gas velocity, and ρ = Σ∕H is the vertically averaged density.

As our simulations are in 2D, we simply use ω to refer to the z-component of ω. The quantity ω₀ corresponds to the value of the (z-component) vortensity for the initial disk profile.

3.2 Lagrangian points L₄ and L₅

We are interested in the dust accumulation in the vicinity of the Lagrangian points L₄ and L₅. To compute their accretion, we need to establish the criterion of dust trapping in the vicinity of those points. A simple approximation can be obtained by solving the three-body problem consisting of two massive bodies of mass M_⋆ (primary) and M_p (secondary),where the secondary is moving in circular orbit around their mutual center of mass, and a third particle (test particle of negligible mass) moves in the same plane. Since we are interested in the motion of the test particle, it is convenient to refer to the coordinate system of the particle. In that case, the system is rotating with an angular speed $\Omega =\sqrt{GM/a^3}$, where M = (M_⋆ + M_p), and a the separation between the two bodies. The dynamic behavior of the third body can be obtained therefore from the effective potential $$\Psi =-G\left(\frac{M_{\star }}{r_1}+\frac{M_{\text{p}}}{r_2}\right)-\frac{1}{2}{r}^2{\Omega }^2,$$(11)

where r₁ and r₂ are the distances between the test particle and M_⋆, and M_p respectively, and r is the distance to the center of mass of the two bodies (in practice the distance to M_⋆). Equation (11) is called the co-rotating effective potential.

In Fig. 1, we plot the co-rotating equipotential contours from Eq. (11) using the parameter of our simulation: M_⋆ = M_⊙ and M_p = 5 M_J, both separated by 7.8 au. The figure shows the Lagrangian points L₄ and L₅ at the planet orbital radius, corresponding to stable equilibrium positions with a zero gradient potential. L₄ and L₅ are located ~π∕3 radians in azimuth ahead of the planet and behind it, respectively.

Following the contour plot in Fig. 1, we consider that a particle will be in tadpole orbit around a Lagrangian point if it belongs to the cylindrical areas drawn in Fig. 1 that surround L₄ and L₅. Each cylindrical area is delimited by the position of the Lagrangian points at ± π∕3 adding (or subtracting) an angle of ± π∕8 rad in the azimuthal direction, and between an effective capture radius given by the range r_p ± R_Hill, that is, $$\begin{array}{l}\pi /3 (L_4)-\pi /8<\varphi <\pi /3 (L_4)+\pi /8,\\ -\pi /3 (L_5)-\pi /8<\varphi <-\pi /3 (L_5)+\pi /8,\\ r_{\text{p}}-R_{\text{Hill}}<r<r_{\text{p}}+R_{\text{Hill}},\end{array}$$

where r_p is the planet radial location and R_Hill is the Hill radius⁴. We show a posteriori that the libration of particles around a Lagrangian point occurs in such a pre-defined area.

At each time-step, we count particles that enter and leave the defined region. This method enables us to calculate the net dust flux in L₄ and L₅. At late evolutionary stages of the disk, when a stationary regime is reached, the Lagrangian points neither accumulate nor lose material.

Fig. 1 Equipotential lines according to Eq. (11) in a rotating frame with the planet, applied to a model with parameters M⋆ = 1 M⊙, rp = 7.8 au, and Mp = 5 MJ. The leading Lagrangian point (L4) is located at + π∕3, while the trailing point (L5) is at − π∕3. A stable librating azimuthal angle is assumed to be ±π∕8 around eachLagrangian point.

4 Results

4.1 Gas and dust evolution

Figure 2 shows the evolution of the disk model 1 for the timescales 46, 230, and 459 orbits (as one planetary orbit takes 21.8 yr, this corresponds to 1000, 5000, and 10 000 yr, respectively), the different subplots include: the gas density, dust density, particle distribution, and gas-to-dust ratio for all the sizes ranging from 10⁻⁴ to 10⁰ cm. All models start with a gas-to-dust ratio equal to 100:1.

During the evolution, we observe the formation of an inner disk and an outer one as the planetary gap develops, along with typical density wakes at the Lindblad resonances. By the end of the simulation, the planet practically depletes the co-rotating zone from gas (bottom of Fig. 2). As expected, the dust evolves differently depending on its coupling to the gas.

At 46 orbits (upper-panel, Fig. 2), a gap is already present in the disk. However, a horseshoe structure of dust with particles of all sizes populates the co-rotating zone. The dust density inside the gap is very low ~ 0.003 gr cm⁻², with a gas-to-dust ratio of ~528 : 1, except at the two Lagrangian points, where the gas-to-dust ratio is reduced to 13:1 (L₄) and 9:1 (L₅), indicating dust accumulation. In the outer disk, two blobs of dust appear, corresponding to two vortices (see Sect. 4.2). The vortices are not noticeable in the gas. The gas-to-dust ratio at the large swarm reach ~ 8 : 1 (Fig. 2: coordinates x ~−7;y ~−8 au), showing efficient dust trapping.

After 230 orbits (middle panel, Fig. 2), the gap is more depleted of gas, while the co-rotation zone is still populated by dust of all sizes. The mm–cm particles inside this region continue to accumulate around the Lagrangian points. The gas-to-dust ratio in these points is reduced to 4.4:1 (L₄) and 3.8:1 (L₅). The two dusty blobs present in the outer disk at 46 orbits have collided to form a single one (located at coordinates x ~ 10; y ~ −5 au), with a total mass of 16.3 M_moon, and a gas-to-dust ratio of about 13.8:1.

At the end of the simulation (orbit 459, bottom panel, Fig. 2), μm and submm particles are librating in horseshoe orbits around L₄ and L₅, while mm–cm particles are confined in tadpole orbits around either L₄ or L₅. From Fig. 2 (orbit 459), the leading swarm (L₄) accumulate more small particles than the trailing Lagrangian point (L₅). However, there is more mass accumulated in L₅ than in L₄ since most of the mass is concentrated in large particles (mm–cm), which are concentrated mostly around L₅.

Fig. 2 Evolution of the gas and dust in model (1), for three evolutionary stages; 46, 230, and 459 orbits. The dust panels show the full size range of the particles, from 10−4 to 100 cm. The panels have been rotated to always keep the planet at the same location x = 7.8 au, y = 0.

Fig. 3 Top: Evolution of the vortensity and dust in model 1, for three evolutionary stages: 46, 230, and 459 orbits. Bottom: Dust distribution of particles with sizes in the range 10−2 − 100 cm. The planet rotates in the counter-clockwise direction. The panels have been rotated to always keep the planet at the same location, namely x = 7.8 au, y = 0.

4.2 Vortices

Figure 3 shows the amplitude of the vortensity deviations from the background for model 1, that is, (ω − ω₀)∕ω₀, where Eq. (10) was used to compute ω for the outputs: 46, 230, 460 orbits. For comparison, we plot the dust distribution by sizes next to the vortensity. For a better analysis, we divide the disk into three regions: (a) the inner disk, (b) the gap, and (c) the outer disk, as described below.

• (a)

  In early stages, the inner disk shows a homogeneous vortensity, with a local minimum (blue color) developing at the edge of the inner disk just next to the gap. The bottom panel of Fig. 3 shows the size distribution of dust, where sub-millimeter to cm particles (yellow color) are observed. As the disk evolves, two dusty rings concentrating submm to cm grains develop. The location of such rings coincides with the location of the minimum in vortensity.

• (b)

  The gap exhibits a local minimum vortensity with a typical horseshoe structure, trapping submm and cm sizeparticles. At early stages, two prominent vortices are identified at the Lagrangian points L₄ and L₅, where L₅ covers a larger area than L₄ (two blue islands observed in the gap). Hence, the Lagrangian points start to collect large particles (mm–cm). The vortices persist during the simulation, trapping more dust particles as the disk evolves. Trojan dust of 10⁻²−10⁰ cm lies in L₄ and L₅ at the end of the simulation.

• (c)

  Initially, the external disk exhibits two notorious vortices (blue islands), which start to trap dust. With time, the two vortices collide to form one single vortex, resulting in a larger banana-shaped vortex, concentrating submm to cm particles. The dust mass of the final outer lobe or satellite reaches ~ 16.6M_moon with a gas-to-dust ratio of about 9:1. One should note that the dust distribution in Fig. 3 only shows submm to cm grains because these sizes are effectively trapped in the local pressure maximum and vortices.

At the end of the simulation, mm to cm particles have been efficiently trapped at the Lagrangian points by small vortices created by the planet–star system. The vortices continue to be present at the end of the simulation although the size of the now-single vortex has reduced over time. If there were more material within the gap, the slightly higher pressure at the center of the vortices would continue to trap dust. Besides the populated Lagrangian points, two concentric rings at the inner region develop. These are dust traps created by planetary wakes.

4.3 Dust accumulation by sizes around Lagrangian points

4.3.1 Gas–dust interaction

Figure 4 shows the evolution for each dust bin size (μm to cm) at different evolutionary stages: 300, 2 000, 4 000, 6 000, and 10 000 yr. We observe that the largest particles in the range of0.3–1 cm (yellow) quickly accumulate at the Lagrangian points. In less than 1000 yr (45 orbits), L₄ and L₅ are already settled, including a forming system of two inner rings (at ~4 au), originated from planetary wakes. On the same timescale, small particles in the range of 10^−1.5−10⁻¹ cm (green color) are also trapped at the Lagrangian points. At the end of the simulation, a relatively massive swarm of dust of ~ 16.6 M_moon is formed bymm–cm particles (yellow/green) located at the outer region showing a banana-shaped dust concentration. The swarm corresponds to the leftover of the extinct external dusty disk.

Particle sizes in the range of 10^−2.5–10^−1.5 cm (green–blue) behave in a slightly different fashion. A horseshoe starts to develop at 1000 yr. As the disk evolves, a gap is carved by the planet. Soon, the gap reduces its number of green–blue particles, which remain mostly trapped at the Lagrangian points. An external disk is also created as the dust evolves. The external swarm of 16.6 M_moon described above is hidden in this outer ring of 10^−2.5–10^−1.5 cm particles. Micrometer to 10^−2.5 cm (blue–violet) particles are distributed all over the disk; in the inner disk, within the gap (librating in horseshoe orbits around L₄ and L₅), and in the outer disk.

Once the simulation ends, we analyze the distribution of particle size and Stokes number around the Lagrangian points, computing how much material is finally concentrated in L₄ and L₅. In Figs. 5 and 6, we plot the mass spectrum of L₄ and L₅ for model 1, and their ratio L₄∕L₅ after 10 000 yr (460 orbits) of evolution as a function of particle size and Stokes number, respectively.

In Fig. 7, we plot the accumulated effective mass per radial bin. We notice that centering an annulus at the planet location r_p, the capture area extends from r_p + 2R_Hill towards the outer region, and r_p − 1.5R_Hill towards the star. However, the effective capture radius, defining tadpole orbits, is mostly inside r_p ± R_Hill as defined by Eq. (14).

We observe the following interesting features from the mass spectrum: (i) There is an asymmetry in the total mass accumulated in each Lagrangian point. L₄ accumulates 2.1 M_moon, while L₅ accumulates 3.0 M_moon (Figs. 5 and 7). (ii) Most of the effective mass trapped in L₄ and L₅ is found in particles in the range of ~0.03 to 1 cm (top panel of Fig. 5). A small amount of micron-sized particles are also trapped in the Lagrangian points. (iii) The mass asymmetry between L₄ and L₅ depends on the particle size. Small particles in the range 10⁻⁴ to 10^−1.5 cm are more abundant in L₄, while particle of 10^−1.5 to 1 cm are more abundant in L₅ (see bottom panel of Fig. 5). The transition (when mass in L₄ is equal to that in L₅) happens for 10⁻² cm particles. However, the effective mass (all sizes) is always larger in L₅ than L₄. (iv) The range of the Stokes number of particles trapped in the Lagrangian points lies in the range St ~ 10⁻⁶−10. Particles trapped in L₄ (small particles) are dominated by Stokes number St < 0.1. Particles with large Stokes numbers, that is with St > 0.1, are commonly trapped in L₅.

Fig. 4 Evolution in time of the dust for model 1 as a function of particle size, ranging from μm to cm. The panels have been rotated to always keep the planet at the same location x = 7.8 au, y = 0.

4.3.2 Mass accretion at Lagrangian points

The dust accumulation rate around a Lagrangian point differs depending on the specific particle size. Figure 8 shows mass as a function of time in the vicinity of Lagrangian points for a specific range of sizes for model 1, that is: 10⁻⁴−10⁰, 10⁻⁴−10⁻³, 10⁻³−10⁻², 10⁻²−10⁻¹, and 10⁻¹−10⁰ cm. Initially, the total dust mass grows exponentially. During the first 50 planetary orbits, L₅ accumulates ~ 2.7 M_moon of dust, while L₄ reaches about ~ 1.9M_moon, where the main contribution of the mass comes from the largest particle size, namely, 10⁻¹−10⁰ cm. After that, the accretion rate is almost completely halted, reaching its final mass of 2.13 M_moon for L₄, and 3.0 M_moon for L₅.

Lagrangian points are almost completely depleted of micron-size particles. Figure 8 shows that particles in the range of 10⁻⁴−10⁻³ cm do not contribute to the mass trapped in L₄ and L₅. The generated vortex at the Lagrangian points is not strong enough to attach stable orbits for this size range. However, a considerable amount of gas is located around L₄ and L₅, namely 4.8 M_moon of gas in L₄, and 5.6 M_moon in L₅ (with mean gas density of ~0.4 and ~0.5 g cm⁻² respectively), which is responsible for the pressure bump trapping mm–cm dust. At the end of the simulation, the gas-to-dust ratio is 2.2:1 for L₄ and 1.9:1 for L₅.

Table 2 summarizes the dust properties at the Lagrangian points at the end of each simulation, for all the models in Table 1. $M_{{\text{L}}_4,{\text{L}}_5}$ represents the accumulated dust mass, ${\sigma }_{{\text{L}}_4,{\text{L}}_5}^{\text{dust}}$ the dust density, Σ_gas∕σ_dust the gas-to-dust ratio, and ⟨T⟩ the mean dust temperature, where we assume that the dust has the same temperature as the gas. All these quantities are evaluated at each Lagrangian point L₄ and L₅, inside the effective capture region given by Eqs. (12)–(14).

Fig. 5 Mass spectrum of model 1 as a function of the size of the dust accumulated in L4 and L5 for thelast evolutionary stage (460 orbits). The bottom panel shows the mass ratio L4 ∕L5 as a function of dust particle size.

4.4 Instabilities at the Lagrangian points

Our fiducial model (model 1) shows a relatively stable behavior of mass accumulation at the Lagrangian points (Fig. 8). We explore the impact of turbulent viscosity α, planetary mass M_p, and stellar irradiation L_⋆ on the stability of trapped dust. We present our findings below.

Fig. 6 Mass spectrum of model 1 as a function of the Stokes number for the dust accumulated in L4 and L5 at the last evolutionary stage (460 orbits). The bottom panel shows the mass ratio L4∕L5.

Fig. 7 Histogram of the effective dust mass accumulated at the end of the simulation inside L4 and L5, distributed by radial bins. This defines a capture region given by rp + 2RHill from the planet towards the outer region, and rp − 1.5RHill from the planet towards the star. An effective capture radius is expected at rp ± Rhill.

Fig. 8 Dust mass accumulated in time around the Lagrangian points L4 (continuous line) and L5 (dotted line) for model 1. The evolution is shown for different sizes, including the effective (all sizes) accumulated mass (black continuous line).

4.4.1 Effect of turbulent viscosity

The situation is somewhat different if the turbulent viscosity is increased. Models 1, 2, and 6 share the same parameters, apart from the viscosity: α = 10⁻⁴ (model 1), α = 10⁻² (model 2), and α = 10⁻³ (model 6) (see Table 1).

In Fig. A.1, we plot the dust-mass evolution around L₄ (solid line) and L₅ (dotted line) for model 2. In this enhanced viscosity model (α = 10⁻²), the effective mass reaches final values of 0.6 (L₄) and 2.2 (L₅) M_moon. This is in stark contrast to the accumulated mass in the low viscous model 1 (α = 10⁻⁴), where the final mass reaches 2.1 (L₄) and 3.0 (L₅) M_moon.

Another remarkable difference from model 1 is that the only particles that effectively accumulate in both Lagrangian points are in the range of 0.3–1 cm. Small particles in the range of 10⁻²−10⁻¹ cm accumulate in L₄ and L₅ for a short period only. An instability is triggered after 200 orbits, evacuating these particles from the Lagrangian points (Fig. A.1). When the gap is being carved, small particles initially located in co-rotating orbits are dragged towards the inner disk by gas accretion. They never follow stable tadpole orbits around the Lagrangian points for such a high-viscosity model. Furthermore, increasing the gas viscosity increases the transport of angular momentum, promoting gas accretion towards the star. The enhanced accretion produces a less active vortensity. The less active vortex does not retain the smaller (10⁻⁴−10⁻² cm) particles efficiently.

There is also the effect of diffusion of dust particles. In our simulations, such diffusion is modeled through the coefficient D_d = αc_sH∕(1 + St²), which is proportional to the turbulent α-viscosity. The enhanced viscosity promotes radial diffusion of particles, helping to trigger the observed instabilities, and making the accumulation of dust particles more difficult in the Lagrangian points. However, despite having the effect of higher turbulent viscosity and a diffusive mechanism for particles, we still have a significant accumulation of mm–cm particles in L₅.

Figure A.5 (top panel) shows the mass spectrum of the dust accumulated in the Lagrangian points for model 2. The bottom panel shows the mass ratio L₄∕L₅ per bin size. The total mass ratio gives 0.3 (compared to 0.7 for model 1). This comparison suggests that a higher turbulence in the disk stimulates the evacuation of dust grains from L₄ rather than L₅ (therefore the mass ratio L₄∕L₅ is reduced for high-viscosity models). In other words, the lower the viscosity, the greater the similarity between L₄ and L₅ in terms of dust mass.

With an intermediate value for the turbulent viscosity (α = 10⁻³, model 6), we find an intermediate situation between model 1 (α = 10⁻⁴) and model 2(α = 10⁻²). See Fig. A.2 for the dust accumulation rate, and Fig. A.6 for the final mass spectrum of model 6.

4.4.2 Effect of planetary mass

In model 3, the mass of the planet is increased to 5 M_J, keeping the same parameter space as in model 1. The massive planet of model 3 carves a broader gap by a factor $M_{\text{p}}^{1/2}=5^{1/2}$ (Kanagawa et al. 2016), and on a shorter timescale compared to model 1. Gravitational perturbations from the planet severely affect the evolution of dust particles around the Lagrangian points.

Figure A.3 shows the dust evolution around L₄ (solid lines) and L₅ (dotted lines) for model 3. In this case, only L₅ accumulates grains, where most of them are large particles in the range 10⁻¹−10⁰ cm.

The rapid formation of the gap for this model leads to an enhanced evacuation of gas inside the co-rotation zone, through a region slightly closer to the leading zone of the planet, thereby removing most of the particles attached to L₄. The larger the planetary mass, the less efficient the dust capture at a Lagrangian point.

4.4.3 Effect of stellar irradiation

Comparing models 3 and 4 reveals some aspects of the influence of stellar irradiation. Model 3 has a central star with 5 L_⊙, while model 4 is a colder model with 1 L_⊙. In Model 3 (described in Sect. 4.4.2, particles at L₄ are completely decoupled after 200 orbits of evolution. For a colder model such as model 4, particles that once accumulated around L₄ also become unbound as a consequence of the huge gap created by the massive 5 M_J planet), but at a later evolutionary stage, after 340 orbits (compared to 200 orbits for the hotter model 3). The broader gap leads to particles decoupling from L₄ anyway. In a colder disk, the pressure scale-height is smaller, making it easier for a planet to open a gap and produce variations in the pressure field with sharper transitions in the radial direction, which translates to the promotion of RWI and hence the formation of vortices that last longer.

4.5 Evolution and trajectory of dust particles

For better visualization, we divide the last evolutionary stage of model 1 into three regions: (a) inner ring, (b) the gap, and (c) the outer disk. Each zone features different morphologies worth studying (see Fig. 4). We randomly select a number of particles of different sizes from the last step of the simulation belonging to specific zones of the above three regions and trace back along their evolutionary path, starting from the initial time-step to the last one.

4.5.1 Inner ring

At the inner rim of the disk, two concentric mm–cm dusty rings produced by planetary wakes are located at r ~ 3.8, and r ~ 5.6 au, revolving at resonances of 3:1 (3.8 au) and 5:3 (5.6 au), respectively, with respect to the planet. The rings are clearly visible at the end of the simulation (top right panel of Fig. 4).

We selected particles belonging to the second dusty ring at r = 5.6 au (for the output 459 orbits, model 1) and traced back their trajectories. Figure 9 shows the evolution of their radial position as a function of time, starting at t = 0 to t = 10 000 yr. The color bar indicates the size of each particle. The two horizontal red lines indicate the position of the Lagrangian libration zone (the planet is located at r_p = 7.8 au). The particles that end up trapped at the inner ring (r = 5.6 au) came from different regions of the disk, initially located between the radii at 5.5 au and 8.8 au (see t = 0 in Fig. 9).

The mm–cm particles beyond 6.5 au migrate due to drag forces to their final position at 5.6 au (Fig. 9), while μm particles of this ring were initially located at the same original radial distance (close to 5.6 au), oscillating with an amplitude of ~ 1 au around r ~ 5.6 au while they travel through the full orbital path around the star. A few μm particles were initially found at the co-rotation zone of the planet; being coupled to the gas, these particles were relocated to the inner ring by gas accretion. It is interesting to note that the ring located at r ~ 5.6 au is rather composed of two close small rings; one accumulates submm particles, while the other accumulates mm–cm particles.

Fig. 9 Radial trajectory of randomly selected particles (μm to cm) that end up trapped in the second inner ring located at r = 5.6 au featured at the end of the simulation in model 1.

4.5.2 The gap: tadpole and horseshoe orbits

We are interested in the libration of particles around L₄ and L₅, including tadpole and horseshoe trajectories between both locations. We selected several particles from the last output in the size range 10⁻⁴−10⁰ cm that end up trapped inside the librationarea of L₄ defined by the azimuthal position: π∕3 − π∕8 < ϕ < π∕3 + π∕8, and radius: r_p − R_Hill < r < r_p + R_Hill (see Sect. 3.2, Eq. (12)). Figure 10 shows the evolution in time of the radial and azimuthal coordinates of each particle. The horizontal red lines in Fig. 10 indicate the defined libration area.

As soonas these particles start their motion, they engage into a damped oscillation mode between the capture region r_p ± R_hill ~ (R_hill ~ 0.5 au for model 1). The larger the particle size, the larger the damping. Millimetric and centimetric particles become immediately constrained to oscillate very close around r_p (under-damped regime). Submm and μm particles also oscillate in an under-dumping mode but with a larger amplitude around r_p (see top panel of Fig. 10). Regarding the azimuthal coordinate, mm and cm particles oscillate with small amplitude around a full orbit around the star; they quickly reach the equilibrium position at L₄ (recall that we selected only particles ending up in L₄; see particlesinside the horizontal red lines at the bottom panel of Fig. 10). However, small particles (μm and submm) librate with larger amplitudes. Some of them become trapped in the other Lagrangian point L₅ for a while, going back and forth from L₄ to L₅, before finally ending up inL₄.

In summary, mm–cm particles that end up in L₄ oscillate around L₄ (tadpole orbits) during the whole simulation with an amplitude defined by the libration amplitude (Eq. (12)). Submm and μm particles oscillate around both Lagrangian points (horseshoe orbits). L₄ catches more smaller particles (μm–submm), and L₅ efficiently traps larger particles (mm–cm).

A major result is that all particles that end up at L₄ (or at L₅), regardless of their size, were initially located inside the libration region defined by the effective capture radius (Eq. (14)), that is, an annulus centered at r_p with a width of Δr = 2 × R_Hill (~ 1 au for model 1). No particles initially located outside this area end up trapped around a Lagrangian point. Particles traveling from outer regions dragged by the gas pass through the gap but without being trapped in tadpole orbits around a Lagrangian point. The capture region suggests that the mass reservoir to be accumulated in a Lagrangian point should be $M_{\text{res}}^{L_4,L_5}={\displaystyle {\int }_{r_{\text{p}}-\Delta r}^{r_{\text{p}}+\Delta r}{\sigma }_{\text{dust}}}(r)2\pi r\text{d}r$, where σ_dust(r) is the initial dust density profile, and r_p the planet position.

This maximum imposes some constraint on the in situ formation scenario of Trojan planets. For instance, from the initial condition defined for model 1 (Sect. 2.1, Eq. (2)), we have σ_dust(r) = Σ_gas(r)∕100, which gives $M_{\text{res}}^{L_4,L_5}~184$ M_moon (2.3 Earth masses). For this model, the final masses around L₄ and L₅ are 2.13 and 3.0 M_moon, respectively (see Table 2), representing a capture efficiency of about 1.2% (L₄) and 1.6% (L₅) of $M_{\text{res}}^{L_4,L_5}$. Depending on the simulation parameters, the capture efficiency may vary. However, the mass reservoir is an independent constraint, which will not change when varying the parameters or if a larger or smaller disk is assumed.

Fig. 10 Trajectories of some particles (μm to cm) from model 1 that end up trapped inside the Lagrangian point L4 only. The top horizontal red lines represent the location L5, while the bottom horizontal lines represent L4. Some submm particles oscillate in horseshoe orbits (between both L4 and L5) before ending in L4. Large mm–cm particles oscillate only in tadpole orbits around L4.

Fig. 11 Radial trajectory of randomly selected particles (μm to cm) that end up trapped in the outer ring located at r = 11 au featured at the end of the simulation model 1.

4.5.3 Outer disk

The outer disk shows a wide ring structure located at r ~ 11 au composed of μm and submmparticles. Within this ring, a large vortex is produced as a consequence of the gravity of the planet, with an over-density of mm–cm particles with a total mass of 16.6 M_moon (bottom-right panel of Fig. 3).

We selected some of the trapped grains in the vortex to trace their trajectories back in time. Figure 11 shows the evolution of those particles starting from t = 0 to t = 10 000 yr. We find that the particles were initially located at radial positions between 8.8 to 15 au (the vortex is located at r ~ 11 au).

Some large (mm–cm) particles migrate outward, driven by forces from planetary wakes, while others migrate inward, influenced by drag. Submm dust particles trapped in the vortex were initially located only in orbits close to its final position inside the vortex, that is, r~ 11.5 au (see Fig. 11).

As shown in Fig. 11, mm–cm particles barely oscillate in the radial direction while they travel through the disk. In contrast, small, submm particles oscillate with decreasing amplitude in the radial direction while they move toward the vortex.

Analyzing the vortex dynamics as a whole, we find that at t = 1000 yr, it revolves with a mean motion resonance (MMR) of 3:2 with respect to the planet, shifting radially to its final position at r ~ 11 au reaching an MMR of 5:3 at the end of the simulation. The resonant vortex, with its dusty banana-shaped structure (bottom panel of Fig. 3), could feature some observational signatures such an asymmetry in the dust continuum of the dusty ring to which it belongs, peaking in scattered infrared light and submm emission (e.g., Bae et al. 2016; Baruteau et al. 2019). However, further radiative analysis is beyond the scope of this work.

In our simulations, the planet was not able to migrate. However, for completeness, we run two dedicated models: (i) one in which a 10 M_J planet ‘feels’ the disk, starting a type II migration regime; (ii) and an identical model, but without migration (see parameter Table 1, models 7 and 8). It is worth mentioning that, in general a massive planet migrates faster than a lower mass planet. For instance, at first order, a reference migration timescale would be ${\tau }_{\text{mig}}\propto M_{\text{p}}^{-1}$ (e.g., Armitage 2010). In our simulation, a 10 M_J planet starts a fast inward type II migration regime from 5.2 to 3.7 au in 460 orbits.

The effect of planet migration slightly reduces the total mass accumulated in both L₄ and L₅, but the asymmetry favoring L₅ over L₄ remains. Migration does not enhance material accumulation in L₄ over L₅, or trigger any destabilization mechanism around L₅ as proposed by Gomes (1998). The only difference we find is on the trajectories of micron particles coupled to the gas. Those particles are not captured in stable tadpole orbits in the horseshoe region, which is probably because of the large planetary mass used (needed to obtain a fast type II migration regime).

Fig. 12 Spectral energy distribution from the disk, and the Lagrangian point L5 from model 1 computed at the last evolutionary stage. Emission from the Lagrangian point comes from an optically thin dusty swarmof ~3 Mmoon with a mean temperature of ~55 K.

4.6 Radio flux from Lagrangian points

From the last outputs of the simulation, noting that the system has reached a quasi-stationary regime, we can estimate the emission from the disk and the Lagrangian points. Assuming that the source is located D = 150 pc away, the total disk flux can be estimated by integrating the Planck function B_ν over the disk surface; $F_{\nu }=(2\pi /D^2){\displaystyle {\int }_{{\text{R}}_{\text{min}}}^{{\text{R}}_{\max }}{B}_{\nu }}(T(r))r\text{d}r$, obtaining a peak flux of about ~450 mJy at 50 μm for model 1. On the other hand, the dust accumulated in a Lagrangian point contributes with a specific localized emission. If it comes from an optically thin region, the emission can be computed from F_ν = (1∕D²)M_dustκ_νB_ν(T), where M_dust corresponds to the dust mass located at the Lagrangian point, and T its temperature (Hildebrand 1983). Using the values from Table 1 (M_dust ~ 3M_moon, T ~55K), we obtain for model 1 an integrated flux of ~20 mJy also peakingat 50 μm. This is an idealized estimation, without taking noise or instrument limitations into account. The spectral energy distribution for this model is shown in Fig. 12.

5 Discussion

We studied the evolution of dust present in a viscous disk with an initial gas-to-dust ratio of 100:1. The disk has an embedded approximately Jupiter-mass planet located at 7.8 au which we follow for 460 orbits. The dust is treated as Lagrangian particles with a full spectrum of sizes ranging from 10⁻⁴ to 10⁰ cm. Our simulations are done over two stages: The first one computes the gas dynamics by solving the Navier–Stokes equations, including a nonstationary energy equation for an irradiated disk. The second stage computes the dust dynamics in which dust particles ‘feel’ the gravity from both the star and the planet, and the drag caused by the gaseous disk.

We mainly focus on dust dynamics around Lagrangian points L₄ and L₅, examining the impact of three parameters that play an important role in the evolution of the dust: the mass of the planet M_p, the turbulent viscosity of the gas α, and the stellar irradiation (to heat the gas) from the star L_⋆. Some general conclusions, independent of the parameter choice, arise from these models: Once the planet has carved a gap in the disk, two vortices appear located at the Lagrangian points L₄ and L₅, revealed by the vortensity minima (Fig. 3). These minima act as dust traps (Crnkovic-Rubsamen et al. 2015; Surville et al. 2016), accumulating material in the librating region around L₄ and L₅. In all the models, we find that L₄ invariably captures a larger number of small (μm-submm) dust grains than to L₅. On the other hand, mm and cm particles are always more abundant in the trailing L₅ Lagrangian point. Larger particles account for more mass; therefore, we always find that at the end of the simulation, L₅ has accumulated more mass than L₄. Typical values of the mass accumulated around L₅ are on the order of 2–3 M_moon.

The observed asymmetry applies to all the models, consistent with past theoretical works suggesting that in the presence of drag, the orbits of Trojans around L₅ are more stable than those around L₄ (e.g. Peale 1993; Murray 1994). The origin of this asymmetry is due to an over-density created at the trailing region of the planet (e.g., de Val-Borro et al. 2006; Lyra et al. 2009). Such an over-density (equivalent to a pressure bump) in L₅ can be understood as follows: the orbital motion of the planet induces a larger gas depletion in the leading direction because the planet excites pressure waves that can remove angular momentum away from Lindblad resonances (Goldreich & Tremaine 1979; Artymowicz 1993), thus producing an asymmetry between the density lobes at the leading and trailing zones. This explains, for instance, why in a pressureless or inviscid (no viscosity) model the Lagrangian points accumulate an equal amount of material. In our case, we show that when lowering the viscosity, the mass ratio L₄ ∕L₅ approaches unity.

A major result concerns the mass reservoir that feeds the Lagrangian points. According to our models, particles ending up trapped either in L₄ or L₅ at the end of the simulation were initially located only at the co-orbital region of the planet r_p within the range ~±2R_Hill (Fig. 7). No particles in further orbits, located in the external regions of the disk (i.e., r > r_p + 2R_Hill), were trapped in stable orbits by the Lagrangian points (Fig. 10). This suggests that the total mass available to be trapped in a Lagrangian point of a planet located at r_p, is present at the beginning of the disk evolution inside the co-orbital area defined by an effective capture range given by Eq. (14). Therefore, the mass reservoir to feed a Lagrangian point will be: $M_{\text{res}}^{L_4,L_5}={\displaystyle {\int }_{r_{\text{p}}-R_{\text{Hill}}}^{r_{\text{p}}+R_{\text{Hill}}}{\sigma }_{\text{dust}}}(r)2\pi r\text{d}r$, where σ_dust(r) is the initial dust density distribution. This constraint holds regardless of the free parameters or the radial extension of the disk; and its evaluation is only a function of the initial dust density distribution and the mass of the planet (to compute R_Hill).

However, not all this mass is effectively captured by a Lagrangian point. Roughly speaking, ~ 1–2% (1–3 M_moon) of the initial co-rotating dust M_res (~ 2.3 Earth mass)will end up trapped as Trojan dust. Interestingly, this estimated mass can be taken to be, for instance, the origin of a swarm of material that will produce a Trojan moon-like planet around a Jupiter planet, as suggested by Beaugé et al. (2007).

The lack of initial material in the co-orbital zone and the dust capturing percentage (~ 1–2%) in L₄ or L₅ makes it challenging to assemble a co-orbiting earth-size planet in a Lagrangian point from a primordial configuration of the protoplanetary disk (e.g., assuming a minimum mass solar nebula model Weidenschilling 1977; Desch 2007) unless a very thick and massive initial dusty disk with larger grains is considered. For instance, based on the parameters of model 1 (one Jupiter mass planet located at 7.8 au) and assuming a dust capturing efficiency of 1%, in order to assemble one Earth mass in a Lagrangian point we need an initial dust density value of about ${\sigma }_L^{\text{dust}}=54\text{g}{\text{cm}}^{-2}$ at that Lagrangian point. This value is consistent with the disk density model presented by Lyra et al. (2009).

Recent observations of Jovian Trojans in the Solar System indicate a large number of objects⁵ located in L₄ than in L₅ (Yoshida & Nakamura 2005; Nakamura & Yoshida 2008; Pitjeva & Pitjev 2020). In contrast, our results show that the trailing (L₅) point is more efficient at accumulating large grains.

We suggest that in the early evolution of a planetary system, in the pre-transitional phase, Trojans form in situ on short timescales (~ 10 000 yr) as a consequence of a massive planet. Hence, the chemical composition of Trojans located at L₄ or L₅ should be similar. Also, far away asteroids, such as those from the main belt in our Solar System, are probably built from different shapes and components than Trojans. Spectroscopic and photometric observations seem to suggest this, as most Jovian Trojans are D- or P-type, while those from the main belt are C- and S-type (Hartmann et al. 1987; Fitzsimmons et al. 1994). However, as the system evolved, they were probably contaminated by collisions, heavy bombardment, and modified by gravitational interactions with other planets (e.g., Pál & Süli 2004; Freistetter 2006), as suggested by the Nice model (Tsiganis et al. 2005).

We also propose that if a protoplanet is found in a transitional disk, for example inside the gap of the circumstellar disk around HD100546, theorized to be the consequence of a Jupiter-mass planet, similar to our model parameters (Bouwman et al. 2003; Tatulli et al. 2011), two asymmetrical swarms of dust should be located at L₄ and L₅. The detection of such swarms could be used to reinforce or infer the presence of an embedded putative planet.

The planetary system around PDS 70 constitutes another interesting astronomical laboratory because recent near-infrared SPHERE and NACO observations revealed the presence of a couple of planets within the disk cavity: PDS 70b (Keppler et al. 2018) and PDS 70c (Isella et al. 2019). The latter was discovered through Hα emission with MagAO and MUSE. Our model suggests that such planetary companions, particularly PDS 70b, should have gathered potentially observable dusty swarms in their Lagrangian points. However, a more detailed study of the dynamical interaction between PDS 70b and PDS 70c is required in order to give an accurate prediction.

In our in situ scenario, intrinsic parameters (M_p, α, L_⋆) play an important role in the primordial structure of Trojans, especially concerning their reported mass asymmetry. We briefly summarize these effects below.

5.1 Viscosity

Increasing the turbulent α viscosity of the disk enhances the accretion rate (by promoting angular momentum transport), making the dust trapping around the Lagrangian points more difficult. This result can be understood as follows: a high accretion rate means higher values of the radial velocity, producing a reduced potential vortensity at the Lagrangian points to act as attractors. Comparing models 1, 2, and 6 (Table 1) reflects the mentioned influence of the viscosity (Fig. 8 next to Fig. A.1). The lower the α, the more similar the mass in L₄ is to that in L₅.

5.2 Planetary mass

Increasing the mass of the planet does not help to increase dust accumulation around a Lagrangian point. For instance, Models 3 and 4 (Fig. A.3) possess a larger planetary mass (5 M_J), which more easily destroys the leading swarm in L₄ by removingits angular momentum through pressure waves from the massive planet (Goldreich & Tremaine 1979; Artymowicz 1993).

5.3 Stellar irradiation

Stellar irradiation that heats gas also plays an important role in the final structure around the Lagrangian point. The instability of a vortex begins when the vorticity at its center is close to the background vorticity (e.g., Surville et al. 2016). A hotter disk would have a higher background pressure, reducing the difference (gradient) from the background to the center of the vortex, thus reducing the strength of the vortex. This is observed in model 3 (5 L_⊙), where an early instability appears when compared to model 4 (1 L_⊙).

5.4 Final remarks

We do not include the back-reaction from dust on gas, which can be relevant when the local gas-to-dust ratio is about ~ 1 : 1. In our simulations, we show that at the Lagrangian points, the gas-to-dust ratio could reach values of ~ 2 : 1 (Table 2). Also, in model 1, a large vortex at the outer disk accumulates a considerable amount of dust, reaching a gas-to-dust ratio of ~9:1 (Fig. 2).

The inclusion of the back-reaction of dust should include a term of the form $\varepsilon \frac{\Omega (r)}{St}\Delta v$ (similar to the drag force in Eq. (6)), where ɛ is the dust-to-gas ratio. Two-dimensional numerical simulations show that this feedback could reduce the lifetime of vortices when dust density becomes comparable to gas density (ɛ ~ 1) within the vortex, reducing dust trapping (Fu et al. 2014). However, when the dust density diminishes, the vortex can appear to decrease again due to the growing dust-to-gas ratio. This ensures efficient dust trapping over the lifetime of the disk (Raettig et al. 2015). Therefore, with or without a dust back-reaction, this should not affect our conservative estimation of 1–2% of the mass reservoir ($M_{\text{res}}^{L_4,L_5}$) being piled up around Lagrangian points. Also, for a long-term evolution model (1000 orbits), the vortices at the Lagrangian points disappear. However, the dust accumulated during the first stages remains trapped, making the in situ gravitational collapse of dust to form km-sized asteroids (or planetesimals) at those locations a plausible scenario.

6 Conclusions

We highlight our main findings as follows:

• When a planet carves a gap, it creates an overdensity in the trailing position (L₅) compared to the leading one (L₄), producing a large vortex at the former (larger pressure bump in L₅). The asymmetry is due to pressure waves from the planet and the action of viscosity, which together promote a slightly larger loss of angular momentum at L₄. This translates to lower density and pressure around L₄ with respect to L₅;

• L₄ accumulates more μm and submm particles, while L₅ efficiently traps larger grains (mm–cm). The total bulk mass is retained in the largest size particles, making L₅ always more massive;

• Most of the particles trapped in L₄ possess Stokes numbers St < 0.1, while particles trapped around L₅ have St > 0.1 up to 10;

• Planet migration does not appear to influence the observed asymmetry between L₄ and L₅. However, the free parameters of the model may affect it: lower viscosity tends to L₄ ~ L₅. Lower stellar irradiation (colder disks) tends to enhance dust accumulation in both L₄ and L₅. Larger planetary mass tends to destabilize both, but especially L₄ by enlarging the planetary gap;

• The mass reservoir around a Lagrangian point is limited to $M_{\text{res}}^{L_4,L_5}={\displaystyle {\int }_{r_{\text{p}}-R_{\text{Hill}}}^{r_{\text{p}}+R_{\text{Hill}}}{\sigma }_{\text{dust}}}(r)2\pi r\text{d}r$, where σ_dust(r) is the initial dust density distribution. The retained dust was initially located only at the co-rotating orbital path of the planet, inside the effective capture radius defined by Eq. (14);

• If Trojans appear to have formed at an early evolutionary stage of the solar nebula on short timescales (10⁵ yr), after which their chemical composition should have been similar. However, their subsequent evolution and composition may have been significantly altered by dynamical instabilities and interactions (as in the Nice model for instance).

Future observations of embedded planets in disks – in particular the hypothetical detection of thermal emission in co-rotation with the planet – will allow us to quantify the role of gas effects for dust trapping around Lagrangian points. This mechanism has profound implications for the formation of Trojans in our Solar System but also in extra-solar systems.

Acknowledgements

M.M. acknowledges financial support from the Chinese Academy of Sciences (CAS) through a CAS-CONICYT Postdoctoral Fellowship administered by the CAS South America Center for Astronomy (CASSACA) in Santiago, Chile. J.O. acknowledges financial support from the Universidad de Valparaíso and from Fondecyt (grant 1180395). C.A.G. acknowledges Mulatona Cluster from CCAD-UNC, which is part of SNCAD-MinCyT, Argentina. AB acknowledges support from FONDECYT Regular 1190748. N.C. acknowledges support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska–Curie grant agreement No. 210 021. M.M., J.G., J.C., J.O., A.B., and M.S. acknowledge support from Iniciativa Científica Milenio via the Núcleo Milenio de Formación Planetaria.

Appendix A Figures models: 2, 3, 4, and 6

For a better visualization, we regroup all the figures concerning models 2 (Figs. A.1, A.5), 3 (Fig. A.3), 4 (Fig. A.4), and 6 (Fig. A.2) in this appendix.

Fig. A.1 Dust mass accumulated over time around the Lagrangian points L4 (continuous line) and L5 (dotted line) for model 2. The evolution is shown for different sizes, including the effective (all sizes) accumulated mass (black continuous line).

Fig. A.2 Dust mass accumulated in time around the Lagrangian points L4 (continuous line) and L5 (dotted line) for model 6. Evolution is shown for different sizes, including the effective (all sizes) accumulated mass (black continuous line).

Fig. A.3 Dust mass accumulated over time around the Lagrangian points L4 (continuous line) and L5 (dotted line) for model 3. Evolution is shown for different sizes, including the effective (all sizes) accumulated mass.

Fig. A.4 Dust mass accumulated over time around the Lagrangian points L4 (continuous line) and L5 (dotted line) for model 4. Evolution is shown for different sizes, including the effective (all sizes) accumulated mass.

Fig. A.5 Mass spectrum of model 2 as a function of particle size for the dust accumulated in L4 and L5 for the last evolutionary stage (460 orbits). The bottom panel shows the mass ratio L4∕L5 as a function of particle size.

Fig. A.6 Mass spectrum of model 6 as a function of particle size for the dust accumulated in L4 and L5 for the last evolutionary stage (460 orbits). The bottom panel shows the mass ratio L4∕L5 as a function of particle size.

References

All Tables

All Figures

	Fig. 1 Equipotential lines according to Eq. (11) in a rotating frame with the planet, applied to a model with parameters M⋆ = 1 M⊙, rp = 7.8 au, and Mp = 5 MJ. The leading Lagrangian point (L4) is located at + π∕3, while the trailing point (L5) is at − π∕3. A stable librating azimuthal angle is assumed to be ±π∕8 around eachLagrangian point.
In the text

	Fig. 2 Evolution of the gas and dust in model (1), for three evolutionary stages; 46, 230, and 459 orbits. The dust panels show the full size range of the particles, from 10−4 to 100 cm. The panels have been rotated to always keep the planet at the same location x = 7.8 au, y = 0.
In the text

	Fig. 3 Top: Evolution of the vortensity and dust in model 1, for three evolutionary stages: 46, 230, and 459 orbits. Bottom: Dust distribution of particles with sizes in the range 10−2 − 100 cm. The planet rotates in the counter-clockwise direction. The panels have been rotated to always keep the planet at the same location, namely x = 7.8 au, y = 0.
In the text

	Fig. 4 Evolution in time of the dust for model 1 as a function of particle size, ranging from μm to cm. The panels have been rotated to always keep the planet at the same location x = 7.8 au, y = 0.
In the text

	Fig. 5 Mass spectrum of model 1 as a function of the size of the dust accumulated in L4 and L5 for thelast evolutionary stage (460 orbits). The bottom panel shows the mass ratio L4 ∕L5 as a function of dust particle size.
In the text

	Fig. 6 Mass spectrum of model 1 as a function of the Stokes number for the dust accumulated in L4 and L5 at the last evolutionary stage (460 orbits). The bottom panel shows the mass ratio L4∕L5.
In the text

	Fig. 7 Histogram of the effective dust mass accumulated at the end of the simulation inside L4 and L5, distributed by radial bins. This defines a capture region given by rp + 2RHill from the planet towards the outer region, and rp − 1.5RHill from the planet towards the star. An effective capture radius is expected at rp ± Rhill.
In the text

	Fig. 8 Dust mass accumulated in time around the Lagrangian points L4 (continuous line) and L5 (dotted line) for model 1. The evolution is shown for different sizes, including the effective (all sizes) accumulated mass (black continuous line).
In the text

	Fig. 9 Radial trajectory of randomly selected particles (μm to cm) that end up trapped in the second inner ring located at r = 5.6 au featured at the end of the simulation in model 1.
In the text

	Fig. 10 Trajectories of some particles (μm to cm) from model 1 that end up trapped inside the Lagrangian point L4 only. The top horizontal red lines represent the location L5, while the bottom horizontal lines represent L4. Some submm particles oscillate in horseshoe orbits (between both L4 and L5) before ending in L4. Large mm–cm particles oscillate only in tadpole orbits around L4.
In the text

	Fig. 11 Radial trajectory of randomly selected particles (μm to cm) that end up trapped in the outer ring located at r = 11 au featured at the end of the simulation model 1.
In the text

	Fig. 12 Spectral energy distribution from the disk, and the Lagrangian point L5 from model 1 computed at the last evolutionary stage. Emission from the Lagrangian point comes from an optically thin dusty swarmof ~3 Mmoon with a mean temperature of ~55 K.
In the text

	Fig. A.1 Dust mass accumulated over time around the Lagrangian points L4 (continuous line) and L5 (dotted line) for model 2. The evolution is shown for different sizes, including the effective (all sizes) accumulated mass (black continuous line).
In the text

	Fig. A.2 Dust mass accumulated in time around the Lagrangian points L4 (continuous line) and L5 (dotted line) for model 6. Evolution is shown for different sizes, including the effective (all sizes) accumulated mass (black continuous line).
In the text

	Fig. A.3 Dust mass accumulated over time around the Lagrangian points L4 (continuous line) and L5 (dotted line) for model 3. Evolution is shown for different sizes, including the effective (all sizes) accumulated mass.
In the text

	Fig. A.4 Dust mass accumulated over time around the Lagrangian points L4 (continuous line) and L5 (dotted line) for model 4. Evolution is shown for different sizes, including the effective (all sizes) accumulated mass.
In the text

	Fig. A.5 Mass spectrum of model 2 as a function of particle size for the dust accumulated in L4 and L5 for the last evolutionary stage (460 orbits). The bottom panel shows the mass ratio L4∕L5 as a function of particle size.
In the text

	Fig. A.6 Mass spectrum of model 6 as a function of particle size for the dust accumulated in L4 and L5 for the last evolutionary stage (460 orbits). The bottom panel shows the mass ratio L4∕L5 as a function of particle size.
In the text