Trojans in the Solar Neighborhood

TL;DR

This paper investigates whether the Hercules stream in the solar neighborhood could be composed of stellar trojans captured at the L4 Lagrangian point of the galactic bar, using high-resolution N-body simulations to analyze their dynamics.

Contribution

It demonstrates that stellar trojans at L4 can explain the Hercules stream's velocity asymmetry and provides testable predictions for Gaia data to confirm this origin.

Findings

Trojans cause the velocity asymmetry in Hercules.

Trojans are only temporarily captured at L4 before escaping.

The Hercules stream's distribution can be explained by trojan dynamics.

Abstract

About 20% of stars in the solar vicinity are in the Hercules stream, a bundle of stars that move together with a velocity distinct from the Sun. Its origin is still uncertain. Here, we explore the possibility that Hercules is made of trojans, stars captured at L4, one the Lagrangian points of the stellar bar. Using GALAKOS--a high-resolution N-body simulation of the Galactic disk--we follow the motions of stars in the co-rotating frame of the bar and confirm previous studies on Hercules being formed by stars in co-rotation resonance with the bar. Unlike previous work, we demonstrate that the retrograde nature of trojan orbits causes the asymmetry in the radial velocity distribution, typical of Hercules in the solar vicinity. We show that trojans remain at capture for only a finite amount of time, before escaping L4 without being captured again. We anticipate that in the kinematic plane the Hercules stream will de-populate along the bar major axis and be visible at azimuthal angles behind the solar vicinity with a peak towards L4. This test can exclude the OLR origin of the Hercules stream and be validated by Gaia DR3 and DR4.

Tables

Table 1. Table 1 Fraction of Hercules stars measured at various patches along the solar circle for a long and a short-fast bar.

$\varphi$	% Hercules	% Hercules
	Long bar	Short bar (R${}_{OLR}{}^{}/$R${}_{\odot }{}^{}=0.9$)
0o	6%	10%
-28o	18%	14%
-90o	25%	16%

Table 2. Table 2 CR: co-rotation resonance, ILR: Inner Lindblad resonance,OLR: Outer Lindblad resonance, IUH: Inner ultra-harmonic, OUH: outer ultra-harmonic.

t [Gyrs]	${\Omega }_P$ [km s-1kpc-1]	CR [kpc]	ILR [kpc]	OLR [kpc]	IUH [kpc]	OUH [kpc]
1.5	55	4.0	1.0	7.8	2.2	6.0
2.5	40	6.0	1.5	10.8	3.7	8.5
3.0	35	6.8	1.5	11.5	4.0	9.1
3.5	30	7.5	1.5	12.0	5.0	10

Equations

$$S_R={\int }^R\dot{R}dR$$

$$S_{\varphi }=L_z\varphi$$

$$J_R=\frac{1}{2\pi }\int \frac{\partial S}{\partial R}dR$$

$$J_{\varphi }=\frac{1}{2\pi }\int \frac{\partial S}{\partial \varphi }d\varphi =L_z$$

$${\Phi }_{eff}=\Phi -\frac{1}{2}{\Omega }_P^2{R}^2$$

$$m({\Omega }_{\varphi }-{\Omega }_P)=\pm l{\Omega }_r$$

$$\vec{v}=-\vec{\nabla }{\Phi }_{eff}+2\vec{\Omega }\wedge \vec{v}$$

$$J_{fast}=J_R-\frac{l}{m}{J}_{\varphi }.$$

$${\rho }_{\mathrm{DM}}=\frac{{\mathrm{M}}_{\mathrm{DM}}}{2\pi }\frac{a}{r(r+a{)}^3}$$

$${\rho }_d(R,z)=\frac{{\mathrm{M}}_{\mathrm{d}}}{8\pi z_0{\mathrm{R}}_{\mathrm{d}}^2}{\mathrm{e}}^{-\mathrm{R}/{\mathrm{R}}_{\mathrm{d}}}{\mathrm{sech}}^2(\mathrm{z}/{\mathrm{z}}_0)$$

$${\sigma }_r(R)=Q\frac{3.36G\Sigma (R)}{\kappa (R)}$$

$${\Phi }_S(u,v)=\frac{U(u)-V(v)}{{\sinh }^2\mathrm{u}+\mathrm{si}{\mathrm{n}}^2\mathrm{v}}$$

Full text

1. Introduction

The Gaia satellite is mapping the phase-space of the Milky Way, revealing proper motions and distances for a few millions of stars in the solar neighborhood. The Gaia DR2 observations confirmed that stars in the solar vicinity do not have a smooth distribution of velocities: about 3.5 million stars within 1.5 kpc of the Sun have a rich kinematic substructure that manifest as ridges in the velocity space (Gaia Collaboration et al. 2018; Ramos et al. 2018).

