When Cold Radial Migration is Hot: Constraints from Resonant Overlap

A particular family of orbits, sometimes called "horseshoe" orbits, underlie the physics that can lead to cold torquing (Goldreich & Tremaine 1982; Sellwood & Binney 2002; Daniel & Wyse 2015). These orbits occur near the corotation resonance, where the orbital frequency of stars, Ω(R), equals the pattern speed, Ω_p, of a non-axisymmetry in the disk. A star with a trajectory that belongs to this orbital family, hereafter called "'trapped," will have periodic changes in orbital angular momentum, L_z, causing its mean orbital radius, R_L, to oscillate about the radius of corotation, R_CR. The mean orbital radius of a star can be defined by using its orbital angular momentum such that

$$\begin{eqnarray}&&{R}_{L}=\displaystyle \frac{{v}_{\phi }}{{v}_{c}}R,\end{eqnarray} \tag{ 2 }$$

where v_c is the orbital circular velocity at R_L, v_ϕ is the instantaneous tangential velocity about the disk, and R is the instantaneous radial position of the star. Critical to the definition for cold torquing, there is little to no change in the star's orbital circularity after a trapped orbital period. Should the non-axisymmetric pattern be transient, a trapped star could have a permanent change in its orbital angular momentum, and thus mean orbital radius, and no change in its orbital circularity (Sellwood & Binney 2002).

An approximation for the behavior of a star in a trapped orbit was derived in Section 3.3.3 of Binney & Tremaine (1987) and invoked to describe limits on cold torquing by Sellwood & Binney (2002). This approximation assumes a smooth disk that is perturbed by a weak bar pattern with potential amplitude, $| {{\rm{\Phi }}}_{b}|$, at the radius of corotation. By making the further assumption that the underlying disk has a flat rotation curve, the maximum radial excursion for a star in a trapped orbit scales as ${\rm{\Delta }}R\propto \sqrt{| {{\rm{\Phi }}}_{b}| }$, and the minimum timescale for the smallest excursions scale as ${T}_{\min }\propto {R}_{\mathrm{CR}}/\sqrt{| {{\rm{\Phi }}}_{b}| }$. While these provide guidelines for approximating the efficiency of cold torquing, it is informative to refine these scaling relations.

Daniel & Wyse (2015) demonstrated that the location of corotation does not occur only at the radius where Ω(R) = Ω_p. Rather corotation can be better described in coordinate space by a 2D region in the plane of the disk. Disk stars with their mean orbital radius within this "corotation region" are, to first order, trapped in stable orbits around corotation. The analytic criterion is derived in action space where vertical action is assumed to be separable in a thin disk. A transformation can be made from 4D action-angle space to 4D phase space for a given radially local rotation curve, spiral strength, and pitch angle. A star is trapped when its mean orbital radius, R_L, is approximately within the 2D coordinate-space corotation region, where there is a higher order dependence on orbital circularity. The shape of the corotation region depends on the morphology of the spiral pattern in that the radial range of the corotation region increases with increasing spiral strength and openness of the spiral arms. There is also a dependence on the rate of divergence with radial distance from corotation between the spiral arm's pattern speed and the orbital frequency for disk stars. For example, spirals that corotate with the disk at all radii (Ω(R) = Ω_p(R)) have the largest corotation regions, possibly spanning the full radius of the disk (Grand & Kawata 2012; D'Onghia et al. 2013), while a spiral with radially constant pattern speed (Ω_p = const) in a Keplerian disk (Ω(R) ∝ R^−3/2) will have a much less extended corotation region (see Equation (32) in Daniel & Wyse (2015), for their generalized analytic expression for the shape of the corotation region).

The inner and outer Lindblad resonances (ILR/OLRs) are where a disk star passes or is passed by the spiral pattern at the star's epicyclic frequency, κ. The LRs and their harmonics occur at radii where

$$\begin{eqnarray}&&\kappa =\pm m(n+1)\,[{{\rm{\Omega }}}_{p}-{\rm{\Omega }}(R)]\end{eqnarray} \tag{ 3 }$$

is satisfied, where n stands for the nth harmonic of the LR and m is the number of spiral arms. The ILR/OLRs have harmonic number n = 0 in this notation. In a disk with a flat rotation curve, the radial locations for the LRs are given by

$$\begin{eqnarray}&&{R}_{\mathrm{LR}}^{(n+1)}=\left(1\mp \displaystyle \frac{\sqrt{2}}{m(n+1)}\right)\displaystyle \frac{{v}_{c}}{{{\rm{\Omega }}}_{p}},\end{eqnarray} \tag{ 4 }$$

where ${{\rm{\Omega }}}_{p}={v}_{c}/{R}_{\mathrm{CR}}$.

