© ESO 2021

1 Introduction

Measurements of the Venera 9 and 10 spacecraft (Krasnopolsky & Krysko 1979) served as first evidence for a co-orbital dust ring around Venus. These measurements were tested later by the Helios mission (Leinert & Moster 2007) and were confirmed with observations from STEREO (Jones et al. 2013, 2017). In Pokorný & Kuchner (2019) the total amount of material within the dust ring was estimated to be about 1.3 × 10¹³ kg, which corresponds to an equivalent of an asteroid of about 2 km in diameter (assuming a mean bulk density of 3 g cm⁻³). According to the authors, the analysis of the measurements of the STEREO spacecraft shows that the ring is spread within a region of a radial width equal 0.06 AU and a height of 0.1 AU. The question now is where it comes from. If the dust ring of Venus is not a temporary phenomenon (e.g., caused by comets passing the inner Solar System), it must be continuously filled up with new material either from the interplanetary medium or from the erosion of some primordial Venus Trojans (Pokorný & Kuchner 2019).

The capture probability of dust in co-orbital resonance with Venus is rather low (Jeong & Ishiguro 2012), but still possible (Jackson & Zook 1989), and it is also known to exist for the other planets, for example, Earth (Dermott et al. 1994) or Jupiter (Lhotka & Celletti 2015). It has been shown that the capture process into mean motion resonance (MMR) between a dust particle and a planet strongly depends on the order of the resonance. It is well known that capture within the orbit of the perturber is not possible but can take place for outer resonances (Sicardy et al. 1993; Beauge & Ferraz-Mello 1994; Liou & Zook 1995). Low capture probabilities are not the only problem that needs to be solved to explain the dust ring of Venus. Capture is also known to be a phenomenon, but it is only temporal (Lhotka & Celletti 2015; Lhotka & Gales, 2019). The reason is found in the effect of nongravitational perturbations, that is, solar wind drag and the force of the so-called Poynting-Robertson effect (Klačka et al. 2014): micron-sized dust particles in the Solar System experience a net drift toward the inner part of the Solar System through loss of orbital energy that leads to a shrinking of the semimajor axes. The drift is canceled in (outer or co-orbital) resonance for finite times, but because of nonlinear perturbations, the particle starts librating with increasing amplitude around the exact resonance condition, until it is removed and again starts drifting (in radial distance) toward the Sun.

We here mainly investigate the timescales of temporary capture of micron-sized dust grains that are in co-orbital resonance with Venus. We also provide a detailed analysis of the location of the exact resonance that depends, as we show, on various parameters: the size and the charge of the dust grains, the interplanetary magnetic field strength and so on. We focus in particular on some nongravitational effect that is usually neglected in this type of studies: the role of charge in the interplanetary magnetic field. The resulting Lorenz force may strongly affect the dynamics, as has previously been shown for micron-sized particles out of resonance (Lhotka et al. 2016) and in outer MMRs with Jupiter (Lhotka & Gales, 2019). While in the former study the authors find suitable conditions under which the drift in semimajor axes of micron-sized dust grains is reduced as a result of the normal component of the interplanetary magnetic field, a reduction in capture time has been found in Lhotka & Gales (2019) from the enhancement of chaotic regions close to resonance because of acharge. We therefore investigate the shift in equilibria and time of temporary capture in the case of co-orbital motions of charged dust grains with Venus, including the important effect of the interplanetary magnetic field.

We state the models and methods that we used in Sect. 2. The dynamical behavior of uncharged particles in 1:1 MMR is studied both analytically and numerically in Sect. 3. The effect of the Lorentz force resulting from the interplanetary magnetic field is described in detail in Sect. 4. We present a summary and the discussion of our study in Sect. 5. Further mathematical details related to the analytical study, including a more general formula for the shift in location of the co-orbital resonance caused by nongravitational effects and the second-order solution to the characteristic equation in linear stability analysis, can be found in the appendices.

2 Model and method

In the absence of nongravitational forces, the only perturbations of the dust grain orbits are due to additional planets and the problem reduces to the standard N-body problem. In the case of only one perturber, that is, Venus, the problem further reduces to the restricted three-body problem, which serves as a prototype model in Solar System dynamics. Depending on the geometry of the perturber’s (Venus) orbit the planar or inclined, circular or elliptic three-body problems may be defined. These models serve as the basis for our investigations (Dvorak & Lhotka 2013). However, when we consider dust grain dynamics in the Solar System, nongravitational effects can usually not be neglected (Burns et al. 1979). The dust grain orbits are strongly affected by the conjugated effects of the solar radiation and solar wind. Let r be the position of a dust grain in the heliocentric coordinate system. Assuming the dust grain to be a spherical object of radius R, let m be its mass, A = πR² the cross-sectional area exposed to sunlight, ρ the mass density, and q the charge. The parameters are related by m = 4πρR³∕3 and q = 4πε₀UR with surface potential U and electric constant ε₀. As we show below, the strength of the nongravitational effects is proportional to 1∕R (obtained from the area-to-mass ratio), and U∕R² (obtained from the charge-to-mass ratio), respectively. In this section we essentially follow the mathematical treatment in Lhotka & Gales (2019), where the equation of motion for a test particle is provided by Eq. (1) therein: $$m\stackrel{\mathrm{..}}{r}=F_{\text{g}}+F_{s/p}+F_{\text{L}}.$$(1)

Here F_g indicates the gravity of the Sun and additional bodies; F_s∕p is the force induced by the interaction of the dust grain with solar radiation and solar wind, and F_L is the Lorentz force. For a complete treatment of the model, see Lhotka & Gales (2019). We only report the formalism related to the nongravitational effects that we used.

2.1 Solarradiation pressure, the Poynting-Robertson effect, and solar wind drag

According to Klačka et al. (2012, 2014), the force caused by solar radiation pressure, the Poynting-Robertson effect, and solar wind drag can be expanded up to first order in v∕c as $$F_{s/p}=-\nabla \frac{\mu m\beta }{r}-\frac{\mu m\beta }{r^2}\left(1+s_w\right)\left[\frac{(\dot{r}\cdot r)r}{c{r}^2}+\frac{\dot{r}}{c}\right],$$(2)

where v is the velocity of the particle and c is the speed of light, ∇ is the gradient operator, μ is the mass of the Sun, β is the ratio between the magnitudes of the solar radiation pressure and gravity, ṙ is the velocity vector of the particle relative to the Sun, r indicates the distance to the Sun, and s_w is the ratio of solar wind drag to Poynting-Robertson effect.

The first term in Eq. (2), which corresponds to the solar radiation pressure, can be expressed as the gradient of a potential proportion to gravitational potential. Therefore we incorporate this term into the Kepler motion by assigning the Sun a reduced mass of (1 − β) μ when we consider the solar radiation pressure for particles. Because the force induced by the solar wind is similar to that from Poynting-Robertson effect but has a different coefficient s_w, we denote the combined effect of them as P-R-S effect hereafter.

2.2 Lorentz force

The interplanetary magnetic field has an Archimedian spiral field line pattern because the solar wind plasma takes the magnetic field radially away from the Sun, and in addition, the solar rotation winds the field lines around the rotation axis. The spiral field pattern is also referred to as the Parker spiral. The magnetic field can be oriented away from the Sun (the “away” sector) or toward the Sun (the “toward” sector), depending on how the coronal field lines are configurated. In interplanetary space near the ecliptic plane, the interplanetary magnetic field typically has a four-sector structure, which comes from an ecliptic section of the curved shaped of the magnetic equator of the Sun. In response to the solar magnetic field, the interplanetary magnetic field roughly exhibits a 22-yr cycle such that the magnetic polarity reverses in every solar cycle with 11 yr.

The Parker spiral model (Webb et al. 2010; Bieber et al. 1987) is adopted in this paper to describe the interplanetary magnetic field. We follow the simplified analytic expression in Lhotka & Gales (2019), $$B=\frac{B_0{r}_0^2}{r^2}\left(\frac{r}{r}-\frac{{\Omega }_s}{u_{\text{sw}}}{g}_{\tilde{z}}\times r\right)\tanh \left(\alpha \frac{r\cdot g_{\tilde{z}}}{r}\right),$$(3)

where B₀ is the background magnetic field strength given at the reference distance r₀, Ω_s indicates the solar rotation rate, and u_sw is the radial speed of the solar wind, α is a positive parameter determining the shape of the sign change function f(ϕ) = tanh(αϕ), and $g_{\tilde{z}}$ is the unit vector pointing to the rotation axis of the Sun and can be defined by $$g_{\tilde{z}}=\sin (i_0)\left[\sin ({\Omega }_0)g_x-\cos ({\Omega }_0)g_y\right]+\cos (i_0)g_z,$$(4)

where i₀ is the angle between solar equator and the ecliptic plane, and Ω₀ is the ecliptic angle difference between the direction of the vernal equinox and the line of nodes between the equatorial and ecliptic planes. Following Lhotka & Gales (2019), {g_x, g_y, g_z} represents the basis of heliocentric reference frame, where g_z points to the ecliptic pole, g_x is aligned with the direction to the vernal equinox, and g_y is determined by the right-handed rule.

Hence in the heliocentric ecliptic reference frame, the Lorentz force F_L is (Lhotka & Gales, 2019) $$\begin{array}{lll}{F}_{\text{L}}\hfill & =\hfill & -\frac{q{B}_0{r}_0^2}{r^2}\left[\frac{1}{r}r\times \dot{r}+\frac{{\Omega }_s}{r}r\times (r\times g_{\tilde{z}})+\frac{{\Omega }_s}{u_{\text{sw}}}(r\times g_{\tilde{z}})\times \dot{r}\right]\hfill \\ \hfill & \hfill & \times \tanh \left(\alpha \frac{r\cdot g_{\tilde{z}}}{r}\right).\hfill \end{array}$$(5)

In addition to the position and velocity, the acceleration caused by the Lorentz force is mainly dependent on the charge-to-mass ratio q∕m of particles.

We note that the interplanetary magnetic field given by Eq. (3) does not include time-dependent effects, such as variations in the axial tilt or the solar cycle. However, it is based on the classical Parker spiral including the axial tilt of the solar dipole with respect to the ecliptic (Lhotka & Narita 2019). The approach is therefore only valid on secular timescales, that is, where time-dependent effects due to solar activity are canceled because of averaging. For additional information, see Lhotka & Narita (2019).

2.3 Parameters and numerical setup

