7 Phase space crowding

Given the orbital structure of phase space, we now want to know how nature populates the available orbits. This is done resorting to test particle simulations with the following two integral initial distribution function (Dehnen 1999a):

$$f_{\circ}(E,L_z)\!=\!\frac{2\Omega(R_{\rm c})}{\kappa(R_{\rm ... ...c{L_z\!-\!L_{\rm c}(E)} {\tilde{\sigma}^2(R_{\rm c})}\right]},$$ (12)

where E and L_z are respectively the energy and the z-component of the angular momentum, $R_{\rm c}(E)$ and $L_{\rm c}(E)$ the radius and angular momentum of the circular orbit with energy E, $\Omega$ and $\kappa$ the circular and epicycle frequencies, and $\tilde{\Sigma}(R)$ and $\tilde{\sigma}(R)$ the approximate surface density and radial velocity dispersion profiles. This is a modified Shu (1969) distribution where the radius of the guiding centre $R_{\rm c}$ is deduced from the energy instead of the angular momentum, with the main advantage that the density function extends smoothly towards negative L_z. Adopting $\tilde{\Sigma}(R)\sim \exp{(-R/h_{\rm R})}$ and $\tilde{\sigma}^2(R)= \tilde{\sigma}_{\circ}^2\exp{(-(R-R_{\circ})/h_{\sigma})}$, where $R_{\circ}$is the galactocentric distance of the Sun and $\tilde{\sigma}_{\circ}$ the approximate velocity dispersion at that distance, the initial distribution function has three free parameters, which unless otherwise specified are set to $h_{\rm R}/R_{\circ }=0.33$, $\tilde{\sigma}_{\circ}/v_{\circ}=0.2$ and $h_{\sigma}/R_{\circ}=1.0$. This results in exactly the same initial conditions as for the default flat rotation curve case in D2000. In the beginning of the simulations, the non-axisymmetric part of the potential is gradually switched on from no contribution at t=0 to its full value at $t=2t_{\rm b}$, where $t_{\rm b}=2\pi/\Omega_{\rm P}$ is the rotation period of the bar, exactly the same way as in D2000 for the simulations with default integration time.

For the time integration, D2000 uses a subtle backward integration technique based on the conservation of the phase space density along the orbits in collisionless systems. The idea is to integrate back in time until t=0 the phase space points on a Cartesian grid of u-v velocities at a given space position $(R_{\circ},\varphi)$ and time $t_{\rm end}$, to determine the energy E and angular momentum L_z of the originating orbits in the initial axisymmetric potential, and from these infer $f(t_{\rm end},R_{\circ},\varphi,v,u)=f_{\circ}(E,L_z)$. The advantages of this technique is that only the orbits strictly necessary to derive the evolved local velocity distributions need to be computed, and that these velocity distributions are not affected by Poisson noise. Moreover, a unique simulation suffices to get u-v distributions for different initial conditions, because $f_{\circ}$ comes in only after the orbit integration.

Unfortunately, the backward integration technique faces two major problems illustrated in Fig. 11, which shows the long term evolution of the planar velocity distribution at $R_{\circ}/R_{\hbox{\tiny OLR}}=1.1$ and $\varphi =25^{\circ }$ using this technique. The integration time in D2000 ranges from 4 bar rotation for most simulations, corresponding to only $\sim$2 orbital periods at $R=R_{\circ}$in the inertial frame, up to 8 bar rotation for some cases, but always matching the growth time of the bar to half the total integration time. The frame at $t=4t_{\rm b}$ is similar to his results, revealing a clear bimodal distribution. However, at $t=10t_{\rm b}$, the valley between the two modes becomes heavily populated, destroying the bimodality, and at later times, incurved waves appear in this valley with a spacing between the maxima decreasing with time. This is a typical signature of ongoing phase mixing in a regular region of phase space, which here corresponds to the eccentric orbit part of the quasi-x₁(1) region according to the previous section. A very similar phenomenon, with similarly incurved waves, can also be barely recognised near the 1/1 resonance. Between the 2/1 and 1/1 resonances, phase mixing operates on a shorter timescale and the orientation of the wave fronts seems to change from nearly-vertical at $t=30t_{\rm b}$ to nearly-horizontal at $t=120t_{\rm b}$. The backward integration technique in fact yields the fine-grained distribution function, which never smoothes out on sufficiently small scales, whereas the physical one to compare with the observations is the coarse-grained distribution. The second problem is related to chaos: at , the u-v distribution becomes noisy in the chaotic regions because the phase space points integrated backwards from these regions sample the initial distribution function more randomly. Hence much longer integration times than in D2000 are required to obtain quasi-equilibrium (coarse-grained) distribution functions and to properly take into account the effect of chaos (remember that the divergence timescales for chaotic orbits is of the order of several orbital periods), and one cannot escape the fate of smoothing the fine-grained distribution.

Figure 11: Time evolution of the velocity distribution in the u-v plane at $R_{\circ}/R_{\hbox{\tiny OLR}}=1.1$ and $\varphi =25^{\circ }$ using the backward integration technique. The bar strength is F=0.10 and the time is given in units of the bar rotation period $t_{\rm b}$.

Therefore, most of the test particle simulations in this paper were done by simple forward integration. The initial phase space density is sampled by N particles which are then integrated forward in time, and the u-vdiagrams at space position $(R_{\circ},\varphi)$ are constructed from all the particles within a distance $d_{\rm max}=R_{\circ}/80$ from that position (corresponding to $d_{\rm max}=100$ pc for $R_{\circ}=8$ kpc), using a Cartesian velocity binning with a bin size of $0.005v_{\circ}$. To increase the particle statistics, the u-v distributions are averaged over 10 bar rotations and then smoothed within a square of $11\times 11$ bins. The time average, which is hardly possible in the backward integration technique, is also very convenient to reduce the phase mixing problem, as the contribution of each particle to a u-v distribution becomes proportional to the time the particle spends within the volume where the distribution is computed. With the forward integration technique, the time evolution of the velocity distribution can be followed within a unique simulation. The test particle simulations of this paper all have N=10⁶ and the velocity distributions have been derived within two time intervals, $25\leq t/t_{\rm b}\leq 35$ and $55\leq t/t_{\rm b}\leq 65$. Since the distance parameters in $f_{\circ}$scale as $R_{\circ}$ and not $R_{\hbox{\tiny OLR}}$, the results at different $R_{\circ}/R_{\hbox{\tiny OLR}}$ require distinct simulations.

Figure 12 shows the u-v distributions at various space positions averaged over the time interval $25\leq t/t_{\rm b}\leq 35$ for a bar strength F=0.10. As for the Liapunov diagrams, the distributions are obviously symmetric with respect to u=0 for $\varphi=0$ and $\varphi=90^{\circ}$ and the distributions at same radius but supplementary angles are anti-symmetric to each other in u. This is clearly not the case in the simulations of D2000 (see his Fig. 2, where F=0.089), providing a further argument that these have not achieved a quasi-stationary regime. Moreover, the traces of the stable x₁(1) periodic orbits away from the 2/1 resonance curve lie closer to the high angular momentum peak of the velocity distributions. For this bar strength and the adopted values of the parameters in the initial distribution function, there is also no clear bimodality with a deep separation valley at and near $\varphi =30^{\circ}$, although a clear density excess remains at low v and positive u. However, all space positions where a low-eccentricity regular x₁(2) region exists (see former section and Fig. 9) present a nice bimodality, with the low angular momentum mode coinciding very well with that region and always peaking inside the H₁₂contour.

While the traces of the non-resonant x₁(1) and x₁(2) orbits are generally embedded within their associated quasi-periodic orbit modes, they do not necessarily exactly coincide with the peak of these modes, especially in the x₁(2) case (e.g. $R_{\circ}/R_{\hbox{\tiny OLR}}=1.0$ and $\varphi=60^{\circ}$- $120^{\circ}$ in Fig. 12). Furthermore, since the quasi-periodic orbits cover a larger space area than the periodic orbits themselves, quasi-x₁(1) and quasi-x₁(2) modes may occur even at positions where no x₁(1) or x₁(2) orbit are passing through.

Increasing the bar strength (Fig. 13) provides a better understanding of how the velocity distributions are affected by chaos. Now, there appears to be an obvious second source producing a low angular momentum mode, which adds to the quasi-x₁(2) orbit flow at the space regions reached by these orbits, and acts alone elsewhere, as for instance at $R_{\circ}/R_{\hbox{\tiny OLR}}=1.1$ and $\varphi =30^{\circ}$. The overdensity in velocity space generated by this second source correlates very well with highly stochastic regions in the Liapunov diagrams (see Fig. 8) and seems to always peak outside the H₁₂ contour. At $\varphi=90^{\circ}$, the overdensity culminates at $u\approx 0$ and is enclosed by the regular arc of eccentric quasi-x₁(1) orbits. At $\varphi=0$, these quasi-x₁(1) orbits occupy the u=0 region and chaotic overdensities happen symmetrically on both positive and negative u sides of this region. At $\varphi \sim 30^{\circ}$, the quasi-x₁(1) region is located at negative u and there is a large chaotic overdensity at positive u.

These properties result from the decoupling between the regular and the chaotic regions of phase space. Since chaotic orbits cannot visit the regions of regular motion and, vice versa, regular orbits avoid the chaotic regions, the distribution function in each of these regions evolves in a completely independent way. In the regular regions, it recovers roughly the initial distribution after phase mixing, whereas in the chaotic regions, it is substantially modified through a process known as chaotic mixing and which operates on the Liapunov divergence timescale (e.g. Kandrup 2001): the particles on chaotic orbits quickly disperse within the easily accessible phase space regions, i.e. not impeded by cantori or an Arnold web, and converge towards a uniform population of these regions. The dominant manifestation of chaotic mixing is a net migration of particles from the inner to the outer space regions. For instance, in the simulations with $R_{\circ}/R_{\hbox{\tiny OLR}}=1.1$, the scale length of the radial particle distribution, which remains very close to an exponential in the range , increases by $\sim$30% for F=0.10 and by $\sim$90% for F=0.20 within this radial range. This migration is particularly marked for particles on hot chaotic orbits because the region inside corotation initially represents a large reservoir of such particles and because these particles can freely pass over corotation. As a consequence, in the explored $R_{\circ}/R_{\hbox{\tiny OLR}}$ range, the chaotic regions in the u-v diagrams are more heavily crowded than the regular regions at , therefore corresponding to local overdensities.

At first glance, there seems to be a rather continuous transition from the quasi-x₁(2) orbit mode to the main chaotic orbit mode when moving across the OLR radius towards increasing $R_{\circ}$, with always a single effective peak showing up and with the involved quasi-x₁(2) region progressively dissolving in the chaotic one. But in some cases, the two mode-generating sources really contribute to distinct peaks in the velocity distribution (see Fig. 20d for an example).

The process of chaotic mixing leads to velocity distribution contours which are parallel to the contours of constant H in the chaotic regions (e.g. Fig. 13). This property is also consistent with Jean's theorem stating that the distribution function in a steady-state system depends only on the integral of motions. The only integral for the chaotic orbits in the present 2D barred models is the value of the Hamiltonian, hence the distribution function and therefore the corresponding velocity distributions at fixed space position should be a function of only H in the chaotic regions. It should be noted that the Jeans theorem does not strictly apply to the hot and disc chaotic orbits. Indeed, these orbits are not energetically bound (in terms of H) and thus the phase space density around such orbits and within the finite space volume of the galaxy should decrease with time, conflicting with the steady-state assumption of the theorem. However, the escaping timescale of these chaotic orbits, which is essentially controlled by Arnold diffusion across the confining cantori, is much longer than the Liapunov divergence timescale and even the Hubble time, and thus the density in the phase space regions covered by these orbits can be considered as almost constant (see also Kaufmann & Contopoulos 1996 and references therein).

Figure 12: Velocity distribution in the u-v plane as a function of space position in the test particle simulations with a bar strength F=0.10. The distributions are time averages within the interval from 25 to 35 bar rotations after the beginning of the simulations. The velocity contours are as in Fig. 1, whereas the axis labels, the Hamiltonian contours, the resonance curves and the periodic orbits are as in Fig. 7.

Secondary chaotic orbit overdensities also occur between the x₁(1) and the 1/1 regular arcs, especially at $R_{\circ}/R_{\hbox{\tiny OLR}}\geq 1.0$and $30^{\circ}\leq \varphi \leq 150^{\circ}$ (Fig. 13). These secondary overdensities and the above described main overdensities connect to each other in phase space, i.e. are traced by the same orbits. Hence it is also expected that the u-v density at constant H-value is the same for all overdensities. This is only roughly the case in Fig. 13, probably because the smoothing of the diagrams lowers the peaks of the tiny secondary overdensities relative to the broader main overdensities. Small chaotic overdensities may sometimes even be noticed beyond the 1/1 resonance curve. However, at high angular momentum, the hot orbits spend most of their time in the outer galaxy, far away from the influence of the bar, and thus become more regular (the energy and angular momentum are more nearly conserved). The eccentric x₁(2) regular regions can also represent density depressions between chaotic regions in the velocity distributions, as can be marginally inferred for example from the $R_{\circ}/R_{\hbox{\tiny OLR}}=0.9$ and $\varphi=90^{\circ}$ frame in Figs. 9 and 12.

The valley between the main high-L_z velocity mode (or LSR mode after Dehnen) and the main chaotic orbit mode is generally close to the H₁₂contour, reflecting the decline of the hot orbit population as $H\rightarrow H_{12}$. Such a valley should in principle also exist between the LSR mode and the secondary chaotic overdensities (see for instance $R_{\circ}/R_{\hbox{\tiny OLR}}=1.1$ and $\varphi =30^{\circ}$ in Fig. 13). The main chaotic orbit mode also seems to always peak between the H₁₂ and H₄₅ contours in Fig. 13, but this is not true for all our test particle simulations, as demonstrated by the $h_{\rm R}/R_{\circ}=0.2$ frame in Figs. 14 and 20a. However, this property might be more generic for self-consistent models (see Sect. 10). For some not fully understood reasons, the symmetry properties mentioned previously for the case F=0.10 are somewhat less evident for F=0.20, despite the longer integration time.

D2000 attributes the valley between the main LSR mode and the Hercules-like mode to stars scattered off the OLR, in the sense that the resonance generates chaotic orbits. In particular, he claims that the unstable x₁^*(2) orbit falls exactly between the two modes. This is not quite correct for a stochastically induced Hercules-like mode, as is best illustrated by Fig. 13: for and $\varphi \sim 30^{\circ}$, the u>0 part of the 2/1 resonance curve passes through the Hercules-like mode and the x₁^*(2) orbit clearly lies within the mode at $R_{\circ}/R_{\hbox{\tiny OLR}}=1.2$, and the low density region below this mode is due to regular resonant orbits. D2000 also claims that in his simulations the extension of the LSR mode to u<0 at $v\approx -0.1v_{\circ}$ is caused by the elongation of the (presumably quasi-x₁(1)) orbits near the OLR. Our results indicate that at least the final part of this extension, corresponding to the secondary overdensities, is produced by chaotic orbits. Such an extension exits in the observations, but is only significant down to heliocentric $u\approx -60$ kms^-1.

The velocity distributions at the two different mean integration times $<\!t\!>\,=30t_{\rm b}$ and $60t_{\rm b}$ reveal some secular evolution. As $<\!t\!>$ increases, the crowding contrast between the regular and chaotic regions becomes more evident, with denser high-H chaotic regions and deeper regular region valleys. The quasi-x₁(1) mode squeezes towards its high angular momentum side for and, especially in the case F=0.10, the peak of the quasi-x₁(2) mode moves closer to the H₁₂ contour for $R_{\circ}=R_{\hbox{\tiny OLR}}$, betraying a longer phase mixing timescale near this resonance.

In a real galaxy, the presence of mass concentrations like giant molecular clouds and of transient spiral arms will prevent the strict conservation of the Jacobi integral and cause the stars to diffuse from the regular regions to the chaotic regions of phase space and vice versa (see also Sect. 10). The chaotic regions should be very efficient in heating galactic discs and the communication between the two kind of regions may even allow to heat regular regions. The quantification of this phenomenon might be an interesting problem to study, but is beyond the scope of this paper.

Finally, Figs. 12 and 13 also suggest that with increasing bar strength, the velocity dispersion ratio $\sigma_{v}/\sigma_{u}$decreases and, as reported by D2000, the u-range and in particular the mean u-velocity of the main chaotic overdensity become larger.

Figure 13: Same as Fig. 12, but for a bar strength F=0.20 and a time average over the interval $55\leq t/t_{\rm b}\leq 65$.