Higher order harmonics of the LRs (n > 0) are alone not expected to have a significant impact on disk kinematics. However, since the corotation region has some finite radial range, it is possible for stars in otherwise stable, trapped orbits that have radially periodic orbits about the corotation radius to temporarily also be in resonance with a harmonic of the LRs. Chirikov (1979) predicted that stochastic⁷ behavior emerges in cases where more than one resonance occupies the same N-dimensional space. He drew particular attention to the case of a pendulum under the influence of an external, periodic perturbation, a well-known example of chaotic behavior. Explorations into the kinematic response of a stellar disk in the presence of multiple, non-axisymmetric patterns with different pattern speeds suggest there is a stochastic response for stars at radii that meet a resonant criterion for each pattern (e.g., Quillen 2003; Jalali 2008; Minchev & Famaey 2010).

The nature and degree of the dynamical response to resonant overlap can be challenging to identify. "Wild" (Martinet 1974) or ergotic (Athanassoula et al. 1983) behavior associated with resonant overlap can be identified in regions of a surface of section (SoS) diagram by irregularly distributed consequents (see Binney & Tremaine 2008, Section 3.7.3). Pichardo et al. (2003) combined SoS and Lyapunov analyses to show that for a sufficiently strong perturbation, orbits in these regions exhibit chaotic behavior. SoS analysis can have limited utility since, in order to to resolve a chaotic signature, it necessitates a significant number of orbits in order to populate the phase space within a given energy slice, and these are best evolved over several dynamical periods.

In a simultaneous study to the present paper, D. Schaffner et al. (2019, in preparation) use PESC analysis to effectively investigate the nature of the dynamical response in cases where the corotation region overlaps an ultraharmonic produced by the same perturbation. Since perturbation theory cannot predict nonlinear effects from resonant overlap, PESC analysis provides an avenue to identify such a dynamical response. The PESC method is well suited to identifying chaos in poorly resolved or limited data sets since it does not need to fill the phase space, it best identifies a chaotic response on shorter, rather than longer, timescales, and it uses only one governing dimension. In the case of resonant overlap, timescales from classical perturbation theory do not apply since nonlinear effects causing the chaotic response happen on arbitrarily short timescales. For some of the same models described in this work, D. Schaffner et al. (2019, in preparation) find a significant chaotic signature using only orbital angular momentum for star particles in small batches of a few ×10² orbits on timescales less than an orbital period and time resolution equal to ∼T_dyn/20. Figure 1 shows an example of the results using PESC analysis for orbits that are resonant with only corotation (squares, Type 1) and orbits that are resonant with both corotation and an ultraharmonic (circles, Type $3\to 2$) in model M6e (described below in Section 3). There are two relevant points. First, there can be a chaotic dynamical response even for a single spiral pattern. Second, that chaotic response is readily quantifiable.

Zoom In Zoom Out Reset image size

The nature of the dynamical response during resonant overlap is not the aim of the present study. Rather, the focus is on whether or not the disk response due to resonant overlap can place an upper limit on the efficiency of cold torquing. In the case where the corotation region overlaps the first harmonic of the LR any stochastic or chaotic dynamical response is of interest. The timescale for a response at the ultraharmonic is reasonable since our tracer particle simulations show that a trapped orbital period is on order a few epicyclic periods. Such a response minimally implies that trapped stars with excursions in their mean orbital radii ($| {\rm{\Delta }}{R}_{L}|$) greater than the distance between the n = 1 harmonic of the LRs

$$\begin{eqnarray}&&{\rm{\Delta }}{R}_{\mathrm{LR}}^{(2)}=\displaystyle \frac{\sqrt{2}{R}_{\mathrm{CR}}}{m},\end{eqnarray} \tag{ 5 }$$

would have a chaotic response. A strong dynamical response at resonant overlap could lead to a change in orbital family from trapped orbits to non-trapped orbital families. A chaotic response is also a likely driver for radial migration that is associated with kinematic heating. Since kinematic heating is not expected to be associated with radial migration around corotation, such resonant overlap is of significant interest.

Initial conditions for each star particle are produced by sampling the distribution function, f_new, from Dehnen (1999). This distribution function resembles a kinematically warm 2D disk and has moments that are similar to observed trends in the Milky Way, namely, a flat rotation curve (e.g., Rubin 1983; Sofue et al. 2009, and references therein) and an exponential surface density profile (e.g., Freeman 1970; van der Kruit 1987; Jurić et al. 2008; de Jong et al. 2010). In energy-momentum space this distribution function is given by

$$\begin{eqnarray}&&{f}_{\mathrm{new}}(E,{L}_{z})=\displaystyle \frac{{\rm{\Sigma }}({R}_{E})}{\sqrt{2}\pi {\sigma }_{R}^{2}({R}_{E})}\exp \left\{\displaystyle \frac{{\rm{\Omega }}({R}_{E})[{L}_{z}-{L}_{c}({R}_{E})]}{{\sigma }_{R}^{2}({R}_{E})}\right\},\end{eqnarray} \tag{ 6 }$$

where R_E is the orbital radius for a star in a circular orbit with energy E, Ω(R) is the circular frequency at a given radial coordinate, L_z represents orbital angular momentum about the z-axis, L_c(R) is the orbital angular momentum for a star in a circular orbit at radius R, and σ_R(R) is the radial velocity dispersion at radius R.