The values we used and the corresponding references for parameters involved with nongravitational effects are summarized in Table 1. The orbital elements of the planets at epoch of JD 245 1545.0 in the ecliptic system of J2000.0 are taken from the JPL HORIZONS system¹ (Giorgini et al. 1996). By convention, the orbital elements with respect to the ecliptic plane in the heliocentric reference frame, which are semimajor axis, eccentricity, inclination, argument of perihelion, longitude of the ascending node and mean anomaly, are denoted by a, e, i, ω, Ω, and M, respectively. We use the subscript “V” to indicate Venus throughout.

We investigated the dynamics of co-orbital particles under different magnitudes of nongravitational effects, which can be evaluated by β and q∕m. Setting ρ = 2.8 g cm⁻³, we have (Beauge & Ferraz-Mello 1994; Lhotka & Gales, 2019) $$\beta =0.205/R,\gamma =0.0094U/R^2,$$(6)

where the radius R is in microns and the surface potential U is in Volt, and γ is the numerical value of q∕m in C/kg. Given a surface potential of 10 V, the dependence of β and γ on the size of particle is depicted in Fig. 1. As we show in Sect. 3.3, the β value that we consider for the Venus dust ring does not exceed 0.34, which corresponds to the minimum particle radius of 0.6 μm.

In numerical simulations, we modified the symplectic massive body algorithm (SyMBA) integrator (Duncan et al. 1998; Levison & Duncan 2000) to include nongravitational effects for dust particles. SyMBA is a second-order symplectic integrator that can handle close encounters between objects and also a nonconservative effect. It is based on a variant of standard mixed variable symplectic (MVS) methods and could switch to an improved multiple time-step method when close encounters occur. Except for solar radiation pressure, which is incorporated into the gravitational potential, we included other nongravitational effects by simply adding corresponding acceleration in every time step.

In our analytical study, we considered the planar circular restricted three-body problem (CRTBP) in the barycentric rotating reference frame where the unit of mass is the total mass of the Sun and planet (Venus in this paper), and the unit of length is the semimajor axis of the planet. When G = 1, the unit of time should be T′∕(2π), where T′ is the period of the mean motion of the planet. The numerical value of the speed of light c in the normalized unit system is denoted by c_v and can be calculated by $$c_{\text{v}}\approx \frac{173T^{\prime }}{2\pi a^{\prime }}=\frac{173}{n^{\prime }{a}^{\prime }},$$(7)

where a′ and n′ are the semimajor axis and mean motion of the planet, and T′ is the corresponding period in days. For Venus, c_v is 8561.

The framework of an elliptic restricted three-body problem (ERTBP) in heliocentric reference frame is mainly adopted in our numerical simulations. Like the CRTBP, the ERTBP here also consists of the Sun, Venus, and a massless particle, but Venus is assumed to move in its real orbit with low eccentricity (0.00673) and inclination (3.39°) in the ecliptic system of J2000.0 centered on the Sun instead of on a circular orbit. We also considered the complete model of the current planetary configuration of our Solar System based on J2000.0 in our numerical simulations to study the effect of the planetary perturbations. All planets in the Solar System are included, and the Earth and Moon are treated as a whole body in their barycenter. All numerical simulations were carried out within the framework of the ERTBP, unless otherwise specified. Although the CRTBP is based on the barycentric reference frame, the orbital elements are always defined with respect to the Sun for all models in this paper.

Fig. 1 Parameter β (blue) and γ (orange) assuming a surface potential of 10 V. The dashed line shows R = 10 μm, for which we find β = 0.0205, γ = 0.00094.

3 Uncharged problem

The triangular Lagrangian points L₄ and L₅ are dynamically stable for all planets in our Solar System in the pure gravitational CRTBP (see e.g. Murray & Dermott 1999). In the vicinity of these points, an important reservoir for celestial objects such as asteroids and dust is therefore located. The investigation of these enormous populations could not only reveal the dynamical mechanisms behind the complicated orbital behavior, but also determine the possible origin of these objects and verify the scenarios proposed to describe the early stage of our Solar System.

The typical phase portraits of the 1:1 MMR with Venus obtained within the framework of the ERTBP on (σ, a) and (Δω, e) planes are shown in Fig. 2. Here, σ = λ − λ_V is defined as the resonant angle of 1:1 MMR where λ = ω + Ω + M is the mean longitude. We note for planar CRTBP that ω and Ω have no definition so that σ = M − M_V. We use Δω = ω − ω_V to denote the difference in argument of perihelion between the particle and Venus, and e is the eccentricity of the dust particle. We integrated 60 orbits for each plot in the framework of the ERTBP. Two stable points in the center of tadpole regimes, which are denoted by L₄ and L₅, reside symmetrically at 60° and 300°, respectively. The saddle point L₃, located at 180°, separates two tadpole regimes and the horseshoe regime. We note that the semimajor axis a in the plot is normalized by $a_{1:1}=a_{\text{V}}{(1-\beta )}^{1/3}$, which is the location of the 1:1 MMR for dust particles (see below). Here β = 0 and so a_1:1 = a_V. The apsidal difference Δω could librate around the center for orbits nearby. The center (Δω_k, e_k) with k = 4, 5 is (± 60°, e_V) for L₄ and L₅. The phase portraits on the (Δω, e) plane are also symmetric for L₄ and L₅ when β = 0, so that we only show the portrait for the former in Fig. 2.

Fig. 2 Phase portrait of the 1:1 MMR with Venus in the ERTBP. Upper panel: (σ, a) plane with resonant angle σ and (normalized) semimajor axis a, while lower panel: (Δω, e) plane (for L4 only) with apsidal difference Δω and eccentricity e.

3.1 Location of triangular Lagrangian points

The theoretical derivation in this section is based on the planar CRTBP. We aim to derive formulae to predict the shift of the Lagrangian points L₄ and L₅ due to the nongravitational effects, that is, radiation pressure and the P-R-S effect. We worked in synodic coordinates where r₁ and r₂ indicate the position vectors from the Sun and planet, respectively. The distance to the Sun r (in Eqs. (2), (3), and (5)) therefore equals r₁. In the rotating reference frame centered on the center of mass of the Sun-planet system with the planet fixed at the + x -axis, the equations of motion including F_s∕p in the planar case (Chernikov 1970; Schuerman 1980) can be written as $$\begin{array}{cc}\ddot{x}& -2\dot{y}+\frac{{\gamma }_c(1+s_w)}{r_1^2}\left\{\dot{x}+\left[\dot{x}(x+{\mu }^{\prime })+\dot{y}y\right]\frac{x+{\mu }^{\prime }}{r_1^2}\right\}\\ & =-\frac{(1-\beta )\mu (x+{\mu }^{\prime })}{r_1^3}-\frac{{\mu }^{\prime }(x-\mu )}{r_2^3}+x+\frac{{\gamma }_c(1+s_w)y}{r_1^2},\end{array}$$(8) $$\begin{array}{cc}\ddot{y}& +2\dot{x}+\frac{{\gamma }_c(1+s_w)}{r_1^2}\left\{\dot{y}+\left[\dot{x}(x+{\mu }^{\prime })+\dot{y}y\right]\frac{y}{r_1^2}\right\}\\ & =-\frac{(1-\beta )\mu y}{r_1^3}-\frac{{\mu }^{\prime }y}{r_2^3}+y-\frac{{\gamma }_c(1+s_w)(x+{\mu }^{\prime })}{r_1^2},\end{array}$$(9)

where μ′ = 1 − μ is the mass of the planet, (x, y) are the coordinates of the dust particle with respect to the barycenter of the Sun and Venus, r₂ indicates the distance to the planet, and γ_c is a parameter given by γ_c ≡ β μ∕c = β (1 − μ′)∕c, which mainly determines the magnitude of the P-R-S effect for particles on specific orbits (see Eq. (2)).

In addition to our analytical solutions to Eqs. (8) and (9), some numerical simulations within the framework of the ERTBP were conducted to verify the theoretical formulae. We also included other planets in the Solar System to study the effect of the planetary perturbations on the co-orbital dynamics for Venus.

3.1.1 Analytical study

When β≠0, nongravitational effects introduce additional forces that are related to the position and velocity of the particles. Therefore the equation of motion in the classic CRTBP should be modified, and the equilibrium solutions were changed accordingly. Because nongravitational effects are dependent on β, which is determined by the size of the particle, the locations of the triangular Lagrangian points rely on the size of particles.

In the rotating reference frame of the CRTBP, the locations of libration points considering solar radiation pressure and the P-R-S effect are determined by setting ẍ = ý = ẋ = ẏ = 0 in Eqs. (8) and (9). Following Schuerman (1980), we defined the following symbols: $$\begin{array}{lll}{r}_1^{(0)}\hfill & =\hfill & {(1-\beta )}^{1/3}=\delta ,\hfill \\ r_2^{(0)}\hfill & =\hfill & 1,\hfill \\ x^{(0)}\hfill & =\hfill & \frac{{(1-\beta )}^{2/3}}{2}-{\mu }^{\prime }=\frac{{\delta }^2}{2}-{\mu }^{\prime }=x_0,\hfill \\ y^{(0)}\hfill & =\hfill & \pm {(1-\beta )}^{1/3}{\left[1-\frac{{(1-\beta )}^{2/3}}{4}\right]}^{1/2}\hfill \\ \hfill & =\hfill & \pm \delta {\left(1-\frac{{\delta }^2}{4}\right)}^{1/2}=y_0,\hfill \end{array}$$

which are also the analytical solution of the locations of L₄ and L₅ when we ignore the P-R-S effect and consider only solar radiation pressure. We note that the plus sign in Eq. (13) corresponds to L₄ while the minus corresponds to L₅.

We multiplied Eq. (8) by y and Eq. (9) by (x + μ′), then subtracted them. Making use of ẍ = ý = ẋ = ẏ = 0 and recalling that $r_1^2={(x+{\mu }^{\prime })}^2+y^2$, we can obtain $$r_2={\left[1-\frac{{\gamma }_c(1+s_w)}{{\mu }^{\prime }y}\right]}^{-1/3}.$$(14)

To keep r₂ > 0 (it is always satisfied for L₅), we should ensure that for L₄ $${\gamma }_c<\frac{{\mu }^{\prime }y}{(1+s_w)}.$$(15)

Or equivalently, $$\beta <\frac{{\mu }^{\prime }yc}{(1+s_w)(1-{\mu }^{\prime })}.$$(16)

Applying Eq. (16) to Venus and taking y ≈ y₀, we can obtain β < 0.01356 for L₄. Substituting the expression of r₂, given by Eq. (14), in the right-hand side of Eq. (9) which is equal to 0, we obtain $$\left[{\gamma }_c(1+s_w)+(1-{\mu }^{\prime })y\right]r_1^3-{\gamma }_c(1+s_w)(x+{\mu }^{\prime })r_1={\delta }^3(1-{\mu }^{\prime })y.$$(17)