In particular, the distribution of stars in the velocity space defined by their radial and azimuthal velocities is dominated by a number of streams, bundles of stars that move together in the same direction with velocities that are distinct from neighboring stars (Antoja et al. 2018).

In addition to the traditional phase-space analysis used in previous works to identify moving groups of stars in the solar vicinity, more recent investigations include the action-space analysis. The phase-space of star motions is six-dimensional and commonly described by stars’ positions and velocities. The radial action, $J_{R}$ quantifies the amount of oscillation that the star exhibits in a radial direction along its orbit. In particular, this quantity describes the orbit’s eccentricity. The azimuthal action $J_{\phi}=L_{z}$ is the angular momentum on the vertical component $z$. In axisymmetric potentials, all of these actions are integrals of motions; that is, they are constant along the orbits and therefore can be used to label different stellar orbits (Binney & Tremaine 2008). In the presence of strong perturbations, such as those caused by spiral arms and bars, the actions are not fully conserved (see e.g. Vera-Ciro & D’Onghia 2016).

At a given point in time, however, the actions can still be calculated and may represent a diagnostic of non-axisymmetric perturbations (Contopoulos et al. 1973). This is because long-lived perturbations like spiral arms and the stellar bar can lead to an orbital diffusion that forms distinct features in orbit space. Therefore, in the solar vicinity, streams and structures in the velocity space manifest as structures and ridges in the action-space (Trick et al. 2018; Monari et al. 2018).

The origin of the structure in the phace-space is uncertain. Some works in the past pointed toward and extragalactic origin. Thus, Helmi et al. (1999) suggested that substructure in the kinematic plane was caused by accreted and disrupted satellites. However, the dynamical origin is the most accepted theory for the substructure in the kinematic plane. As previously noted, these features cannot simply have arisen from groups of stars that were born with similar kinematics, because these streams have various ages and metallicities, as confirmed in detailed studies (Famaey et al. 2005; Kushniruk et al. 2017), and most likely have a dynamical origin, e.g. caused by non-axisymmetric features of the stellar disk like the stellar bar and spiral arms.

While the Milky Way is a barred Galaxy (de Vaucouleurs 1964; Weinberg 1992) with a spiral structure, there is still debate on the structure and morphology of the stellar bar, its orientation with respect to the Sun, and its rotation rate, usually referred to as its pattern speed. Also uncertain is the number and strength of the spiral arms that propagate in the Galaxy (Reid et al. 2014). The existence of perturbations in the Milky Way in the form of the stellar bar and spiral arms ensures that stars may have their orbits perturbed when they move at the same frequency as these features, a process called resonance. For a stellar bar, there are three strong resonances: i) the co-rotation resonance characterized by star’s orbit moving at the same azimuthal frequency as the bar; ii) the inner-; and iii) the outer-Lindblad resonances, characterized by two radial oscillations for every azimuthal oscillation in the frame rotating with the Galactic bar pattern. The same resonances can be defined for spiral patterns.

Previous models of the inner Galaxy suggest that the radius of the outer-Lindblad resonance (OLR) of the Galactic bar lies in the vicinity of the Sun and affects the velocity distribution of nearby stars (Kalnajs 1991; Dehnen 2000). These studies have shown that for a model resembling the stellar disk with a bar with pattern speed $\Omega_{P}=55$ km s*-1* kpc*-1*, the OLR causes a bimodal distribution of the velocity between a moving group of low-velocity stars centered on the Local Standard of Rest (LSR) and an association of stars, the Hercules stream, moving outward and rotating more slowly than the Sun. The scenario proposed for the formation of the Hercules stream of stars matches the data if the Galactic bar has a radius of 3 kpc and is a fast rotator (Monari et al. 2017; Fragkoudi et al. 2019).

Recent measurements based on the three-dimensional density of red-clump giants, however, show that the Galactic bar extends to $\sim$5 kpc from the Galactic center (Wegg et al. 2015), in agreement with previous estimates (Benjamin et al. 2005) with a pattern speed of $\Omega_{P}\sim 39$ km s*-1* kpc*-1* (Portail et al. 2017). These values were recently confirmed by using a modified version of the Tremaine & Weinberg (1984) method by combining proper motions of Gaia DR2 and the VVV survey (Sanders et al. 2019; Clarke et al. 2019). If this is the case, then the OLR is placed at  10.5 kpc from the Galactic center, too far for the Hercules stream of stars to be caused by orbits migrating outwards at this resonance (Pérez-Villegas et al. 2017). The former study showed that for a bar with pattern-speed of  39 km s*-1kpc-1* the Hercules stream consists of stars trapped at co-rotation resonance with the bar.

Other studies reproduce distributions of stars with a similar degree of substructure in the kinematic space by invoking a sequence of short-lived spiral transients (De Simone et al. 2004; Quillen & Minchev 2005; Quillen et al. 2018), but the Hercules stream is not always well matched.