The Python-based galaxy modeling package galpy v1.4 (Bovy 2015)⁸ is used to produce a set of initial conditions for each model. The galpy command dehnendf is set to numerically refine the initial distribution function given by Equation (6) over 20 iterations so that the final distribution function better approximates an exponential surface density profile and flat rotation curve (see Bovy 2015, for more details). This is done by setting β = 0 and niter=20. galpy uses natural length unit h₀ and natural velocity unit is v₀. Both the radial velocity dispersion profile (σ_R(R)) and the radial surface density profile, Σ(R), are assumed to be exponential. The scale length, h_σ, for σ_R(R) is set to h_σ = h₀ and is assumed to be three times the scale length, h_R, for Σ(R). The radial velocity dispersion profile is normalized so that ${\sigma }_{R}({h}_{0})=0.16\,{v}_{0}$.

The galpy function sample is used to produce initial conditions for a set of test particles in several models. The number of initial conditions, N, minimum and maximum sampling radii, R_min and R_max respectively, for each model are given in Table 1. The radial ranges for the sampling annuli are chosen to span approximately 4 kpc centered on corotation in order to provide a complete sampling of the corotation region with minimal inclusion of initial conditions for stars that are never trapped at the corotation resonance.

4D phase-space coordinates were recovered for each initial condition by setting parameter returnOrbit=True. Each phase-space coordinate is transformed from natural to physical coordinates using v₀ = 220 km s⁻¹ and h₀ = 8 kpc.

Each orbital trajectory is calculated from the initial phase-space coordinate using the second order leapfrog orbital integrator with fixed δt = 10⁵ yr time steps initially described in Daniel & Wyse (2015). Orbits are evolved through the potential given by

$$\begin{eqnarray}&&{\rm{\Phi }}(R,\phi ,t)={{\rm{\Phi }}}_{0}(R)+{{\rm{\Phi }}}_{1}(R,\phi ,t),\end{eqnarray} \tag{ 7 }$$

where Φ₀(R) is the radially dependent axisymmetric disk potential, and Φ₁(R, ϕ, t) is a time-dependent non-axisymmetric spiral perturbation to the underlying potential.

An analytic form for the underlying potential is chosen such that it is approximately consistent with the moments produced by the adopted distribution function (f_new, Equation (6)) for a 2D disk, but it is not strictly self-consistent through Poisson's equation and the collisionless Boltzmann equation. The underlying disk potential is set to

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{0}(R)={v}_{c}^{2}\mathrm{ln}(R/{R}_{p}),\end{eqnarray} \tag{ 8 }$$

where the circular velocity is set so that v_c = v₀ and the scale length for the potential is R_p = 1 kpc, in order to reproduce a flat rotation curve in two dimensions. For the sake of analytic simplicity the spiral perturbation to the potential is assumed to be a density wave given by (Lin et al. 1969),

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{1}(R,\phi ,t)={{\rm{\Phi }}}_{s}(R,t)\cos \left[\alpha \mathrm{ln}(R/{R}_{\mathrm{CR}})+m({{\rm{\Omega }}}_{p}t-\phi )\right],\end{eqnarray} \tag{ 9 }$$

where ϕ is the azimuthal coordinate, and the constant $\alpha =m\cot \theta$ depends on the spiral pitch angle, θ.

The amplitude of the spiral perturbation,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{s}(R,t)=\displaystyle \frac{2\pi {GR}{\rm{\Sigma }}(R)\epsilon (t)}{\alpha },\end{eqnarray} \tag{ 10 }$$

depends on the radial surface density profile and the time-dependent fractional amplitude of the surface density, (t).

The adopted time dependence for (t) has the Gaussian form

$$\begin{eqnarray}&&\epsilon (t)={\epsilon }_{0}\,{e}^{-{\left(t-2{T}_{\mathrm{Dyn}}\right)}^{2}/2{\sigma }_{t}^{2}},\end{eqnarray} \tag{ 11 }$$

so that the peak spiral amplitude occurs two orbital periods, 2T_dyn, after the simulation begins, where the standard deviation is set to ${\sigma }_{t}={T}_{\mathrm{Dyn}}$. These assumptions ensure slow growth and decay for the perturbation, thus avoiding a non-adiabatic dynamical response. The value for the maximum fractional amplitude is set to ₀ = 0.3.

This study explores two forms for the radial surface density profile for the disk, Σ(R). The distribution function, f_new, produces an exponential surface density profile in a 2D disk, which is expressed analytically by

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)={{\rm{\Sigma }}}_{0}\,{e}^{-R/{R}_{d}},\end{eqnarray} \tag{ 12 }$$

where the disk scale length for the surface density R_d is set to match the prescription used for the distribution function ${R}_{d}={h}_{R}={R}_{0}/3$.

An alternate set of models uses an inverse radial form for the surface density

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)=\displaystyle \frac{{{\rm{\Sigma }}}_{0}}{R}.\end{eqnarray} \tag{ 13 }$$