Since $$x+{\mu }^{\prime }=\frac{r_1^2-r_2^2+1}{2},$$(18)

we can rewrite Eq. (17) as $$\begin{array}{ll}\hfill & \left[\frac{{\gamma }_c(1+s_w)}{2}+(1-{\mu }^{\prime })y\right]r_1^3+\frac{{\gamma }_c(1+s_w)}{2}\hfill \\ \hfill & \times \left\{{\left[1-\frac{{\gamma }_c(1+s_w)}{{\mu }^{\prime }y}\right]}^{-2/3}-1\right\}{r}_1={\delta }^3(1-{\mu }^{\prime })y.\hfill \end{array}$$(19)

The Nth order solution to Eqs. (19) and (14) can be written as a series expansion in γ_c : $$\begin{array}{ll}\hfill & x^{(N)}={\sum }_{m=0}^N{p}_{4,5}^{(m)}(\beta ,{\mu }^{\prime }){\gamma }_c^m{(1+s_w)}^m,\hfill \\ \hfill & y^{(N)}={\sum }_{m=0}^N{q}_{4,5}^{(m)}(\beta ,{\mu }^{\prime }){\gamma }_c^m{(1+s_w)}^m,\hfill \end{array}$$(20)

where $p_{4,5}^{(m)}(\beta ,{\mu }^{\prime })$ and $q_{4,5}^{(m)}(\beta ,{\mu }^{\prime })$ are coefficients of the mth order term for L₄ or L₅. Obviously $p_4^{(0)}=p_5^{(0)}=x_0$ and $q_4^{(0)}=q_5^{(0)}=y_0$. Solutions up to the third order are derived in Appendix A.1. Because the solutions are complicated, we tried to simplify them by assigning β = 0 in $p_{4,5}^{(m)}(\beta ,{\mu }^{\prime })$ and $q_{4,5}^{(m)}(\beta ,{\mu }^{\prime })$ for m > 0 terms. Therefore, the coefficients $p_{4,5}^{(m)}$ and $q_{4,5}^{(m)}$ with m > 0 are only dependent on μ′. This is feasible because in the case of Venus, even the original formulae shown in Appendix sectionlinking A.1 are well consistent with the numerical solutions to the equations of motion only when β ≲ 0.01 for both L₄ and L₅ (see below). In this region we have small differences in the resonant angles of libration centers σ_k (k = 4, 5) between original and simplified formulae, up to 0.12° and 0.007° for L₄ and L₅, respectively. Up to the third order in γ_c, we have $$\begin{array}{ccccccc}{p}_4^{(1)}& =& -p_5^{(1)}& =& & -\frac{2-{\mu }^{\prime }}{3\sqrt{3}(1-{\mu }^{\prime })\mu },& \\ p_4^{(2)}& =& p_5^{(2)}& =& & \frac{(64-9{\mu }^{\prime }){\mu }^{\prime }-44}{162{(1-{\mu }^{\prime })}^2{{\mu }^{\prime }}^2},& \\ p_4^{(3)}& =& -p_5^{(3)}& =& & \frac{u[-3(102-{\mu }^{\prime }){\mu }^{\prime }+520]-248}{729\sqrt{3}{(1-{\mu }^{\prime })}^3{{\mu }^{\prime }}^3},& \\ q_4^{(1)}& =& q_5^{(1)}& =& & \frac{2-3{\mu }^{\prime }}{9{\mu }^{\prime }(1-{\mu }^{\prime })},& \\ q_4^{(2)}& =& -q_5^{(2)}& =& & \frac{-{\mu }^{\prime }(16-13{\mu }^{\prime })+4}{54\sqrt{3}{(1-{\mu }^{\prime })}^2{{\mu }^{\prime }}^2},& \\ q_4^{(3)}& =& q_5^{(3)}& =& & -\frac{3{\mu }^{\prime }[{\mu }^{\prime }(27{\mu }^{\prime }-58)+16]+40}{2187{(1-{\mu }^{\prime })}^3{{\mu }^{\prime }}^3}.& \end{array}$$(21)

Equations (12) and (13) show that the locations of L₄ and L₅ are symmetric for γ_c = 0. The distance from both L₄ and L₅ to the Sun is between 0 and 1 for 0 < β < 1. As the effect of solar radiation pressure, the triangular Lagrangian points approach the Sun for increasing values of β, while the distance from the planet remains unchanged. This causes a larger resonant angle for L₄ (> 60°) and a smaller resonant angle for L₅ (< 300°). However, when we take the P-R-S effect into consideration, Eq. (14) tells us that L₄ recedes from the planet while L₅ approaches (y > 0 for L₄ while y < 0 for L₅). This leads to the asymmetry between L₄ and L₅. Therefore, the combined effect of solar radiation pressure and the P-R-S effect always drives L₄ away from the planet. At the same time, we always have r₁ < δ from Eq. (19) for L₄. This causes the resonant angle of L₄ to increase rapidly with β. The situation is much more complicated for L₅. To the first order (see Eq. (A.3)), r₁ is always larger than δ, however, which causes the resonant angle to increase. In this case, solar radiation pressure competes with the P-R-S effect to slow down the increase or decrease of the resonant angle, which depends on the value of parameter β and c_v. We recall Eq. (7), and we find that c_v is inversely proportional to v_p∕c, where v_p is the velocity of the planet and is close to the velocity of particles v. Therefore a larger c_v corresponds to a smaller v∕c and thus a weaker P-R-S effect (see Eq. (2)). When c_v is large, solar radiation pressure is more likely to dominate so that the resonant angle of L₅ still decreases with β. For small enough c_v, however, the P-R-S effect could counteract the effect of solar radiation pressure or even dominate the motion, leading to an increase in resonant angle.

3.1.2 Application to Venus

We applied our formulae to Venus. For L₄ and L₅, in addition to the simplified formulae shown in Eq. (21) and the original formulae (first-order to third-order solutions) shown in Appendix A.1, we can also obtain (x, y) for L₄ and L₅ by either numerically solving Eqs. (19) and (14) via Newton’s method (denoted by “numerical solution” hereafter) or numerically integrating the equation of motion directly (denoted by “simulation” hereafter). Then the resonant angle σ_k is calculated with the formulae describing the transformation from the barycentric (CRTBP) or heliocentric (ERTBP) Cartesian coordinate system to heliocentric polar coordinate system. We implemented the Newton method in a Python code to compute the numerical solutions and adopted the SyMBA integrator for simulations. In numerical simulations, the positions of libration centers (σ_k, a_k) are determined by locating the centers of libration island in the phase portraits on the (σ, a) plane for different β, where the libration amplitudes Δσ (max[σ] −min[σ]) are minimum. Concretely, we sampled hundreds of orbits with initial conditions on the (σ₀, a₀) plane for each β and then integrated them for thousands of years. The other orbital elements were set to be the same as those of Venus. The intervals of β were set to be 10⁻⁴ and 0.002 for L₄ and L₅, respectively. Then we delimited the region of libration centers by comparing the libration amplitudes of these orbits. Afterward, the region enclosing the orbits with minimum libration amplitudes was located. Iteratively, we were able to sample and integrate some orbits in this reduced region and repeated the above steps until we reached a desired resolution. Because there is almost no difference in value between a_k and a_1:1, we always used a_1:1 to evaluate a_k for simplicity.

The resonant angles for different Lagrangian points σ_k from various solutions are shown in Fig. 3. Apparently, σ₄ and σ₅ always increase with β. The numerical solutions to Eqs. (19) and (14) agree perfectly with the result of numerical simulations based on the ERTBP. Both of them show that there are no solutions for β > 0.01135 for L₄ and β > 0.33865 for L₅ because they have collided with L₃ and L₁, respectively (see below). However, even the third-order approximation formula cannot determine the real libration center for large β, especially for L₅. Intuitively, this is because the parameter c_v (≈ 8561) in this case is not large enough, leading to a relatively large γ_c and slow convergence of the expansion. In other words, the high orbital velocity of Venus results in a relatively large v∕c, which corresponds to a relatively small c_v according to Eq. (7). This means that it is not very effective to treat the perturbation from the P-R-S effect as an expansion in γ_c for large β. At the same time, small c_v could also cause an increase of the resonant angle of L₅ because the P-R-S effect is more effective here, as we mentioned before.

Given corresponding initial positions (x, y), we can also numerically solve Eqs. (19) and (14) via Newton’s method and then calculate the resonant angle for L₁, L₂, and L₃. As shown in Fig. 3, σ₁ and σ₃ decrease with β. Finally, L₃ coincides with L₄ at β = 0.01135 with σ₃ = σ₄ = 108.4°. Similarly, L₁ coincides with L₅ at β = 0.33865 with σ₁ = σ₅ = −5.57°. For larger β, these Lagrangian points disappear and the corresponding equilibrium solutions no longer exist. L₂ is always located in the forth quadrant forany nonzero β and σ₂ varies in an extremely small range (<5 × 10⁻⁵ deg for 0 < β < 1).

Fig. 3 Resonant angle of the libration centers (σk) for L4 (top) and L5 (bottom) of Venus from different solutions. The positions of the libration centers for L3 and L1 obtained from the numerical solutions to the equations of motion are also presented in the upper and lower panels, respectively. “SIMP” represents the results calculated from the simplified approximation formulae Eqs. (20) and (21). The results from approximation formulae for L5 are only shown for 0 < β < 0.01 (we insert a window to zoom in) because their deviations from the numerical solutions are too large for β > 0.01.

3.2 Linear stability

Consider a dust particle that is given a small displacement (X, Y) from the equilibrium point (X₀, Y₀). The linearized equations of motion (see Appendix A.2) have the form $$\stackrel{.}{\text{X}}=\text{AX},$$(22)

where A is the matrix (a_ij), with a_ij (i, j = 0...3) indicating the linear coefficients with respect to X, Y, Ẋ, and Ẏ in Eqs. (A.15), (A.16), and $$\text{X}=\left(\begin{array}{c}X\\ Y\\ \dot{X}\\ \dot{Y}\end{array}\right).$$(23)

The characteristic equation of A can be written as $${\lambda }^4+a_3{\lambda }^3+(1+a_2){\lambda }^2+a_1\lambda +9{\mu }^{\prime }(1-{\mu }^{\prime })b(1-a_0)=0,$$(24)

where a_k (k = 0...3) are expressions obtained from a_ij and $$b=1-\frac{{\delta }^2}{4}.$$(25)