In an attempt to better reproduce the observed positions and velocities (namely the phase-space) of the stars in the Hercules stream, more recent models have included both the stellar bar and the spiral structure of the Milky Way. These studies, however, in some cases are based on N-body simulations that reproduce the solar vicinity by stacking several snapshots owed to the lack the spatial resolution of a few hundred parsecs (Fragkoudi et al. 2019). Other studies adopted semi-analytical models in which the spiral structure and the bar are modeled as externally applied perturbations with assumed properties (Hunt et al. 2018; Monari et al. 2018).

Here, following Pérez-Villegas et al. (2017), we approach the question of the formation and evolution of the Hercules stream in the long-bar scenario. We set up a high-resolution N-body simulation of the Milky Way with structural parameters that reproduce the currently observed properties of our Galaxy (Bland-Hawthorn & Gerhard 2016) and with enough spatial resolution to resolve the solar vicinity. This includes a self-consistent spiral pattern and bar whose properties are in agreement with current observations. The advantage of this methodology is that the spiral waves and the bar arise from the stars themselves as self-consistent collective disturbances (D’Onghia et al. 2013; D’Onghia 2015). In fact, if self-consistency is not an integral part of the model, the amplitude of the perturbation and the time-dependence assumptions made in the model may bias the resonant interactions between the bar, the spiral arms, and the stars, which are responsible for the formation of the Hercules stream.

The structure of this paper is as follows: we describe our methodology in Section 2. In Section 3 we present our results by showing the evolution of the Hercules stream by time. An example of the orbit of a star trapped in co-rotation with the bar and its escape are also illustrated. Finally, validating tests of our model are provided.

2. Methods

The simulations were carried out with GADGET3, a parallel TreePM-Smoothed particle hydrodynamics (SPH) code developed to compute the evolution of stars and dark matter, which are treated as collisionless fluids. The six-dimensional phase space is discretized into fluid elements that are computationally realized as particles in the simulations. A detailed description of the code is available in the literature. Here we note its essential features.

GADGET3 is a cosmological code in which the gravitational field on large scales is calculated with a particle-mesh algorithm, while the short-range forces are computed using a tree-based hierarchical multipole expansion, resulting in an accurate and fast gravitational solver. This scheme combines the high spatial resolution and relative insensitivity to clustering of tree algorithms with the speed and accuracy of the particle-mesh method to calculate the long-range gravitational field.

Pairwise particle interactions are softened with a spline kernel of scale-length $h_{s}$, so that they are strictly Newtonian for particles separated by more than $h_{s}$. The resulting force is roughly similar to traditional Plummer softening with scale-length $\epsilon=h_{s}/2.8$. For our applications the gravitational softening length is fixed to 40 pc for the dark halo and 28 pc for the stellar disk and 80 pc for the bulge. The total number of N-Body particles employed in the simulation (stellar disk, bulge, and halo) is approximately 90 million (see the Appendix for the details).

2.1. Solar vicinity and velocity definitions

In this simulation, named GALAKOS, the Sun is located at 8.1 kpc from the Galactic center at the azimuthal angle of -28o with respect the semi-major axis of the bar (see Figure 1). In cartesian coordinates its in-plane location is identified at (x,y)=(7.15, -3.8) kpc. Within 300 pc at the Sun’s location there are more than 4,000 stars with -200 $\leq$ z $\leq$ 200 pc. The stellar bar rotates counter-clockwise, and thus the Sun moves behind it. After 2.5 Gyrs of the life of the Milky Way, the stellar bar is in place with a length of approximately 4.5 kpc and a pattern speed of 40 km s*-1kpc-1* (see the Appendix for details), in full agreement with current data (Bland-Hawthorn & Gerhard 2016). Note that in the disk there are two strong arms extending from the edge of the bar to 5 scale lengths in radius ($\sim 11$ kpc). As shown in Figure 1 the stellar disk shows deviations from being axisymmetric with an asymmetry in the amplitude of the two spiral arms dominating the disk.

2.2. Action-angle space

In an axisymmetric disk, the energy $E$ and the angular momentum $L_{z}$ are the classical integrals for the in-plane motion. Actions are an alternative set of integrals of motion. For motion in the symmetry plane of an axisymmetric potential $\Phi(R,z)$ the equation of the action function $S$ separates in two terms, one azimuthal and one radial respectively: $S=S_{\phi}(\phi)+S_{R}(R)$ (Contopoulos et al. 1973).

Then

[TABLE]

and

[TABLE]

The radial action is:

[TABLE]

where the integral is taken over the full period. The azimuthal action is:

[TABLE]

$J_{\phi}=L_{z}$ is the specific angular momentum of the star around the Galactic center. Angles, $\theta_{\phi}$ and $\theta_{R}$ can be conveniently defined and are the conjugate variables to the actions. Angles increase at uniform rate: $\dot{\theta}_{\phi}=\Omega_{\phi}$, the azimuthal frequency and $\dot{\theta}_{R}=\Omega_{R}$, the radial frequency.