This form is self-consistent for the adopted axisymmetric potential, Φ₀(R) (Equation (8)), but not with the assumed distribution function (Equation (6)). In all models Σ₀ is set so that the the surface density at R₀ is 50 M_⊙ pc⁻², similar to conditions in the solar neighborhood of the Galaxy (e.g., Kuijken & Gilmore 1991).

Test particle trajectories in each model are evolved for four orbital periods, where T_dyn is given in Table 1.

Each model is named to signify its assumed radius of corotation, ${R}_{\mathrm{CR}}=\{4,5,6,7,8,9,10\}$ kpc, and adopted surface density profile. For example, model M6e assumes R_CR = 6 kpc and an exponential, "e," surface density profile.

Stars in the plane of a non-axisymmetric potential have an invariant energy defined by the Jacobi integral

$$\begin{eqnarray}&&{E}_{J}=E-{{\rm{\Omega }}}_{p}{L}_{z},\end{eqnarray} \tag{ 14 }$$

where the orbital energy, E, and angular momentum about the vertical axis, L_z, are taken to be in the non-rotating frame. The orbital energy can be defined as

$$\begin{eqnarray}&&E={E}_{c}({R}_{L})+{E}_{\mathrm{ran}},\end{eqnarray} \tag{ 15 }$$

where E_c(R) is the orbital energy for a star in a circular orbit with radius R in the unperturbed potential, R_L is the radial coordinate in the unperturbed potential for a star in a circular orbit with unit mass angular momentum L_z = R_L v_c, and E_ran is the energy associated with noncircular motions, hereafter called "'random energy."

Trapped orbits are stable orbits that are resonant with corotation (see discussion in Section 2). A star with zero random energy (E_ran = 0) is trapped when the parameter

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{c}=\displaystyle \frac{{E}_{J}-{h}_{\mathrm{CR}}}{| {{\rm{\Phi }}}_{s}{| }_{\mathrm{CR}}}\end{eqnarray} \tag{ 16 }$$

has absolute value equal to or less than unity (−1 ≤ Λ_c ≤ 1) (Contopoulos 1978), where the invariant h_CR is the Jacobi integral for a star in a circular orbit at the radius of corotation (R_CR) in the unperturbed underlying axisymmetric potential. For a star with some finite noncircular energy (E_ran) to be trapped in a resonant orbit around corotation, the invariant Λ_c must be replaced with the time-dependent quantity, Λ_nc(t) (Daniel & Wyse 2015, their Equation (22)). For a rotationally supported disk, v_ran/v_c ≲ 0.2 where v_ran is the velocity associated with E_ran, with a flat rotation curve, the value of Λ_nc(t) can be expressed as (Daniel & Wyse 2015, their Equation (33))

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\mathrm{nc}}(t)={{\rm{\Lambda }}}_{c}-\left(\displaystyle \frac{{R}_{L}(t)}{{R}_{\mathrm{CR}}}\right)\left(\displaystyle \frac{{E}_{\mathrm{ran}}(t)}{| {{\rm{\Phi }}}_{s}{| }_{\mathrm{CR}}}\right),\end{eqnarray} \tag{ 17 }$$

and the criterion for resonant trapping around corotation is satisfied when

$$\begin{eqnarray}&&-1\leqslant {{\rm{\Lambda }}}_{\mathrm{nc}}(t)\leqslant 1.\end{eqnarray} \tag{ 18 }$$

The expression in Equation (17) assumes the epicyclic approximation in order to convert from action space to energy–angular momentum space.

Figure 2 shows an example of an orbit that meets the criterion for a stable orbit trapped at the corotation resonance. The left panel shows the trajectory (dashed, red) of the star particle evolved over 0.5 Gyr in a spiral potential with R_CR = 6 kpc. The mean orbital radius (R_L, Equation (4)) is indicated by the solid (black) curve and is in the corotation region (shaded) throughout the simulation. The radial coordinate (dotted, red), guiding center radius (solid, black), their relative distances from the LRs (dashed, green), and ultraharmonics (dotted, green) are also shown in the bottom of the right-hand panel. Since this is a trapped orbit the mean orbital radius radially oscillates about the corotation radius. The top right-hand panel of Figure 2 shows the value for Λ_nc(t) (Equation (17)), which satisfies the corotation criterion (Equation (18)) at all times. The middle right-hand panel is the square root of the absolute value for the change in energy associated with noncircular motions ($| {\rm{\Delta }}{E}_{\mathrm{ran}}{| }^{1/2}$), which varies throughout the simulation, yet the orbit remains trapped with the corotation resonance.

Zoom In Zoom Out Reset image size