Neglecting the P-R-S effect, that is, setting γ_c = 0, we have a_i (i = 0...3) = 0. Then Eq. (24) becomes $${\lambda }^4+{\lambda }^2+9{\mu }^{\prime }(1-{\mu }^{\prime })b=0$$(26)

with the solutions $${\lambda }_0=\pm \text{i}{\left\{\frac{1}{2}\pm {\left[\frac{1}{4}-9{\mu }^{\prime }(1-{\mu }^{\prime })b\right]}^{1/2}\right\}}^{1/2}=\text{i}{\lambda }_r.$$(27)

In the absence of nongravitational effects except for solar radiation pressure, L₄ and L₅ are stable if $$36{\mu }^{\prime }(1-{\mu }^{\prime })b\le 1.$$(28)

This constraint is satisfied for all b in (0, 1) when $${\mu }^{\prime }\le \frac{1}{6}\left(3-2\sqrt{2}\right)\approx \mathrm{0.0286.}$$(29)

In the presence of the P-R-S effect, we can also solve Eq. (24) by expanding the parameters a_k in γ_c. The second-order solution is presented in Appendix A.2. We note that there always exists a positive real part for λ when γ_c > 0, which means that L₄ and L₅ are unstable when we consider the P-R-S effect.

Taking μ′≈ 0, the e-folding time of dust particles, which is the time interval in which the oscillation around the Lagrangian points increases by a factor of e, can be given by $$\frac{1}{T}\equiv \text{Re}(\lambda )=\frac{3{\gamma }_c(1+s_w)}{{\delta }^2}-\frac{7{\gamma }_c{}^2{(1+s_w)}^2\left(14-5{\delta }^2\right)}{2y_0{\delta }^2}.$$(30)

The first-order term in γ_c has been derived in Schuerman (1980), and we developed it to the second order. The second-order term is opposite for L₄ and L₅ because they have opposite y₀. Therefore the e-folding time for L₄ is slightly longer than that for L₅. We also extended the solution to the third order. However, μ′ = 0 is a singularity for the third-order solution, so that we cannot give a similar simplified formula for the e-folding time. The linear stability analysis suggests that the stability of triangular Lagrangian points, which is reflected by the e-folding time, declines with β. This qualitative result can be confirmed by numerical simulations. Quantitatively, however, the e-folding time is not necessarily equivalent to the lifespan of particles, which is the time duration of staying in the 1:1 MMR region because the linear stability analysis is conducted in a small neighborhood around the equilibrium point. For Venus, the lifespans are much longer than the e-folding time (see Sect. 3.3).

In Fig. 4 we present an example describing the orbital evolution of a particle in L₅ with β = 0.07. The initial conditions were set to be in the libration center (a_k, e_k, i_V, Δω_k, Ω_V, σ_k) = (0.706 AU, 0.0064, 3.39°, −38.45°, 76.68°, −25.93°) in the ERTBP (see Sec. 3.3). The orbit librated around the libration center with extremely small amplitudes at the beginning. After about 30 kyr, the variations in semimajor axis and resonant angle started to increase. Then the eccentricity changed dramatically, and finally, the orbit escaped at about 45 kyr. As mentioned in Lhotka & Gales (2019), the inclination and longitude of ascending node stay more or less unchanged in the uncharged problem before the orbits escape from the co-orbital region. According to Eq. (30), the e-folding time is about 3 kyr for β = 0.07 for L₅. It is much shorter than the lifespan (45 kyr), as we mentioned before.

3.3 Numerical study. ERTBP and N-body

To better illustrate how solar radiation pressure and P-R-S effect affect the dynamical behavior of orbits in 1:1 MMR, we constructed the phase portraits for different values of β on the (σ, a) plane. Figure 5 presents a comparison between β = 0.006 and β = 0.012. Sixty orbits initially equally spaced (Δσ₀ = 6°) along a = a_1:1 were integrated for 10 kyr for each plot within the framework of the ERTBP. The typical phase portrait on the (σ, a) plane without any nongravitational effect is shown in Fig. 2, where the locations of L₄ and L₅ are symmetric and L₃ is located at 180°. However, the presence of β brings in an asymmetry to the phase portrait. There exists a displacement of L₃, L₄ and L₅ for 0 < β < 0.01135, as the example of β = 0.006 shows. The value of σ corresponding to L₃ decreases to 153.7°, while L₄ (72.8°) and L₅ (308.0°) shift to the right. These data are obtained from numerical simulations (L₄ and L₅) and numerical solutions to the equations of motion (L₃). Apparently, they perfectly depict the locations of libration points and saddles in the phase portraits (indicated by vertical lines). The perturbation caused by the P-R-S effect gives rise to chaos because the dispersion of orbits in the phase portraits is visible, especially for those far away from triangular Lagrangian points. As mentioned in Sect. 3.1, for β > 0.01135, L₃ and L₄ merge and disappear. Along with this, the tadpole region of L₄ vanishes. This is clearly shown by the example of β = 0.012, in which the value of σ for L₅ is 313.6°. For all β, the positions of equilibrium points on the a axis are very close to a_1:1.

In Sect. 3.1 we have also shown a good agreement between the results from the analysis of the CRTBP and numerical simulations in the ERTBP (see Fig. 3) only with a slight difference between them (denoted by Δσ_k) that is due to the eccentricity of Venus. We calculated the difference in σ_k between the numerical solutions (in the CRTBP) and simulations (in the ERTBP). Figure 6 shows that σ_k in the CRTBP is always larger than that in the ERTBP, and this kind of difference increases with β generally. The eccentricity of Venus has a larger effect for L₄ (up to 0.13°) than L₅ (up to 0.025°). We also note smooth trends of Δσ_k versus β, especially for small β. We assume that an expansion similar to that in Sect. 3.1 could be conducted to solve the equations of motion within the framework of the ERTBP. Some additional term involved with β should appear due to the eccentricity of Venus.

We onlyconsidered the case of three-body problem so far. For 1:1 MMR with Venus, the perturbation from the other planets in the Solar System could play an important role in the dynamical evolution of dust orbits. In our paper about Venus (Trojans, in prep.), we conduct a comprehensive investigation of the resonance mechanism in the complete model that includes all planets of the Solar System. In Fig. 6 we also present the deviations in σ_k between thecomplete model and ERTBP. In contrast to the CRTBP, the planetary perturbations giving rise to the onset of chaos in the whole vicinity of the equilibrium points cause disordered deviations to σ_k. However, the planetary perturbations result in a comparable shift Δσ_k as the CRTBP.

As revealed by Fig. 5, in addition to the resonant angle of libration centers σ_k, the libration width of tadpole regions (δa) also varies with β. Obviously, the tadpole region of L₄ decreases with increasing values of β until it disappears completely. In contrast, thetadpole cloud of L₅ expands with β. In practice, we classified the orbits starting in the vicinity of L₄ that never exceedL₃ during the integration into the tadpole region of L₄. Similarly, the tadpole region of L₅ consists of the orbits starting in the vicinity of L₅ that never go below L₃. Many orbitswere integrated to precisely delimit the boundary of the tadpole region. In our simulations we used the values of σ for L₃ obtained from the numerical solutions to the equations of motion (blue curve in the top panel of Fig. 3) in order to classify the tadpole regions of L₄ and L₅. As we mentioned in Sect. 3.1, L₃ and L₄ coincide and disappear for β > 0.01135. In this case, we can no longer define the above-mentioned libration width. Finally, the libration width of each triangular Lagrangian point is defined as the maximum difference in semimajor axis of the orbits belonging to the corresponding tadpole region.

According to Murray & Dermott (1999), massless particles are thought to move on tadpole orbits if $\delta r⩽{(8{\mu }^{\prime }/3)}^{1/2}{a}_{\text{V}}=2.555\times {10}^{-3}{a}_{\text{V}}$, where δr is the radial separation of particles from Venus. This condition holds for both L₄ and L₅ if β = 0, as shown in Fig. 7. This figure also shows that the libration width of L₄ decreases with β, while the width of L₅ increases with β. The corrections are almost linear, especially for L₅, where the slope is lower (0.158 v.s. − 0.184 by linear fitting). The libration width of L₄ decreases to zero at β = 0.01135 and the tadpole region around L₄ vanishes. At the same time, δa of L₅ increases to 4.360 × 10⁻³ a_V.

The planetary perturbations could affect the size of tadpole regions by both changing the location of L₃ and affecting the dynamical behavior near the boundary of tadpole regions. The difference is unremarkable in most cases, however. The maximum deviation is less than 5 × 10⁻⁷ a_V. Considering that the libration width tends to 0 for L₄ when β tends to 0.01135, however, the planetary perturbations play an increasingly important role with increasing values of β.

When we include the horseshoe orbits, determining the extended libration width of the 1:1 MMR region (Δa) becomes more complicated because such orbits are chaotic due to the P-R-S effect. We set different criteria to describe the boundary of the 1:1 MMR region of Venus. Three different thresholds (340°, 350°, and 353°) indicating the maximum allowed libration amplitudes were chosen. The results shown in Fig. 8 suggest that the libration width of the 1:1 MMR increases with β first and then decreases. It is obvious that a higher threshold value of libration amplitudes produces a smaller turning point. In the region before the turning point, the libration width is heavily dependent on the threshold. This is due to the chaotic nature near the boundary, which also causes the roughness of the curve for a high threshold such as 353°. Meantime, the gradual rise of the width is caused by the expansion of the stable region around L₅ (see Fig. 5). After the turning point, the increasing perturbation due to the P-R-S effect gradually destabilizes the orbits in the 1:1 MMR region. Only the particles closely orbiting L₅ could survive, so the dependence of the libration width on the threshold is much weaker.

After we determined the positions of libration centers (σ_k, a_k), we studied the change in the phase plane (Δω, e) by modifying the value of β. In Fig. 9 we show on the (Δω, e) plane the evolution of 60 orbits with initial conditions σ = σ_k, a = a_k. According to Lhotka & Celletti (2015) and Lhotka & Gales (2019), (i, Ω) of L₄ and L₅ are always located at (i_V, Ω_V). These two elements remain unchanged during the integration and are not affected by solar radiation pressure and the P-R-S effect. In comparison with the typical phase portrait shown in Fig. 2 for β = 0, solar radiation pressure and the P-R-S effect exert additional forces in the orbital plane on dust particles and cause a deviation of the center (Δω_k, e_k). In addition, the orbits evolve spirally toward the center in the phase space regardless of the initial conditions. In the two cases we show in Fig. 9 (β = 0.004 for L₄ and β = 0.05 for L₅), the apsidal difference between the libration centers and Venus Δω_k approaches 0 and the eccentricity e_k gets lower than e_V in comparison with the case of β = 0, where Δω_k = ±60° and e_k = e_V.