When the angle and action variables are used $J_{\phi}$ may be considered as the radial area of the orbit, $J_{R}$ as the extent of the oscillation, $\theta_{\phi}$ as pointing at the angular position of the epicenter around the Galaxy and $\theta_{R}$ as pointing at the phase of the epicyclic oscillation about that epicenter.

We used the AGAMA software package (Vasiliev 2019) to extract the potential from each snapshot of our simulation and compute the actions and angles (see Appendix for details).

3. Results

3.1. The velocity plane

An inspection of the radial and azimuthal velocity of nearby stars in Galactocentric coordinates shows that the solar neighborhood is far from smooth in its velocity distribution. Figure 2 displays an underlying smooth component of stars but the entire stellar distribution shows evidence of substructures in kinematics, which are streams of stars with coherent velocities. Note that the Hercules stream is reproduced in our model as a structure extending with azimuthal velocity Vϕ=180-220 km s*-1and radial velocity $-50\leq V_{r}\leq 100$ km s-1*. According to this definition 18% of stars at the solar vicinity in GALAKOS populate the Hercules stream. Thus, as compared to the Sun’s motion, Hercules is moving slower and shows an asymmetry pronounced in the radial velocity distribution with more stars moving away from the Galactic center as compared to the ones moving inwards as indicated also by the data (Dehnen 2000; Antoja et al. 2018). In what follows, we investigate whether this kinematic feature in the velocity space is the result of dynamical processes, such as the effects of the bar and spiral arms, which are fully captured in our simulation.

3.2. Donkey orbits

In the frame of a rotating bar, the effective potential is given by:

[TABLE]

with $\Phi$ the total potential and $\Omega_{P}$ the pattern speed of the bar. In barred galaxies, the Lagrangian points–defined as points in the plane where the effective potential is minimal or maximal–play an important role. In general, bars have two unstable Lagrangian points, L1 and L2, at the minimum potential shortly beyond their ends, along the major axis, and two other points, L4 and L5, along the intermediate axis at the maximum of potential (labeled in Figures 3-4). The situation is similar to that encountered in the restricted three-body problem, except that there is no equilibrium point L3. Such a characteristic means that rotation is an essential factor in bar dynamics (de Vaucouleurs 1964; Pfenniger 1990).

The disk of the Galaxy, as described by an axisymmetric exponential radial profile, is dominated by orbits that appear in the x-y plane to be rosettes of varying eccentricities, and with the morphology dictated by the radial energy. Such regular orbits have constant phase-space coordinates and are considered integrable. These orbits can be represented by a combination of three fundamental frequencies: $\Omega_{\phi}$, $\Omega_{r}$, $\Omega_{z}$ the azimuthal, the radial, and the vertical frequency, respectively. The Jeans equations have been used to perform simplistic assessments of the orbital structure of real galaxies under the assumption that galaxies can be described by the classical integrals of motion, dictated by $E$ and angular momentum $L_{z}$. Unfortunately this condition does not apply to dynamical studies of realistic galaxies. Indeed, commensurate (or resonant) orbits arise when a perturbation, represented by a stellar bar or spiral arms, is applied (Petersen et al. 2019).

For nearly circular orbits the radial frequency is the epicycle frequency, $\Omega_{r}=\kappa(R)$, while $\Omega_{\phi}=\Omega(R)$, the circular motion. Resonances occur when the frequency at which a star encounters the perturbation potential resonates with one of the natural frequencies $\kappa$, $\Omega$ of the stellar orbit. When the stellar orbits are displayed in the frame rotating with the angular velocity of the perturbation, represented by the arms or the bar, then the perturbation will affect only the orbits which close in this frame. When this condition occurs, the bar and arms perturb the star for long time to change the stellar orbit significantly (Contopoulos et al. 1973). Resonances occur when:

[TABLE]

where $l$ is a positive, negative or zero integer and $m=2$ in presence of the bar. In the co-rotation resonance $l=0$ stars move in phase with the bar $\Omega_{\phi}=\Omega_{P}$. First we compute the fundamental frequencies of the stars identified in the velocity plane as Hercules in Figure 2 and verified that their $(\Omega_{\phi}-\Omega_{P})$/$\Omega_{P}<<1$. After 2.5 Gyrs of evolution all stars in the Hercules stream satisfy the condition to be near co-rotation with the stellar bar.

Once we identified the Hercules stream, we traced its position back in time and followed its stars’ orbits in the frame co-rotating with the bar (Figure 3). Stars that at the final time are in the solar vicinity (marked with black points in the right lower panel) were initially distributed following the exponential distribution of disk stars (left top panel).

Note that at 1 Gyrs, the bar is already in place and the majority of Hercules stars circulate the disk, crossing through the edge of the bar (right top panel).