Figures 3 and 4 show contours for the corotation region in models that use the exponential surface density profile in Equation (12) (M4e, M6e, M8e, and M10e) and the inverse radial surface density profile in Equation (13) (M4i, M6i, M8i, and M10i) for four different corotation radii. These models were selected in order to illustrate how the distribution of resonances in the disk depends on the radius of corotation and assumed surface density profile when all other parameters are equal. The corotation regions are shaded and outlined by a solid, thick (black) curve, where minima in the spiral potential are shown as short-dashed (magenta) curves. LRs are indicated by long-dashed (green) curves, ultraharmonics by dotted (green) curves, and radii of corotation by a solid (green) curve. With decreasing corotation radius, the distances between the LRs and ultraharmonics decrease, while the area of the corotation region depends on the assumed surface density profile. For models assuming an exponential disk, the corotation region is broader toward the galactic center and therefore trends toward a greater degree of resonant overlap for spiral patterns with a smaller corotation radius (and higher pattern speed). Models adopting the inverse radial surface density profile have radially independent spiral strength and therefore the degree of resonant overlap depends only on the fractional amplitude for the spiral pattern for constant number of arms and pitch angle.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The expectation for a star trapped at the corotation resonance is that its mean orbital radius (R_L) will oscillate about corotation (R_CR) indefinitely with negligible change to its circularity. The primary focus of this work is to critically examine that expectation in cases when a trapped star simultaneously meets another resonant criterion.

Each test particle's trajectory is used to calculate a time series array for its orbital angular momentum (L_z), the associated mean orbital radius (R_L), the random orbital energy (E_ran), and Λ_nc(t). These arrays give a quantitative measure for whether or not a star instantaneously meets any dynamical resonant criteria. A particle is trapped in a resonant orbit about corotation when Equation (18) is satisfied via the instantaneous value for Λ_nc(t) (Equation (17)). It is resonant with a Lindblad or ultraharmonic resonance when R_L(t) = R_LR⁽ⁿ⁺¹⁾ (Equation (4)).

Each of these quantities is calculated using the time-dependent spiral potential described in Section 3.2, with the exception of the value for Λ_nc(t). The value for Λ_nc(t) is calculated using a time-independent spiral amplitude set to equal the time-dependent potential's maximum amplitude, ${{\rm{\Phi }}}_{s}(R)={{\rm{\Phi }}}_{s}(R,2{T}_{{\rm{d}}\mathrm{yn}})$. This choice was made for the following reason. In a potential that includes a transient spiral pattern, all stars begin and end in orbits that are not trapped at corotation since corotation only exists in the presence of a non-axisymmetric perturbation. Using constant Φ_s(R) to evaluate Λ_nc(t) identifies stars that would remain trapped for the duration of the simulation when the spiral amplitude is large enough for them to be trapped. By adopting this definition, it is possible to distinguish between star particles that are no longer trapped after an episode of resonant overlap from those that are temporarily trapped in Type 1 orbits due to spiral amplitude growth and decay. A comparative study showed that by using a time-dependent spiral potential to calculate Λ_nc,2(t) there was a small degree of contamination of Type 4 orbits with visually identified Type 1 orbits.

Each trajectory is categorized into one of five types. Star particles that remain in a stable, resonant orbit about corotation for the duration of the simulation are categorized as Type 1. Type 2 orbits are those that never satisfy this criterion and are therefore never in trapped orbits, circulating about the galactic center in the frame rotating with the spiral pattern for the duration of the simulation. Type 3 orbits begin in trapped orbits and have at least one episode when the mean orbital radius crosses an ultraharmonic, thus simultaneously meeting two dynamical resonant criteria. Type 3 orbits are subdivided into Type $3\to$ 1, which remain in trapped orbits after crossing an ultraharmonic, and Type $3\to 2$, which are not in trapped orbits at the end of the simulation after crossing an ultraharmonic. All other cases are categorized as "Other." In no case does a trapped orbit have a mean orbital radius that crosses an LR. A summary of these categories is given in Table 2.

Table 2.  Orbital Categorization Scheme

Category	Description
Type 1	Trapped in resonant orbit about corotation
Type 2	Never trapped at corotation
Type $3\to 1$	In a trapped orbit before and after crossing
	an ultraharmonic
Type $3\to 2$	In a trapped before crossing an ultraharmonic,
	but not after
Other	Any other scenario

Download table as:  ASCIITypeset image

Figure 5 shows an illustrative example of a Type $3\to 2$ orbit, where the star particle does not remain in a trapped orbit after crossing the outer ultraharmonic. In this example, the orbital trajectory is evolved over 0.5 Gyr and begins in a trapped orbit. At t ∼ 0.25 Gyr, the star particle is in resonance with both corotation ($1\leqslant | {{\rm{\Lambda }}}_{\mathrm{nc},2}|$) and the outer ultraharmonic (${R}_{L}={R}_{\mathrm{LR}}^{(2)}$), resulting in an increase in noncircular orbital energy (E_ran). This kinematic heating causes the orbit to no longer be resonant with corotation (Λ_nc,2 < −1), and the star particle is therefore able to cross a spiral arm. The star is briefly trapped starting at t ∼ 0.4 Gyr, but ends in an orbit that is resonant with the outer ultraharmonic only. The orbital response to resonant overlap, being chaotic, is highly variable for very similar initial conditions. In practice, this means that Type 3 orbits commonly have intervals when the orbital trajectory oscillates between being trapped and not trapped with the corotation resonance before being permanently non-trapped. As such, it is likely that the distinction between Type $3\to 1$ and Type $3\to 2$ is somewhat unresolved, as the determination is time-dependent, and given a long enough timescale, many Type $3\to 1$ orbits eventually evolve into Type $3\to 2$ orbits. However, the timescale for the transient spiral arms in the current simulations (σ_t = T_dyn) is several times less than the timescale for the simulations (T = 4 T_dyn), and so the classifications presented here are expected to be qualitatively robust. We have not accounted for changes to the spiral perturbation's pattern speed. Analyses of the pattern speeds in simulations with transient spiral arms typically show a time dependence for spiral pattern speed (e.g., Roškar et al. 2012; D'Onghia et al. 2016). Further exploration is required in order to understand how an evolving pattern speed would affect trapped orbits.