Following the method we used to determine (σ_k, a_k) in Sect. 3.1, the centers of (Δω, e) were also precisely determined by minimizing the libration amplitudes in an iterative way. Figure 10 presents the result. For β = 0, the phase portrait is symmetric for L₄ and L₅ and the center (Δω_k, e_k), where the amplitudes of e and Δω are minimum, is located at the point of (e_V, ±60°). In the presence of a non-zero β, we find an asymmetry between L₄ and L₅. The eccentricity e_k with the smallest amplitudes monotonously decreases with β for L₄, at a much faster rate than for L₅ in the same range of β ∈ (0, 0.01135). Therefore a maximum β value of 0.01135 corresponds to a minimum e_k of 0.00235. For dust particles around L₅, although e_k could increase with β for 0.03< β < 0.245, it is still smaller than e_V in the whole range to be considered (0, 0.33865). In short, the combined effect of solar radiation pressure and the P-R-S effect always tends to reduce the eccentricity of dust particles of any size we considered in the triangular Lagrangian points. Despite a wider range of β, e_k of dust particles in L₅ varies in a narrower range with a minimum value of 0.00630.

Figure 10 shows that the apsidal difference Δω_k varies with β in a more sophisticated way than σ_k does (cf. Fig. 3). It always decreases first and then increases with β (for L₅ it even decreases again for β > 0.32). The (first) turning point is β ≈ 0.008 for L₄ and β ≈ 0.005 for L₅. We note that it is possible to observe dust particles of specific sizes moving on orbits with Δω_k = ±60°, just like the case where β = 0, but e_k < e_V. For L₄, this size is 21 μm, corresponding to a β of 0.009963, and it is 28 μm (β = 0.00734) for L₅.

Although the tadpole region of L₅ and the 1:1 MMR region may be enlarged under some conditions (see Figs. 7 and 8), the stability of triangular Lagrangian points is thought to be reduced for increasing β. To verify this, we integrated the orbits in the libration centers (a_k, e_k, i_V, Δω_k, Ω_V, and σ_k) for different β. The integration of one particle was ended when it escaped the 1:1 MMR region of Venus, and then we recorded the lifespan. After several test simulations, we applied a criterion that if an orbit dissatisfies |a −a_k| < 0.0075 AU at any time, it was regarded as escaped. It is feasible to set a threshold of 0.0075 AU, which is slightly higher than the maximum half-width of the 1:1 MMR region (see Fig. 8), because the orbits escaping the 1:1 MMR are scattered away or fall into the star in a short time. Generally, the lifespans shown in Fig. 11 decrease with β, indicating a diminishing stability. The trends of lifespans for L₄ and L₅ seem to be similar in the range of β < 0.01135. In the CRTBP, the triangular Lagrangian points are dynamically stable for all planets of the Solar System. When we included solar radiation pressure and the P-R-S effect, L₄ and L₅ were no longer stable. This was shown in the linear stability analysis in Sect. 3.2. As we also mentioned there, the lifespans here are several times longer than the e-folding time, but these two time scales follow a similar trend with β. Figure 11 also shows that the lifespan is about 38 Myr for orbits in L₄ with β = 10⁻⁴ and about 1.8 Myr for orbits in L₅ with β = 0.002. As β increases, the lifespan drops monotonously in general and decreases to 0.1 Myr when β = 0.03 (only for L₅).

Fig. 4 Orbital evolution of particles in L5 with β = 0.07 (R = 2.9 μm) based on the ERTBP. The orbit is initially placed in the libration center (see more details in the text).

Fig. 5 Phase portrait on the (σ, a) plane. Upperpanel: obtained for β = 0.006, while lower panel: obtained for β = 0.012. The vertical lines indicate σk obtained from numerical simulations (L4 and L5) and numerical solutions to the equations of motion (L3), which are 72.8° (L4), 153.7° (L3) and 308.0° (L5) for β = 0.006, 313.6° (L5) for β = 0.012.

Fig. 6 Difference of resonant angles of libration centers (Δσk) between different models and the ERTBP for L4 (top) and L5 (bottom). “pp” means the planetary perturbations.

Fig. 7 Libration width of tadpole regions for L4 (blue) and L5 (orange). The dashed line represents the result from the theoretical formula for β = 0. The interval of β for L5 is set to be the same as that of L4 in this case (10−4).

Fig. 8 Libration width of 1:1 MMR with Venus for different thresholds. See text for more details.

Fig. 9 Phase portrait on the (Δω, e) plane. Upperpanel: β = 0.004 for L4, and lower panel: β = 0.05 for L5. The horizontal and vertical lines indicate ek and Δωk obtained by minimizing the libration amplitudes, which are 0.00644 and 55.90° for the upperpanel, 0.006355 and − 43.24° for the lower panel.

Fig. 10 Eccentricity (blue) and apsidal difference (orange) with the smallest amplitudes for different values of β for L4 (top) and L5 (bottom). The blue horizontal line indicates the eccentricity of Venus (0.00673), and the orange horizontal line indicates ± 60°.

Fig. 11 Lifespans of dust particles in the libration centers of L4 and L5. A window is inserted to zoom in the region of 0 < β < 0.012.

4 Charged problem

In this section we reexamine the results of the previous section including the charge, that is, we focus on the role of the interplanetary magnetic field on captured charged dust grains in co-orbital resonance with Venus. Given the background magnetic field strength at specific distance, solar rotation rate and solar wind speed (see Table. 1), the Lorentz force shown in Eq. (5) is determined by the charge-to-mass ratio, position, velocity, and the deviation of the magnetic and orbital angular momentum axis.

4.1 Shift of resonance and capture time

In consideration of the important effect of the interplanetary magnetic field, especially for small particles (see e.g. Mann et al. 2014; Lhotka & Gales, 2019), we introduced the Lorentz force in our simulations. The charge-to-mass ratio q∕m (or γ) is related to β through Eq. (6), which can be written as γ ≈ 0.224 Uβ² assuming ρ = 2.8 g cm⁻³. We adopted the value of the maximum surface potential (U = 10 V) used in Mann et al. (2014) and Lhotka & Gales (2019) in this section, and thus we have γ ≈ 2.24 β².

A comparison of the orbital evolution between charged and uncharged particles is shown in Fig. 12. We switched on the Lorentz force for the orbit shown in Fig. 4. The corresponding value of γ is 0.011 for β = 0.07. We observe a shorter lifespan for the charged particle (~30 kyr). Except for this, however, the semimajor axis and resonant angle have a similar evolution for both particles, although the Lorentz force causes larger libration amplitudes at the very beginning. A close inspection also suggests that under the effect of the Lorentz force, deviations in the locations of libration centers in the phase space exist, especially in e and Δω. Furthermore, i and Ω are no longer constant during the integration because of the additional normal force exerted by the interplanetary magnetic field. The inclination can be pumped up to ~10°. The longitude of the ascending node follows a similar evolution as the inclination, librating with decreasing amplitudes. i and Ω are coupled, and their variation is confirmed to be related to the deviation of the magnetic axis from the angular momentum axis of the system, the initial inclination, and the ratio Ω_s∕n (Lhotka et al. 2016; Lhotka & Gales, 2019).

The resonant angle of libration centers σ_k depends on the Lorentz force. Because we fixed the surface potential, the Lorentz force is larger for smaller particles, whose β are also larger, according to γ ≈ 2.24 β². We repeated the iterative method in Sect. 3.1 to locate σ_k of libration centers. The upper panel of Fig. 13 shows that for L₄ the Lorentz force has a negligible effect on σ_k for β < 0.0085, corresponding to R > 24 μm and γ < 1.6 × 10⁻⁴. Even a higher value of γ can only cause a deviation nogreater than 0.02° in σ₄. This is principally due to the low limit of β (≤ 0.01135) required by the existence of L₄ (see Sect. 3.1), which corresponds to a small limit of γ (≤ 2.9 × 10⁻⁴) and thus a small Lorentz force. However, when a higher value of β is allowed, as in the case of L₅ where β ≤0.33865, the Lorentz force could become larger and affect the position of L₅ in a more effective way (see the lower panel of Fig. 13). The maximum difference Δσ_k could reach 2.4°, which corresponds to a distance of 4.5 × 10⁶ km from the location of the exact libration point in uncharged problem. In addition, for β >0.276, no stable points can be found in the phase portrait because of the chaos that is induced by the interplanetary magnetic field, which also causes the chaotic variation in Δσ_k with γ. Figure 13 also shows that in general, σ_k in the charged problem seems to be smaller than that in the uncharged problem for L₄ (negative Δσ₄) and larger than that in the uncharged problem for L₅ (positive Δσ₅) for most β.

Planetary perturbations might clearly enhance the chaotic behavior in the whole system. The isolated effect of the Lorentz force on the position of L₄ (see Fig. 13) is slightly smaller than that of planetary perturbations (see also Fig. 6). However, for L₅, of which the β range is much wider, the deviation coming from planetary perturbations is negligible in comparison of that from the Lorentz force for β > 0.054 (see the lower panel of Fig. 13). We also took into account the planetary perturbations and Lorentz force in our simulations. When these two forces are comparable in magnitude, the combined effect may produce a larger deviation in σ_k, just like in the case of L₄ (see the top panel in Fig. 13). However, when one of these two forces dominates, the combined effect follows the dominant effect, just like in the case of L₅ (see the lower panel in Fig. 13).

The Lorentz force could further enlarge the tadpole region of L₅ while reducing the tadpole region of L₄ (not shown here). In comparison with a maximum deviation of 3.5 × 10⁻⁷ a_V in the libration width of L₅, the Lorentz force has a much stronger effect on the tadpole region in L₄, which can be reduced by up to 8 × 10⁻⁶ a_V. Although the deviations are very small, this reduction could be really remarkable (up to 15%) because the libration width of L₄ for uncharged particles is also very small for large β.

We show in Fig. 14 the apsidal difference from Venus (Δω_k) and the eccentricity (e_k) of the libration centers in a charged problem (shown in dots) obtained with the same method mentioned in Sect. 3.3. Apparently, they differ from the results in the uncharged problem (shown in curves) under the effect of the Lorentz force. We conclude that the Lorentz force tends to increase Δω_k but decreases e_k for both L₄ and L₅ in most situations. For L₄, γ is limited to 3 × 10⁻⁴ due to the maximum value of β 0.01135. Therefore the perturbation from the interplanetary field is small enough for us to determine (Δω_k, e_k) precisely. Conversely, the orbits inL₅ with β > 0.05 are much more chaotic. The lifespans of such orbits are reduced, and the libration centers cannot be located precisely, especially for β > 0.244, for which (Δω_k, e_k) cannot be determined because the lifespans are too short. However, we can still demonstrate that the Lorentz force can greatly change e_k and Δω_k. The maximum differences could be up to 0.0022 and 25° for e_k and Δω_k, respectively.