Those stars feel the torques from the pattern potential in their orbits around the Galaxy. Since they are not aligned with the bar these torques add up over several orbital times. Only after long time those stellar orbits will respond to the perturbation and become trojans. Then they are trapped in resonant orbits in co-rotation with the bar but librate around the equilibrium point, L4, which is located 90o away from the bar. Although those orbits feel the torques from the bar pattern potential, they are incapable of aligning with it. Donald Lynden-Bell called them donkey because of their inability to cooperate with a potential well of the perturbation and deepen it by aligning with it. These stars lose and gain angular momentum depending on their location with respect the bar. When close to the bar and moving ahead of the bar, stars are being pulled backwards by the bar and lose angular momentum (Ceverino & Klypin 2007). At the same time, stars that are behind the bar are being pushed forwards gaining angular momentum. However those orbits are never cooperative and will never align with the bar potential (Earn & Lynden-Bell 1996).

The bar torque with a two-fold symmetry is the mechanism that induces the Hercules stars to librate around L4. Therefore, Hercules stars have two-fold symmetric orbits with symmetry of 180o in an axisymmetric potential of a stellar bar. This implies that by symmetry there is an anti-Hercules stream consisting of stars librating around L5, that we name spartans by analogy to moon and minor planets in the solar system. Figure 4 illustrates the situation with the Sun that will encounter in 100 Myrs anti-Hercules, the twin stream located 180o from the solar vicinity.

Note that trojans around L4 cannot librate around L5. This is a difference with previous studies claiming that the Hercules stream consists of stars librating either around L4 or L5 or circulating with the bar (Pérez-Villegas et al. 2017). However, the former study does not account for a time-varying potential, a limitation that might have affected their stellar orbits.

We cannot display the orbits of all Hercules stars in this paper. We shall therefore provide an example of the orbit of a Hercules star in Figures 5-6 followed for 4 Gyrs of evolution in the frame co-rotating with the bar. The star is moving on a retrograde orbit with a period of approximately 108 yrs as illustrated in Figure 5. After 2 Gyrs the star on a periodic orbit gets trapped in a librational motion around L4, with a period of 250 Myrs. Thus, stars in co-rotation resonance with the bar are on retrograde orbits. Their motion can be described by (see Pichon, the Knight prize 1992, Pichon & Lynden-Bell 1993):

[TABLE]

Note that the retrograde periodic orbit is stable even if the orbit is far from L4. However, as the energy increases the amplitude of the orbit around the equilibrium points grows until the star escapes from L4 and continues in its circulation around the center of the Galaxy. As an example to illustrate this point, Figure 6 shows the case of the same star on a periodic orbit in the frame co-rotating with the bar followed for 4 Gyrs of evolution. It becomes a trojan after 2 Gyrs of the disk evolution and at current time librates around L4. In the next 500 Myrs of evolution the star escapes from L4 and moves to a more energetic orbit continuing its circulation across the stellar disk (displayed as the more external retrograde orbit in Figure 6).

As the co-rotation moves outward passing the solar circle we expect the Hercules stream to still be visible at the solar vicinity but its stars will be on orbits with larger guiding radii, hence greater angular momentum than the Sun. Thus, in the kinematic plane they are expected to populate the region with greater Vϕ and will appear as they move inward the Galactic center ($V_{r}\leq 0$) in Figure 2.

The amount of trapped matter near L4 that reaches the solar vicinity after 2.5 Gyrs of evolution is  18%. Other numerical experiments suggest an estimate of  4% (M. Weinberg, private communications) for a barred galaxy. It is unclear at this stage whether differences in the initial conditions lead to different estimate of the trapped matter in numerical simulations.

3.3. Escape from L4

A question concerns the eccentricities of the orbits of Hercules stars. The observational data suggest that Hercules stars have eccentric orbits (Dehnen 2000; Trick et al. 2018). However according to the perturbative theory applied to stellar orbits, when the motion of stars is perturbed the angles can be split in slow and fast variables $\theta_{slow}=\theta_{\phi}-\Omega_{P}t$ is the angle from the Galactic center to the orbit’s epicenter as measured in axes rotating with the bar. The period of an orbit is small compared with its precession rate relative to the rotating perturbation. This assumption allows averaging over the fast motion of the star around its orbit. Formally, averaging over the fast motion of the stars produces an exactly constant fast action: (Lynden-Bell & Kalnajs 1972; Contopoulos et al. 1973; Lynden-Bell 1979):

[TABLE]

Thus, in co-rotation (l=0), $J_{fast}=J_{R}$ is the adiabatic invariant of the motion around that epicenter. Why do trojans have eccentric orbits?