Zoom In Zoom Out Reset image size

The top panels in Tables 3 and 4 summarize the distribution of orbits among each of the orbital categories for each model.

Table 3.  Distribution of Orbits for Exponential Surface Density Profile

	M4e	M5e	M6e	M7e	M8e	M9e	M10e
Type 1	842 (3.4%)	862 (5.8%)	3910 (7.8%)	2616 (17.4%)	4782 (19.1%)	4999 (33.3%)	8799 (35.2%)
Type 2	12,761 (51.0%)	7293 (48.6%)	29,179 (58.4%)	6631 (44.2%)	14,306 (57.2%)	6861(45.7%)	13,080 (52.3%)
Type $3\to 1$	3027 (12.1%)	1794 (12.0%)	4547 (9.1%)	1169 (7.8%)	952 (3.8%)	463 (3.1%)	314 (1.3%)
Type $3\to 2$	4462 (17.9%)	2733 (18.2%)	6357 (12.7%)	2322 (15.5%)	2110 (8.4%)	740 (4.9%)	265 (1.1%)
Other	3908 (15.6%)	2318 (15.5%)	6007 (12.0%)	2262 (15.1%)	2849 (11.4%)	1937 (12.9%)	2542 (10.2%)
Type 3	7489 (30.0%)	4527 (30.2%)	10,904 (21.8%)	3491 (23.3%)	3062 (12.3%)	1203 (8.0%)	579 (2.3%)
Type 1 + 3	8331 (33.3%)	5389 (35.9%)	14,814 (29.6%)	6107 (40.7%)	7844 (31.4%)	6202 (41.4%)	9378 (37.5%)
Type $3\to 2$/Type 3	59.6%	60.4%	58.3%	66.5%	68.9%	61.5%	45.8%
Type 3/Type 1 + 3	89.9%	84.0%	73.6%	57.2%	39.0%	19.4%	6.2%
Type $3\to 2$/Type 1 + 3	53.6%	50.7%	42.9%	38.0%	26.9%	11.9%	2.8%

Note. Number and percentage of trajectories that satisfy the orbital categorization scheme outlined in Table 2 for models that use an exponential surface density profile (Σ(R) ∝ e^−R).

Download table as:  ASCIITypeset image

Table 4.  Distribution of Orbits for Inverse Radial Surface Density Profile

	M4i	M6i	M8i	M10i
Type 1	4598 (18.4%)	8968 (17.9%)	5199 (20.8%)	5062 (20.3%)
Type 2	12,596 (50.4%)	29,074 (58.2%)	14,415 (57.7%)	13,195 (52.8%)
Type $3\to 1$	1245 (5.0%)	2447 (4.9%)	742 (3.0%)	508 (2.0%)
Type $3\to 2$	782 (3.1%)	2097 (4.2%)	1758 (7.0%)	1934 (7.7%)
Other	5779 (23.1%)	7414 (14.8%)	2886 (11.5%)	4301 (17.2%)
Type 3	2027 (8.1%)	4544 (9.1%)	2500 (10.0%)	2442 (9.8%)
Type 1 + 3	6625 (26.5%)	13,512 (27.0%)	7699 (30.8%)	7504 (30.0%)
Type $3\to 2$/Type 3	38.6%	46.2%	70.3%	79.2%
Type 3/Type 1 + 3	30.6%	33.6%	32.5%	32.5%
Type $3\to 2$/Type 1 + 3	11.8%	15.5%	22.8%	25.8%

Note. Number and percentage of trajectories that satisfy the orbital categorization scheme outlined in Table 2 for models that use the inverse radial surface density profile (Σ(R) ∝ R⁻¹).

Download table as:  ASCIITypeset image

The middle panel in each table gives the total number and percentage of orbits that experience resonant overlap (Type 3) in each simulation as well as the number and percentage of all orbits that have initial conditions for trapped orbits (Type 1 + 3). The trend in the initial fraction of trapped orbits roughly follows the ratio between the area of the corotation region and the area of the annulus of sampled initial conditions for each model. (Also see Figures 11 and 14(d) for Model W—and relevant discussion—in Daniel & Wyse 2018, which show the maximum width for the corotation region and fraction of stars initially in trapped orbits for model parameters equivalent to those in Table 3).