Because many secular resonances control the dynamical behavior of particles in 1:1 MMR with Venus, the eccentricity and apsidal difference could experience complicated secular evolution. We therefore not show the plot of (Δω_k, e_k) under the effect of the planetary perturbations.

The lifespans of charged particles in the libration centers shown in Fig. 15 again confirm that the perturbation from the interplanetary magnetic field further destabilizes the co-orbital region. Considering the uncertainties in determining the lifespans, we find a significant decrease in lifespans for β >0.0105 for both L₄ and L₅. The lifespan difference is about 10⁵ yr for β = 0.01135, although the Lorentz force is smaller for smaller γ, but the longer lifespan makes it easier to accumulate the effect. On the other hand, the Lorentz force can also significantly reduce the stability within the short lifespan of large dust particles. The planetary perturbations can only cause a very limited difference to the lifespans compared to the Lorentz force. The most important dynamical effects caused by other planets on the co-orbital region are secular resonances (see, e.g., Brasser & Lehto 2002; Zhou et al. 2009, 2011, 2019, 2020; Dvorak et al. 2012), which need a long timescale to be effective. Before this, however, dust particles are very likely to have escaped from the 1:1 MMR as a result of nongravitational effects.

Fig. 12 Orbital evolution of particles initially located in L5 for β = 0.07. Orange curves indicate the model including the Lorentz force with γ = 0.011 (surface potential = 10 V), and the blue curves indicate the same orbit as in Fig. 4, where γ = 0. Both orbits are initially placed in the libration center of the uncharged problem.

Fig. 13 Differences of resonant angles of libration centers Δσk between different models and the ERTBP for L4 (top) and L5 (bottom). “Charged” refers to the model including the Lorentz force.

Fig. 14 Same as Fig. 10, but dots indicate the model including the Lorentz force given a surface potential of 10 V.

4.2 Parameter space study

In this section we investigate the effect of charge on the full parameter space (β, γ) with 0 ≤ β ≤ 0.01135 for L₄ and 0 ≤ β ≤ 0.33865 for L₅ and γ derived from Eq. (6), with 0 ≤ U ≤ 10 Volt (for L₄ and L₅) with the aim to provide the domain in β × γ for which charge cannot be neglected in dynamical studies of dust in the Solar System. The results are based on the ERTBP with the standard Parker model for the interplanetary magnetic field. The initial conditions for the numerical simulations are as follows: for fixed pairs of parameters (β, γ), we chose initial conditions centered at the equilibria of L₄ and L₅, taken from our previous analysis, together with initial conditions with slight deviations of δa = ±10⁻³ (normalized units), δσ = ±1°, δΔω = ±1°, or δe = ±e_k (with equilibrium value e_k for L₄ and L₅, respectively). This gives nine initial conditions in total for each pair (β, γ), but this is sufficient to investigate the effect of the charge parameter γ on (i) the time of temporary capture and (ii) the maximum libration amplitudes compared to the uncharged case defined with γ = 0. Let T_γ be the time of temporary capture for a fixed pair of values (β, γ) normalized by the time of temporary capture for the same β, but with γ = 0.

The results for the case L₅ are shown in Fig. 16, where we report the value of T_γ with 0 ≤ T_γ ≤ 1 in the parameter space (β, γ). As expected, for lower values of γ, we find T_γ > 0.5 up to β≃ 0.3 with T_γ close to one for small β and decreasing with increasing β. We also provide test cases within the green region in Fig. 17 with three orbits each for values (β, γ) in the dark green region (a), the light green region (b), and in the white region (c). We note that even low values of γ are sufficient to obtain T_γ ≃ 0 for β >0.3. With increasing values of γ, we find decreasing T_γ for the full range of β ≤ 0.3 up to the threshold γ ≃ 0.01 where T_γ < 0.5. Within the light red regime of Fig. 16 we still find capture, with much shorter capture times, while in the dark red regions the perturbations due to charge dominate the other effects. For comparison we also provide orbits within the red regions (d, e) in Fig. 17. We clearly see the reduction of the capture time caused by charge and perturbations of the interplanetary magnetic field. We conclude that the effects of the interplanetary magnetic field further destabilize motion in 1:1 resonance with Venus compared with vanishing charge-to-mass ratios in case of L₅.

The situation is somewhat different for the leading Lagrange point L₄. Already the standard nongravitational effect, that is, the P-R-S effect, destabilizes the orbits of captured dust in the uncharged case, and no capture persists for values with β > 0.01135. For surface potentials U ≤ 10 V the effect of γ is therefore not sufficiently large to act on the dynamics, that is, the capture processes. A similar study in the case of L₄ reveals that in the full parameter space (β, γ), we find T_γ ≃ 1 for all parameters and initial conditions (however, charge, i.e., the effect of the interplanetary magnetic field, still affects the maximum libration width; this is not shown here).

Fig. 15 Lifespans of dust particles in the libration centers of L4 and L5 in the charged problem. A window is inserted to zoom in the region of 0 < β < 0.012. The gray dots indicate the lifespans in the uncharged problem that we showed in Fig. 11.

Fig. 16 Effect of charge-to-mass ratio γ on stability (i.e., capture) time Tγ in the parameter space (β × γ).

Fig. 17 Effect of different charge-to-mass ratios on the stability time Tγ. Parameters for orbits a–e are taken from Fig. 16 with initial conditions at the libration point L5.

5 Summary and conclusions

We investigated the role of nongravitational effects on the dynamics of charged micron-sized dust grains in the co-orbital region of Venus. Our study includes the investigation of planetary perturbations, solar radiation pressure, the Poynting-Robertson effect, and solar wind drag (P-R-S effect). In addition, we performed a detailed study of the effect of the Lorentz force as a result of theinteraction with the interplanetary magnetic field (Parker spiral model), which is usually neglected in this type of studies (e.g., Beauge & Ferraz-Mello 1994; Liou & Zook 1995; Jeong & Ishiguro 2012; Lhotka & Celletti 2015).

Our results are based on analytical estimates in the framework of simplified dynamical models, that is, the circular restricted three-body problem, and numerical studies in the elliptic restricted three-body problem, and the planetary problem including all nongravitational perturbations. We used a Fortran code that has been used in Lhotka & Gales (2019) as well as SyMBA (Duncan et al. 1998; Levison & Duncan 2000), which we modified to include the additional perturbations, that is, the Lorentz force term.

We mainly focused on studying the position of the triangular Lagrange points L₃ (collinear), L₄ (leading), and L₅ (trailing) that define the geometry of the co-orbital regime of motions. We provided a detailed analysis of (i) the shift of the Lagrange points and (ii) the time of temporary capture of charged and uncharged dust grains in the co-orbital regime of motions, both due to the nongravitational effects. The shift and capture time strongly depend on the area-to-mass ratio of the dust grains, parameterized by the parameter β and the charge-to-mass ratio q∕m parameterized by the parameter γ.

As expected, we find that in the uncharged problem (γ = 0), the shift due to solar radiation pressure is symmetric with respect to L₄ and L₅, while the P-R-S effect acts asymmetrically on the location of L₄ and L₅. This phenomenon has previously been investigated in case of Jupiter (Lhotka & Celletti 2015). We derived an analytical formula to predict the shift for uncharged dust grains, based on the circular restricted three-body problem following the line of argumentation in Schuerman (1980). We generalized the results to higher orders in the small parameters and also performed a linear stability analysis that provides an order-of-magnitude estimate of the e-folding time. We confirmed that the nongravitational effects render the equilibria unstable.

These results serve as the basis for finding good parameters and initial conditions for massive numerical simulations in the framework of the inclined elliptic restricted three-body problem and the full planetary problem. The aim is to find orbits with minimum libration amplitudes close to L₄ and L₅ that define the geometry of the co-orbital regime of motion also in these more realistic dynamical models. Denoting by a_k, σ_k, e_k, and Δω_k, with k = 4, 5 the semimajor axes, resonant angle, orbital eccentricity, and difference in argument of perihelion (with respect to Venus) with minimum libration amplitudes, close to the triangular Lagrange points, we provided their dependence on β as well as γ≠0. We find that nongravitational perturbations destroy librational orbits in the vicinity of L₄ already at β ≃0.01135 (collision of L₃ with L₄), while librational orbits survive close to L₅ up to β ≃ 0.33865 (collision of L₁ with L₅). This difference, caused by the P-R-S effect, results in an asymmetry in the shift of L₄ and L₅ (with respect to the pure gravitational problem) and in different extensions of the regimes of librational motion of dust close to L₄ and L₅. The tadpole region expands with β for L₅ almost linearly, but for L₄ it shrinks. For the width of the co-orbital region, which also contains the horseshoe orbits, we can always find some turning point (β ~ 0.07) where the width is the largest.

The interaction of charged dust particles (γ≠0) with the interplanetary magnetic field contributes to the shift of the tadpole points as well, but to a much lesser extent. However, for large enough γ, we also finda notable effect that is comparable to the effect of the planetary perturbations in the uncharged case. The role of charge is also more important for the trailing Lagrange point L₅ that ensures librational motion up to much higher values in β that correspond to much smaller particles with much higher charge-to-mass ratios for fixed values of density and surface potential. The Lorentz force therefore also affects the libration points in the co-orbital regime of motions in an asymmetric way, supporting librational motion of tiny dust grains close to L₅. Even more importantly, the Lorentz force term will act on the orbital planes of motions. In the uncharged problem, the dust grains stay within the orbital plane of Venus, neglecting the additional perturbations of the planets. However, the normal force component due to the Lorentz force will bring the dust particles out of their initial plane with latitudinal libration amplitudes up to several degrees already in the circular restricted problem. This phenomenon has been studied in full detail in Lhotka & Gales (2019) and is related to the deviation of the magnetic axis of the sun from the axis of orbital angular momentum. Deviations from the orbital plane of Venus may reach several degrees. An additional important effect of the interplanetary magnetic field is the strong effect on the capture time of particles in the co-orbital regime of Venus.