First, we define Hercules stars in the kinematic plane and those are near co-rotation with the bar. Second, we notice that some stars that will be trojans are born already with eccentric orbits (see Figure 7 top panels). Before being trapped in co-rotation resonance with the bar some stars are increasing the random energy of their orbits by the passage of four spiral arms that in our simulation form before the bar. The increase of Hercules eccentricity follows the in-plane heating of the disk that occurs when the non-axisymmetric features develop. In fact the sequence of panels of Figure 7 (top and middle) shows that in the inner part of the disk there is a sharp increase of the radial action of stars that does not significantly change their guiding radius, indicating in-plane disk heating. This occurs when the spiral structure propagates and subsequently the bar starts forming.

Once trojans complete their exchange of angular momentum with the bar, they can produce a very large libration around L4 that allows stars to reach the solar vicinity (see also Pérez-Villegas et al. 2017). As they are on very energetic orbits around L4 they escape from the Lagrangian point. As L4 moves out to larger radii, trojans that are trapped move out to larger guiding radii. In fact, as the bar slows in its pattern speed, the co-rotation region moves outwards and the stars that were trapped in co-rotation migrate radially outwards too.

To illustrate this situation, Figure 7 displays the radial action, $J_{r}$, of the Hercules stream (yellow circles) as a function of their guiding radii (represented by $L_{z}$) as a function of time and as compared to the radial and azimuthal actions of the stars of the disk. At the initial time, the Hercules stars have the same actions as the stars in the background (top left panel). After 2 Gyrs of the evolution of the stellar disk, more and more Hercules stars have guiding radii locked at the co-rotation radius with the stellar bar, located approximately at 5.5 kpc (middle left panel). After another 500 Myrs (at the current time) of evolution the pattern speed of the bar lowered by 5 km s*-1kpc-1* and the co-rotation radius moved out to a galactic radius of 6.0 kpc and Hercules stars that were trapped with the guiding radii at the co-rotation of the bar migrate to larger radii (see Table 2 in the Appendix for the evolution of Resonances locations over time).

Note that the right panels in Figure 7 show the time evolution of the radial action of the stellar disk, displayed against the cylindrical radius with the Hercules stream overlapped and consisting of stars that at the current time have cylindric radii in the solar vicinity (at 8.1 kpc from the Galactic center).

3.4. Limitation of the Model

This numerical experiment suffers from some limitations that we plan to improve in future. First, while the dark halo, the stellar disk and bulge are live, no gas component is included in the Galactic disk. The lack of gas likely affects the evolution of the bar as gas should take part in the angular momentum exchange between the dark halo and the disk. Athanassoula et al. (2013) demonstrated how the gas fraction affects the bar properties. Galaxies with large fractions tend to have shorter and weaker bars. Second, the bulge is described with a Hernquist model and does not rotate. This limitation accelerates the exchange of angular momentum between the bar and the other components of the Galaxy leading to a fast lowering of the bar pattern speed and consequent rapid movement of the resonances outward. The description of the evolution of the trojans hold nevertheless, although the details of the times at capture at L4 and escape will depend on the Galaxy potential adopted.

3.5. Validating tests

The experimantal evidence in support of this model concerns the kinematic plane, in particular the estimate of the Vϕ-Vr at the solar circle at various azimuthal angles that will be possible to validate with Gaia DR3 and DR4. Figure 8 (top row) illustrates the predictions for V${}_{\phi}-$Vr of stars selected in patches of 300 pc each at $\phi=0;-28^{o};-90^{o}$ in the long-bar model. Trojans will de-populate in patches ahead of the Sun along the bar semi-major axis in proximity of the Lagrangian point L2 (left top panel for $\phi=0$). Because trojans will mostly populate the regions around L4, for a patch of stars selected along the bar semi-minor axis (at $\phi=-90^{o}$), more Hercules stars will have low Vϕ in the $V_{\phi}$-Vr plane (right top panel of Figure 8) as compared to the values in the solar vicinity (middle top panel). Note that this behavior seems not to be found in models of a short and fast-rotating bar with the OLR located near the solar radius as displayed in Figure 8 (bottom panels). In this situation Hercules consists of stars changing their orbits from x${}_{1}(1)$ to x${}_{1}(2)$ (Dehnen 2000). Hercules does not increase significantly at $\phi=-90^{o}$. In this case the velocity distributions are evaluated after 1.5 Gyrs of evolution of the disk, when the bar pattern speed is 55 km s*-1kpc-1*, and the bar length is 3 kpc. The OLR is located at 7.8 kpc (see Table 2 in the Appendix). The Sun is placed at 8.6 kpc in order to compare to the fiducial case of Figure 2 of Dehnen (2000) with R${}_{OLR}/R_{\odot}$=0.9.

A quantitative estimate is displayed in Figure 9. Azimuthal, Vϕ, and radial velocity, Vr, distributions are shown for the Hercules stream in the long-bar scenario (left side, dashed lines) and short-bar scenario (right side, dashed lines). The fraction of stars in Hercules, defined as Vϕ=180-220 km s*-1and -50 $\leq$ V${}_{r}\leq$100 km s-1*, is 18% in the solar vicinity, 24% towards the bar semi-minor axis and decreases to 6% for a patch of stars along the bar semi-major axis (see Table 1). While in both scenarios the fraction of stars in Hercules is similar in the solar vicinity and at patches selected along the bar major axis, Hercules is expected to be more prominent along the bar semi-minor axis in the case where the bar is long as compared to the case where it is short and fast rotating.