The bottom panel in each table quantifies several scaling relations that are illustrated in Figure 6. In all models, between 40%–80% of trapped stars have an episode of resonant overlap and are not in trapped orbits at the end of the simulation (Type $3\to 2$/Type 3) as indicated by a dotted (black) line in Figure 6. This suggests that the dynamical response of trapped stars to resonant overlap alters a significant number of orbital classifications from those that lead to well-behaved radial migration due to cold torquing to orbits that inhabit other regions of the phase space and could be kinematically heated. The fraction of the orbits that are initially trapped (Type 3 + 1) and experience resonant overlap (Type 3) during the total time evolved (4T_dyn) is indicated by a solid (red) line (Type 3/Type 3 + 1). For the series of models that uses an inverse radial surface density profile (Σ ∝ R⁻¹), the radial size of the corotation region closely scales with the distance between the ultraharmonics (${\rm{\Delta }}{R}_{\mathrm{LR}}^{(n+1)}$), so the degree of resonant overlap is nearly constant for these models (∼30%). For the series of models that uses an exponential surface density profile (Σ ∝ e^−R) the degree of resonant overlap increases toward the galactic center from ∼6% for M10e to ∼90% for M4e. The shaded (salmon) region shows the ratio of the number of stars initially in trapped orbits that experience resonant overlap and end the simulation in non-trapped orbits to the number of stars that begin in trapped orbits (Type $3\to 2$/Type 1 + 3). These trends are similar to those for Type 3/Type 1 + 3 since the fraction for Type $3\to 2$/Type 3 is nearly constant.

Zoom In Zoom Out Reset image size

Figures 7 and 8 show the time evolution for the changes in the orbital angular momentum and random orbital energy for populations of star particles in three classification categories. Figure 7 shows models using an exponential surface density profile to determine the spiral amplitude, and Figure 8 uses an inverse radial profile. Orbital classifications shown are Type 1 (always trapped—solid, black), Type $3\to 1$ (resonant overlap occurs, but remain trapped—dashed, orange), and Type $3\to 2$ (resonant overlap occurs and no longer trapped—dotted, teal). Top panels show rms changes in orbital angular momentum, ${\left\langle {\left({\rm{\Delta }}{L}_{z}\right)}^{2}\right\rangle }^{1/2}$, using mean orbital radius, R_L, as a proxy. Horizontal lines indicate half the distance between ultraharmonics as expressed in Equation (5). Bottom panels quantify kinematic heating as the square root of the absolute value of the sum of the changes in random orbital energy, $| {\rm{\Delta }}{E}_{\mathrm{ran}}{| }^{1/2}$, so that these changes are expressed in units of velocity. Time is expressed in units of orbital periods, T_dyn, where the vertical gray line indicates the moment of peak spiral amplitude.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Trapped orbits (Type 1) have an increase in ${\left\langle {\left({\rm{\Delta }}{L}_{z}\right)}^{2}\right\rangle }^{1/2}$ over time and with increasing spiral amplitude in each model. In most cases, the maximum value for ${\left\langle {\left({\rm{\Delta }}{L}_{z}\right)}^{2}\right\rangle }^{1/2}$ occurs nearly concurrently with the peak spiral amplitude (vertical, gray line). One exception is the offset of this peak for Type 1 (always trapped) orbits, which likely reflects a timescale difference between the imposed spiral lifetime governed by ${\sigma }_{t}$ (Equation (11)) in these simulations and the timescale for self-gravitational transient spiral structure. Cold torquing likely plays a role in the self-regulation of transient spiral amplitude (Sellwood & Binney 2002). Further, the timescale for the radial oscillations of a trapped orbit depends on the maximum radial excursion for that orbit. For a population of trapped stars, the peak in ${\left\langle {\left({\rm{\Delta }}{L}_{z}\right)}^{2}\right\rangle }^{1/2}$ likely occurs on the timescale that maximizes the phase mixing of trapped orbits, where that timescale changes with the artificially imposed spiral growth. Relevant to this study is the amplitude of the peak in ${\left\langle {\left({\rm{\Delta }}{L}_{z}\right)}^{2}\right\rangle }^{1/2}$, which is robust whether or not the spiral lifetime is artificially imposed.

The degree of kinematic heating for trapped orbits is not negligible, but is relatively insignificant when compared to other orbital classifications. The low grade heating of Type 1 orbits likely arises from a combination of factors, including interactions with spiral arms away from corotation and near passes with the ultraharmonics without crossing. Type 1 orbits experience the greatest degree of heating in models with corotation closer to the inner disk when the surface density profile is exponential. This same demographic describes models with a greater degree of resonant overlap, and therefore, correlates with the fraction of Type 1 orbits that have near encounters with an ultraharmonic.