It is commonly accepted that the orbital lifespan of (uncharged) dust grains in the Solar System is limited by the Poynting-Robertson effect. Outside of MMR with the planets, the effect leads to circularization and shrinking in semimajor axes, bringing the dust grains toward the Sun on spiral-like motions (Burns et al. 1979). However, the net effect that leads to the decrease in semimajor axis can vanish in presence of outer and co-orbital resonance with a planet, with the result of temporary capture of dust at specific distances from the Sun. This capture may strongly extend the lifespan of dust grains in the Solar System. In a detailed study of time of temporary captures of charged and uncharged dust grains in the co-orbital regime with Venus, we find a strong dependence on the area-to-mass ratio, related to parameter β, that is already known. However, we also see a strong dependence of the capture time on the parameter γ; the charge-to-mass ratio of the dust particle as well. The reason is found in the faster increase of the libration amplitudes of charged dust in comparison to uncharged one. While for 0 ≤ β ≤ 0.3 the capture time (in the vicinity of L₅) is about the same for γ = 0 and γ ≤ 0.01, even low values of γ≠0 may lead to release from co-orbital resonance at very short times. For higher charge-to-mass ratios γ > 0.01 the time of temporary capture may be reduced by a factor 2–10 depending on the actual value of γ. Finally, we find that the effect of the interplanetary magnetic field on capture times close to L₄ and L₅ is asymmetric as well; it is much more important in the vicinity of L₅.

We remark that our study is based on a simplified model of the interplanetary magnetic field and assuming time-independent (equilibrium) charge-to-mass ratios of the dust grains. This is an oversimplification of the real situation in our Solar System. Time-dependent effects on the polarity and magnitude of the interplanetary magnetic field and the role of local perturbations due to coronal mass ejections have not been investigated in the present study at all. Time-dependent charge due to changes in the local plasma environment in space may result in additional effects that we have not touched on either.

However, the motivation of this study was to provide details on the role of the mean interplanetary magnetic field on the dynamics that is usually neglected by other authors. It may also serve as a basis for future studies that include more realistic models of the heliosphere. To conclude, the role of charge and the interplanetary magnetic field on the orbital motion of dust in the Solar System has been shown to be important to better understand the distribution and timescales of dust in our Solar System.

Acknowledgements

This work has been supported by the Austrian Science Fund (FWF) in the framework of the stand alone project P30542, National Natural Science Foundation of China (NSFC, Grants No. 11473016, No. 11933001) and National Key R&D Program of China (2019YFA0706601). L.Z. acknowledges the financial support from the program of China Scholarship Council (No. 201906190106) for his visit in University of Vienna.

Appendix A: Locationshift of co-orbital resonances in the Solar System

A.1 Locationof triangular Lagrangian points

To the first order in γ_c, Eqs. (19) and (14) can be simplified as $$\left[\frac{{\gamma }_c(1+s_w)}{2}+(1-{\mu }^{\prime })y\right]r_1^3={\delta }^3(1-{\mu }^{\prime })y,$$(A.1) $$r_2=1+\frac{{\gamma }_c(1+s_w)}{3{\mu }^{\prime }y}.$$(A.2)

Assuming that $y=y_0+ϵ{\gamma }_c+o({\gamma }_c^2),$ where ϵ is the coefficient of the first-order term in γ_c, we can obtain r₁ and r₂ in the form $$\begin{array}{lll}{r}_1^{(1)}\hfill & =\hfill & \delta \left[1-\frac{{\gamma }_c(1+s_w)}{6(1-{\mu }^{\prime })y_0}\right],\hfill \\ r_2^{(1)}\hfill & =\hfill & 1+\frac{{\gamma }_c(1+s_w)}{3{\mu }^{\prime }{y}_0}.\hfill \end{array}$$

Substituting r₁ and r₂ with Eqs. (A.3) and (A.4) in Eq. (18), we obtain $$x^{(1)}=x_0\left[1-\frac{(1-{\mu }^{\prime })+{\mu }^{\prime }{\delta }^2/2}{3(1-{\mu }^{\prime }){\mu }^{\prime }{x}_0{y}_0}(1+s_w){\gamma }_c\right].$$(A.5)

We note that Eq. (A.5) is a corrected version of Eq. (19) in Schuerman (1980), and the first equation in Eq. (4) of Liou & Zook (1995) because both of them missed (1 − μ′), which is the mass fraction of the star, in the denominator. This only causes a negligible difference because the planet mass is far lower than the stellar mass, however, and thus (1 − μ′) ≈ 1. Recalling $r_2^2={(x+{\mu }^{\prime }-1)}^2+y^2$ and making use of Eqs. (A.4) and (A.5), we can obtain y to first order in γ_c: $$y^{(1)}=y_0\left\{1+\frac{{\delta }^2\left[(1-{\mu }^{\prime })-{\mu }^{\prime }\left(1-{\delta }^2/2\right)\right]}{6(1-{\mu }^{\prime }){\mu }^{\prime }{y}_0^3}(1+s_w){\gamma }_c\right\}.$$(A.6)

Equation (A.6) agrees perfectly with Eq. (20) in Schuerman (1980) and the second equation in Liou & Zook (1995).

Similarly, assuming $y=y_0+{ϵ}_1{\gamma }_c+{ϵ}_2{\gamma }_c^2+o({\gamma }_c^3)$ with coefficients ϵ_1,2 for the corresponding terms and making use of Eq. (A.6), we can solve Eqs. (19) and (14) up to the second order in γ_c, $$r_1^{(2)}=r_1^{(1)}+\frac{{\gamma }_c^2{(1+s_w)}^2\left[{\delta }^2(2-{\mu }^{\prime })-4(1-{\mu }^{\prime })\right]}{9\left(4-{\delta }^2\right)\delta {(1-{\mu }^{\prime })}^2{\mu }^{\prime }{y}_0^2},$$(A.7) $$r_2^{(2)}=r_2^{(1)}-\frac{{\gamma }_c^2{(1+s_w)}^2\left[6-4{\mu }^{\prime }-(2-{\mu }^{\prime }){\delta }^2\right]}{9\left(4-{\delta }^2\right)(1-{\mu }^{\prime }){{\mu }^{\prime }}^2{y}_0^4},$$(A.8)

or, equivalently, $$\begin{array}{lll}{x}^{(2)}\hfill & =\hfill & x^{(1)}-{\gamma }_c^2{(1+s_w)}^2\hfill \\ \hfill & \hfill & \times \frac{{\delta }^4{{\mu }^{\prime }}^2-4{\delta }^2\left(2{{\mu }^{\prime }}^2-4{\mu }^{\prime }+5\right)+16\left({{\mu }^{\prime }}^2-5{\mu }^{\prime }+4\right)}{72\left(4-{\delta }^2\right){(1-{\mu }^{\prime })}^2{{\mu }^{\prime }}^2{y}_0^2}\hfill \end{array}$$(A.9) $$\begin{array}{cc}{y}^{(2)}=y^{(1)}-& \frac{{\gamma }_c^2{(1+s_w)}^2}{144{(1-{\mu }^{\prime })}^2{{\mu }^{\prime }}^2{y}_0^5}\cdot \{{\delta }^4\left(5{{\mu }^{\prime }}^2{y}_0^2+4\right)\\ & +2{\delta }^2\left[\left(-11{{\mu }^{\prime }}^2+8{\mu }^{\prime }-10\right)y_0^2+4\right]+8{{\mu }^{\prime }}^2{y}_0^2\left(2y_0^2-1\right)\}.\end{array}$$(A.10)

Up to the third order in γ_c, we have $$\begin{array}{lll}{r}_1^{(3)}\hfill & =\hfill & r_1^{(2)}+\frac{{\gamma }_c^3{(1+s_w)}^3}{10368\delta {(1-{\mu }^{\prime })}^3{{\mu }^{\prime }}^2{y}_0^7}\hfill \\ \hfill & \hfill & \times [3{\delta }^{10}{{\mu }^{\prime }}^2-6{\delta }^8\left(5{{\mu }^{\prime }}^2-8{\mu }^{\prime }+10\right)\hfill \\ \hfill & \hfill & +48{\delta }^6\left({{\mu }^{\prime }}^2-4{\mu }^{\prime }+3\right)\hfill \\ \hfill & \hfill & -96{\delta }^4\left(1-{{\mu }^{\prime }}^2\right)-192\left(3{{\mu }^{\prime }}^2-8{\mu }^{\prime }+5\right)y_0^4\hfill \\ \hfill & \hfill & -32{\delta }^2{y}_0^2\left(12{\mu }^{\prime }+{{\mu }^{\prime }}^2{y}_0^2-12\right)],\hfill \end{array}$$(A.11) $$\begin{array}{lll}{r}_2^{(3)}\hfill & =\hfill & r_2^{(2)}+\frac{{\gamma }_c^3{(1+s_w)}^3}{5184{(1-{\mu }^{\prime })}^2{{\mu }^{\prime }}^3{y}_0^9}\times \hfill \\ \hfill & \hfill & \{3\left(4-{\delta }^2\right){\delta }^6{\left[\left(4-{\delta }^2\right){\mu }^{\prime }-2\right]}^2\hfill \\ \hfill & \hfill & -3{\delta }^2{y}_0^2[{\delta }^6{{\mu }^{\prime }}^2+16{\delta }^2{(2-{\mu }^{\prime })}^2-32\left(1-{{\mu }^{\prime }}^2\right)\hfill \\ \hfill & \hfill & -2{\delta }^4\left[{\mu }^{\prime }(5{\mu }^{\prime }-8)+10\right]]+896{(1-{\mu }^{\prime })}^2{y}_0^6\hfill \\ \hfill & \hfill & +192{\delta }^2(1-{\mu }^{\prime })y_0^4\left[\left(4-{\delta }^2\right){\mu }^{\prime }-2\right]\},\hfill \end{array}$$(A.12) $$\begin{array}{lll}{x}^{(3)}\hfill & =\hfill & x^{(2)}-\frac{{\gamma }_c^3{(1+s_w)}^3}{10368{(1-{\mu }^{\prime })}^3{{\mu }^{\prime }}^3{y}_0^7}\hfill \\ \hfill & \hfill & \times \{-3{\delta }^{10}{{\mu }^{\prime }}^3+6{\delta }^8{\mu }^{\prime }\left(6{{\mu }^{\prime }}^2-9{\mu }^{\prime }+10\right)\hfill \\ \hfill & \hfill & +12{\delta }^6(1-{\mu }^{\prime })\left[-{\mu }^{\prime }(20-9{\mu }^{\prime })+10\right]\hfill \\ \hfill & \hfill & +64(1-{\mu }^{\prime })\left[-{\mu }^{\prime }(35-4{\mu }^{\prime })+40\right]y_0^4\hfill \\ \hfill & \hfill & +48{\delta }^4\left[{\mu }^{\prime }\left({\mu }^{\prime }\left(y_0^2-10\right)+16\right)-6\right]\hfill \\ \hfill & \hfill & +32{\delta }^2[6{(1-{\mu }^{\prime })}^2(1+{\mu }^{\prime })+{{\mu }^{\prime }}^3{y}_0^4-\hfill \\ \hfill & \hfill & 6(5-{\mu }^{\prime })(1-{\mu }^{\prime })y_0^2]\},\hfill \end{array}$$(A.13) $$\begin{array}{lll}{y}^{(3)}\hfill & =\hfill & y^{(2)}-\frac{{\gamma }_c^3{(1+s_w)}^3}{81{\left(4-{\delta }^2\right)}^4{\delta }^4{(1-{\mu }^{\prime })}^3{{\mu }^{\prime }}^3}\hfill \\ \hfill & \hfill & \times [{\delta }^8{{\mu }^{\prime }}^3-2{\delta }^6\left(3{{\mu }^{\prime }}^3+45{{\mu }^{\prime }}^2-120{\mu }^{\prime }+80\right)\hfill \\ \hfill & \hfill & +8{\delta }^4\left(3{{\mu }^{\prime }}^3+78{{\mu }^{\prime }}^2-204{\mu }^{\prime }+115\right)\hfill \\ \hfill & \hfill & -32{\delta }^2\left(5{{\mu }^{\prime }}^3+33{{\mu }^{\prime }}^2-96{\mu }^{\prime }+56\right)\hfill \\ \hfill & \hfill & +384\left({{\mu }^{\prime }}^3-4{\mu }^{\prime }+3\right)].\hfill \end{array}$$(A.14)