In the solar vicinity Hercules shows an asymmetry in the $V_{r}$ distribution in both scenarios as displayed in Figure 9. However, in the case of the long bar the asymmetry can be interpreted by having an excess of trojans with V${}_{r}\geq$0. In fact those stars are on donkey orbits in the co-rotating frame of the bar and librate from the center of the Galaxy on a retrograde motion with respect the Sun. Therefore more Hercules stars come from the inner part of the Galaxy as compared to the ones that move towards the Galactic center. The radial velocity distribution in Figure 9 (left side) breaks in the proximity of L4 ($\phi=-90^{o}$) because at that location the same amount of trojans move inwards or outwards.

For a perfectly axisymmetric disk the asymmetry in the $V_{\phi}$-Vr plane will be inverted in the quadrant for $\phi=-90^{o}$ to $\phi=-180^{o}$. Symmetric predictions hold for anti-Hercules formed by stars librating around L5. Finally, Hercules is expected to be less visible at larger Galactic radii, as fewer and fewer trojans will be captured at L4 on such energetic orbits.

4. conclusions

In this paper, we examined the properties of the Hercules stream at the solar vicinity in the long-bar scenario. We carried out GALAKOS, a high-resolution N-body simulation with enough spatial resolution to resolve the solar neighborhood. The simulation consists of the Milky Way model that accounts for a time-varying potential that after 2.5 Gyrs of evolution forms a bar with a length of 4.5 kpc and a pattern speed of 40 km s*-1kpc-1*. While our analysis has much in common with earlier studies of the Hercules stream, it differs in the following respects.

In our N-body simulation the Galactic potential is time-varying. Density waves that characterize the spiral structure and the bar arise from the stars themselves as self-consistent collective disturbances with amplitude that changes by time. We confirm previous studies that Hercules stars are on trapped orbits in co-rotation with the bar (Pérez-Villegas et al. 2017; Monari et al. 2018). However, our results indicate that Hercules is entirely made of trojans, stars that librate only around L4, the stable Lagrangian point located 90 degrees perpendicular to the major axis of the bar. Stars trapped around L5, the other stable Lagrangian point, that also end up as part of Hercules in Pérez-Villegas et al. (2017), belong in our model to anti-Hercules, which is expected to be located at 180o from the Sun. Those stars are not visible in the solar vicinity. 2. 2.

Trojan evolution is dictated by the evolution of the bar and the Galactic potential. In our model they last at capture at L4 only 500 Myrs on very energetic orbits, before escaping without being captured again. 3. 3.

Our theory explains the asymmetry of the radial velocity distribution of Hercules. Because these stars are on donkey orbits, they are retrograde in the frame co-rotating with the bar. At the Sun’s location of -28o of inclination with respect the bar, Hercules will appear as a stream with the majority of stars coming from the inner part of the Galaxy, with positive radial velocity Vr. The asymmetry is expected to break for patches of stars selected towards the bar semi-minor axis (L4). 4. 4.

Finally, we suggest a test in the kinematic plane that allows us to distinguish between Hercules being stars at the outer-Lindblad resonance (OLR) of the bar or trojans in co-rotation with the bar. Trojans will de-populate in patches of stars along the solar circle but align with the major axis of the bar in the proximity of L2 ($\phi=0$). At -90o towards L4 the fraction of Hercules will be 24% of stars of the patch and will have lower azimuthal velocity as comparted to the solar neighborhood. This increased population of Hercules stars towards the semi-minor axis is not expected in models of a short and fast-rotating bar with the OLR located near the solar radius. This anticipation is a valid test for the forthcoming Gaia DR3 and DR4.

E.D.O is grateful to A. Kalnajs and K. Freeman for much generous advise during her visit at the ANU as Stromlo Distinguished Visitor of the Research School of Astronomy and Astrophysics (RSAA). The authors are grateful to C. Pichon, M. Weinberg and J. Navarro for many insightful suggestions. E.D.O. acknowledges support from the Sierra Foundation at the Instituto de Astrofisica de Canarias (IAC) in Tenerife, the Center for Computational Astrophysics at Flatiron Institute, the Institute for Theory and Computation (ITC) at Harvard, and the KITP program for the hospitality during the completition of this work. This research was supported in part at KITP by the National Science Foundation under Grant No. NSF PHY-1748958, and the grant from the Spanish Ministerio de Economía y Competitividad (MINECO) AYA2017-83204-P. GALAKOS is run at the La Palma supercomputer center.

.1. Setting the Initial Conditions for GALAKOS