Trapped orbits that experience resonant overlap (Type 3) also have a rise in ${\left\langle {\left({\rm{\Delta }}{L}_{z}\right)}^{2}\right\rangle }^{1/2}$, but these changes are associated with kinematic heating. In all cases, orbits that are not trapped after experiencing resonant overlap (Type $3\to 2$) have larger changes in both ${\left\langle {\left({\rm{\Delta }}{L}_{z}\right)}^{2}\right\rangle }^{1/2}$ and $| {\rm{\Delta }}{E}_{\mathrm{ran}}{| }^{1/2}$ than those orbits that remain in trapped orbits (Type $3\,\to$ 1). Nonetheless, trapped orbits that experience resonant overlap are kinematically heated regardless of whether or not they remain in trapped orbits. With the exception of M10e, Type 3 orbits have larger ${\left\langle {\left({\rm{\Delta }}{L}_{z}\right)}^{2}\right\rangle }^{1/2}$ by the end of the simulation than their Type 1 counterparts. The outstanding case of M10e is likely to do with the very small fraction of Type 3 orbits (2.3%).

The large fraction of orbits that are not trapped after an episode of resonant overlap (Type $3\to 2$/Type 3) suggests that this dynamical response could be important to understand how disks evolve around corotation. The distance between ultraharmonics increases with increasing radius of corotation and decreasing number of spiral arms (Equation (5)). This study is limited to spirals with m = 4 symmetry and so explores the consequence of the spacing between the ultraharmonics as a function of radius of corotation only.

Radial migration from a single episode of cold torquing depends on the radial range of the corotation region (Daniel & Wyse 2015), which scales with the strength of the spiral perturbation (Sellwood & Binney 2002). The corotation region is radially more broad toward the galactic center for spirals of the same fractional amplitude of exponential surface density profile. In models that use an inverse radial surface density profile, the radial range for the corotation region closely follows the trend in the radial distance between ultraharmonics. The following discussion does not argue a preference for an underlying model. Rather, it explores the role resonant overlap has on limiting the radial extent of cold radial migration from cold torquing.

Figures 3 and 4 illustrate how the difference between these two models for spiral amplitude affect the fraction of the corotation region that overlaps with the ultraharmonics. The horizontal lines in Figures 7 and 8 show the distance between the corotation radius and an ultraharmonic (solid, gray). The curves for ${\left\langle {\left({\rm{\Delta }}{L}_{z}\right)}^{2}\right\rangle }^{1/2}$ for Type 1 orbits never exceed ΔR⁽²⁾_LR/2 (Equation (5)) even when the radial range for the corotation region is greater than the annular range between the ultraharmonics. There is therefore a limit to the rms change in orbital angular momentum, and thus radial changes, for a population of stars migrating from cold torquing that is set by the spacing between ultraharmonics. The spacing between ultraharmonics is more restrictive toward the galactic center. The linear approximation that maximum changes in orbital angular momentum from cold torquing are set by spiral strength must be further constrained by the nonlinear response at resonant overlap.

Trapped orbits that experience resonant overlap but remain trapped (Type $3\,\to$ 1) are also limited by the constraint that ${\left\langle {\left({\rm{\Delta }}{L}_{z}\right)}^{2}\right\rangle }^{1/2}$ $\leqslant \,{\rm{\Delta }}{R}_{\mathrm{LR}}^{(2)}/2$. However, orbits that begin trapped at corotation but end in non-trapped orbits after an episode of resonant overlap (Type $3\to 2$) can have changes in orbital angular momentum that exceed this limit. The implication is that strong transient perturbations can induce large changes in ${\left\langle {\left({\rm{\Delta }}{L}_{z}\right)}^{2}\right\rangle }^{1/2}$, but these changes could be dominated by a large fraction of stars that are kinematically heated due to resonant overlap, even when their changes in orbital angular momentum are around corotation.

It is worth considering that the trends shown in Figure 6 include regions at lower galactocentric radii where, in the Milky Way, the kinematics would be dominated by the bar. Assuming a circular velocity of 220 km s⁻¹ and bar pattern speed 41 ± 3 km s⁻¹ kpc⁻¹ (Bovy et al. 2019), the Milky Way's bar has a corotation radius equal to ∼5.4 kpc (see also Wegg et al. 2015; Portail et al. 2017). The shape of the underlying potential, and resulting shape of the corotation region, for a bar is rather different from the case of a spiral pattern. Additionally, the vertical component cannot be taken to be separable as can be done with a spiral in a thin disk and as is assumed in the formulation of the capture criterion (Daniel & Wyse 2015). Nonetheless, there is a significant degree of resonant overlap between the ultraharmonics and the corotation region for a bar. The regions around the L4 and L5 points, which are at the corotation radius and coincident with the azimuthal line through the bar's minor axis, are typically considered stable. However, in a given barred disk there may be additional resonances between the galactic center and the corotation radius of the bar. In this scenario orbits that would be considered stable in the approximation that corotation is at a particular radius would be subject to resonant overlap when recognizing that the corotation resonance fills some volume in phase space. Such resonant overlap would presumably induce a chaotic response, where chaotic regions in phase space are expected to be depopulated (Pfenniger 1990). Indeed, Buta (2017) uses the evacuation of orbits in these regions in the so-called "gap method" to infer the location of corotation and thus other resonances of the bar. These dynamics are currently under further investigation. This discussion is relevant in the current study only insofar as to recognize that the kinematic trend does not reflect the kinematics from resonant overlap from a bar.