To further verify the approximation formulae, we also applied them to Jupiter, where c_v ≈ 22 947 (the figure is not shown here). We find a perfect consistency with a maximum deviation of 0.007° for 0 < β < 0.5 for all formulae from first order to third order. The simplified formulae (see Eqs. (20) and (21)) are even less accurate than the first-order formulae in this case because the range of β is larger than in the case of Venus, which causes large errors in the derivation of $p_{4,5}^{(m)}$ and $q_{4,5}^{(m)}$ by assigning β = 0. In addition to the simplified formulae, the first- to third-order formulae are therefore still necessary in a study of nongravitational effects, especially for high values of β. Furthermore, for Jupiter, the resonant angle for L₄ σ₄ is also found to increase with β, while the resonant angle for L₅ σ₅ decreases with β as a result of dominating solar radiation pressure, as we mentioned in Sect. 3.1.

A.2 Linear stability of triangular Lagrangian points

We considered a small displacement (X, Y) from the equilibrium point (X₀, Y₀) in the planar CRTBP including solar radiation pressure and the P-R-S effect. Substituting (x, y) by (X₀ + X, Y₀ + Y) in Eqs. (8) and (9) and expanding in a Taylor series up to the first order of (X, Y, Ẋ, Ẏ), we deduce the linearized equations². $$\begin{array}{lll}\hfill & \hfill & \ddot{X}+\dot{X}{\gamma }_{c^{\prime }}\left[\text{}\frac{{(X_0+{\mu }^{\prime })}^2}{r_1^4}+\frac{1}{r_1^2}\right]+\dot{Y}\left[\frac{{\gamma }_{c^{\prime }}{Y}_0(X_0+{\mu }^{\prime })}{r_1^4}-2\right]\hfill \\ \hfill & \hfill & =X[1-\frac{{\delta }^3(1-{\mu }^{\prime })}{r_1^3}+\frac{3{\delta }^3(1-{\mu }^{\prime }){(X_0+{\mu }^{\prime })}^2}{r_1^5}-\frac{2{\gamma }_{c^{\prime }}{Y}_0(X_0+{\mu }^{\prime })}{r_1^4}\hfill \\ \hfill & \hfill & -\frac{{\mu }^{\prime }}{r_2^3}+\frac{3{\mu }^{\prime }{(X_0+{\mu }^{\prime }-1)}^2}{r_2^5}]\hfill \\ \hfill & \hfill & +Y\left[\text{}\frac{{\gamma }_{c^{\prime }}}{r_1^2}+\frac{3{\delta }^3(1-{\mu }^{\prime })Y_0(X_0+{\mu }^{\prime })}{r_1^5}-\frac{2{\gamma }_{c^{\prime }}{Y}_0^2}{r_1^4}+\frac{3{\mu }^{\prime }{Y}_0(X_0+{\mu }^{\prime }-1)}{r_2^5}\text{}\right],\hfill \end{array}$$(A.15) $$Y\frac{{\gamma }_{c^{\prime }}{Y}_0(X_0+{\mu }^{\prime })}{r_1^4}+2]+\dot{Y}{\gamma }_{c^{\prime }}\left(\frac{Y_0^2}{r_1^4}+\frac{1}{r_1^2}\right)\begin{array}{lll}\hfill & \hfill & =X[-\frac{{\gamma }_{c^{\prime }}}{r_1^2}+\frac{2{\gamma }_{c^{\prime }}{(X_0+{\mu }^{\prime })}^2}{r_1^4}+\frac{3{\delta }^3(1-{\mu }^{\prime })Y_0(X_0+{\mu }^{\prime })}{r_1^5}\hfill \\ \hfill & \hfill & +\frac{3{\mu }^{\prime }{Y}_0(X_0+{\mu }^{\prime }-1)}{r_2^5}]\hfill \end{array}\begin{array}{lll}\hfill & \hfill & +Y[1-\frac{{\delta }^3(1-{\mu }^{\prime })}{r_1^3}+\frac{2{\gamma }_{c^{\prime }}{Y}_0(X_0+{\mu }^{\prime })}{r_1^4}+\frac{3{\delta }^3(1-{\mu }^{\prime })Y_0^2}{r_1^5}\hfill \\ \hfill & \hfill & -\frac{{\mu }^{\prime }}{r_2^3}+\frac{3{\mu }^{\prime }{Y}_0^2}{r_2^5}].\hfill \end{array}$$(A.16)

We note that r₁ and r₂ above should be evaluated at the equilibrium points. The linear coefficients with respect to X, Y, Ẋ, Ẏ are denoted by a_ij (i, j = 0...3).

Up to second order in γ_c, we define the parameters $a_i(i=\mathrm{0..3})={\eta }_1{\gamma }_c+{\eta }_2{\gamma }_c^2+o({\gamma }_c^3)$ where η_1,2 are coefficients of the first and second terms (see Sect. 3.2). When γ_c = 0, a_i (i = 0..3) = 0, and thus Eq. (24) degenerates to Eq. (26) with solutions equal to λ₀ (see Eq. (27)). Therefore the solution to Eq. (24) can be written as λ = λ₀ (1 + c₁ + ic₂), where $c_i(i=1,2)={\zeta }_1{\gamma }_c+{\zeta }_2{\gamma }_c^2+o({\gamma }_c^3)$ with coefficients ζ_1,2 for corresponding terms. By substituting a_i (i = 0...3) and λ in Eq. (24) by series in γ_c, we can derive solutions to c₁ and c₂, $$\begin{array}{lll}{c}_1\hfill & =\hfill & \frac{{\gamma }_{c^{\prime }}}{4{\lambda }_r^4-2{\lambda }_r^2}\cdot \{-\frac{3{\left[a_{21}{\lambda }_r^2+9b{a}_{01}(1-{\mu }^{\prime }){\mu }^{\prime }\right]}^2}{2{\left(1-2{\lambda }_r^2\right)}^2}\hfill \\ \hfill & \hfill & +\frac{{\left[a_{21}{\lambda }_r^3+9b{a}_{01}{\lambda }_r(1-{\mu }^{\prime })\mu \right]}^2}{{\left(2{\lambda }_r^2-4{\lambda }_r^4\right)}^2}\hfill \\ \hfill & \hfill & +\frac{a_{21}\left[a_{21}{\lambda }_r^2+9b{a}_{01}(1-{\mu }^{\prime })\mu \right]}{2{\lambda }_r^2-1}\hfill \\ \hfill & \hfill & \frac{3a_{31}{\lambda }_r^2\left(a_{31}{\lambda }_r^2-a_{11}\right)}{2-4{\lambda }_r^2}-\frac{{\left(a_{11}-a_{31}{\lambda }_r^2\right)}^2}{4{\left(1-2{\lambda }_r^2\right)}^2}\hfill \\ \hfill & \hfill & +\frac{3{\left(a_{11}{\lambda }_r-a_{31}{\lambda }_r^3\right)}^2}{2{\left(1-2{\lambda }_r^2\right)}^2}+\frac{a_{11}\left(a_{11}-a_{31}{\lambda }_r^2\right)}{2-4{\lambda }_r^2}\hfill \\ \hfill & \hfill & +9b{a}_{01}(1-{\mu }^{\prime })\mu +a_{21}{\lambda }_r^2\hfill \\ \hfill & \hfill & -9b{a}_{02}{{\mu }^{\prime }}^2+9b{a}_{02}{\mu }^{\prime }+a_{22}{\lambda }_r^2\},\hfill \end{array}$$(A.17) $$\begin{array}{lll}{c}_2\hfill & =\hfill & \frac{{\gamma }_{c^{\prime }}}{2{\lambda }_r{\left(2{\lambda }_r^2-1\right)}^3}\cdot [a_{11}(-18b{\gamma }_{c^{\prime }}{a}_{01}{{\mu }^{\prime }}^2+18b{\gamma }_{c^{\prime }}{a}_{01}{\mu }^{\prime }\hfill \\ \hfill & \hfill & +{\gamma }_{c^{\prime }}{a}_{21}-4{\lambda }_r^4+4{\lambda }_r^2-1)+4a_{31}{\lambda }_r^6\hfill \\ \hfill & \hfill & +9b{\gamma }_{c^{\prime }}{a}_{01}{a}_{31}{{\mu }^{\prime }}^2-9b{\gamma }_{c^{\prime }}{a}_{01}{a}_{31}{\mu }^{\prime }\hfill \\ [-11pt]\hfill & \hfill & \hfill \\ \hfill & \hfill & -{\gamma }_{c^{\prime }}{a}_{12}{\left(1-2{\lambda }_r^2\right)}^2+2{\gamma }_{c^{\prime }}{a}_{21}{a}_{31}{\lambda }_r^4\hfill \\ \hfill & \hfill & -2{\gamma }_{c^{\prime }}{a}_{21}{a}_{31}{\lambda }_r^2-4a_{31}{\lambda }_r^4+a_{31}{\lambda }_r^2\hfill \\ \hfill & \hfill & +4{\gamma }_{c^{\prime }}{a}_{32}{\lambda }_r^6-4{\gamma }_{c^{\prime }}{a}_{32}{\lambda }_r^4+{\gamma }_{c^{\prime }}{a}_{32}{\lambda }_r^2].\hfill \end{array}$$(A.18)