This section explains our mass models and describes the realization of the initial conditions. We begin with the details of our fiducial simulation. The Milky Way in our study (GALAKOS) consists of a dark matter halo, a rotationally supported disk of stars, and a spherical stellar bulge. The parameters describing each component are independent and models are constructed in a manner similar to the approach described in previous work Springel et al. (2005).

.2. Dark Matter Halo

We model the dark matter mass distribution of each galaxy with a Hernquist profile (Hernquist 1990):

[TABLE]

with cumulative distribution M($<r$)=M${}_{\rm{DM}}r^{2}/(r+a)^{2}$ where $a$ is the radial scale length. This choice for the halo profile is motivated by the fact that in its inner parts, the shape of this profile is identical to the Navarro-Frenk-White fitting formula for the mass density distribution of dark matter halos in cosmological simulations. The total mass of the halo is M${}_{\rm{DM}}=$1x1012 M*⊙* and it is sampled with 60 million particles. Moreover, the radial scale length for the halo is $a=30$ kpc.

.3. Stellar Disk

For the simulation presented here, we adopted an exponential radial stellar disk profile with an isothermal vertical distribution:

[TABLE]

where Md is the disk mass, Rd is the disk scale length, and z0 is the disk scale height. We set Md=4.8x1010M*⊙, Rd*=2.67 kpc, and z0=320 pc throughout in reasonable agreement with the values quoted in the literature for the observed structural parameters of the Milky Way (Bland-Hawthorn & Gerhard 2016). The disk mass is discretized with 24 million particles. The total baryonic mass (stellar disk and bulge) is Md=5.6x1010M*⊙. We do not include a gas component. The rotation curve has a disk fraction $f_{D}=V^{2}_{c}/V^{2}_{\rm{tot}}=0.54$ at R=2.2 Rd*, the radius at which the exponential disk reaches the maximum circular velocity.

disk velocities are chosen by solving the Jeans equations in cylindrical coordinates in the combined disk–halo potential. We set the radial velocity dispersion from the Toomre stability equation:

[TABLE]

where G=1, $\Sigma(R)$ is the stellar surface density, and $\kappa$ the epicyclic frequency. The Q parameter reaches a minimum value of 1.5 between 3 and 7 kpc as displayed in Figure 10 (right panel).

.4. Stellar Bulge

The stellar bulge is described by the Hernquist model (Hernquist 1990). The total mass of the bulge adopted in the simulation is MB=8x109M*⊙* with a scale length $a_{B}$=120 pc.The number of particles sampled in the bulge is 8.4 million. The bulge in this numerical experiment does not rotate.

.5. Action and Angles Estimates

We used the AGAMA software package (Vasiliev 2019) to extract the potential from each snapshot of our simulation and compute the actions and angles. We adopt the Stäckel fudge method (Sanders & Binney 2016) for computing actions and angles under the assumption that the motion is integrable and is locally well-described by a Stäckel potential, separable into prolate spheroidal coordinates. Actions and angles can be computed exactly for a general potential separable in a confocal ellipsoidal coordinate system. AGAMA includes the implementation of this approach. Here we note the essential procedure. A Stäckel potential assumes the form:

[TABLE]

and is defined by two one-dimensional functions, U(u) and V(v), instead of being an arbitrary function of two coordinates. An orbit in such a potential respects three integrals of motion: energy (E), angular momentum (Lz), and a non-classical third integral. Given the three integrals, the equations defining the turn-around points (min/max values of u, v) can be numerically solved and then the actions, angles, and frequencies can be computed by one-dimensional numerical integration. The Stäckel fudge method uses the same expressions, but substitutes the real potential instead, which is not of a Stäckel form. Actions computed in this way are approximate, i.e., not conserved along a numerically integrated orbit. The accuracy of approximation depends on the only free parameter in the method: the focal distance. This parameter is obtained from the combination of second derivatives of the potential; hence this equation is implemented in the real potential at each input point.

.6. Spectral Analysis of the Stellar Disk

The Galactic bar changes the internal structure of disk galaxies by rearranging the angular momenta and energy of star orbits that would otherwise be conserved. After 2.5 Gyrs of evolution, the bar is in place with a length of 4.5 kpc. Figure 11 displays the power spectra for $m=2$ harmonic mode, measured from snapshots sampled every 5 Myrs for the stellar disk in the 2,000$\leq$ t/Myr $\leq$ 2,500 interval. The solid line shows $m$ times the circular angular frequency $\Omega=V_{c}/R$, which determines the co-rotation radius of the bar located around 6 kpc. The ratio of co-rotation to the bar length in this model is 1.3. This value is similar to the estimates of other barred galaxies and indicates that the MW bar is fast-rotating Aguerri et al. (2015). At t=1.5 Gyrs the bar has a length of 3 kpc and a pattern speed of 55 km s*-1kpc-1*, while at 2.5 Gyrs its pattern speed is $\sim$40 km s*-1kpc-1*.