GLOBAL SPIRAL ARMS FORMATION BY NON-LINEAR INTERACTION OF WAKELETS

The present work focuses on the formation of spiral arms in disk galaxies. Many researchers have investigated the formation and evolution of spiral arms using analytic methods and numerical simulations. Despite these efforts, our understanding of the formation mechanism of spiral arms is not yet complete.

The first major approach for this topic is the density waves theory advocated by Lindblad (1960) and Lin & Shu (1964). This theory suggests that spiral structure is not a material arm but is a quasi-stationary density wave that propagates through the galactic disk with a constant pattern speed that does not depend on galactocentric radius. It is easy for this theory to explain the long-lived spiral structure, although its origin cannot be answered.

Another possibility is the swing amplification advocated by Toomre (1981). Swing amplification theory suggests that spiral arms evolve when a density-enhanced structure is wound up by differential rotation of the galactic disk (Goldreich & Lynden-Bell 1965; Julian & Toomre 1966; Toomre 1981). The key element of this mechanism is that the rotational direction of epicycle motion coincides with that of the winding spiral. Through this coincidence, stars in spiral arms remain in high-density regions for a long time, and are influenced strongly by gravitational force from spiral arms. This effect causes the rapid growth of spiral arms.

Maximum amplification of density enhancement by this mechanism is examined in previous works. Toomre (1981) analyzed the maximum amplification factor as a function of X and Q for three disk models. X is the ratio of the azimuthal wavelength of the spiral pattern and the critical wavelength for local instability,

$$\begin{eqnarray}&&X=\displaystyle \frac{{\lambda }_{\phi }}{{\lambda }_{c}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\lambda }_{c}=\displaystyle \frac{4{\pi }^{2}G{\rm{\Sigma }}}{{\kappa }^{2}},\end{eqnarray} \tag{ 2 }$$

where Σ and κ denote surface density and epicycle frequency. Q is the local disk instability criterion (Toomre 1964). Toomre (1981) showed that spiral arms are most developed when X ∼ 2.

Dobbs & Baba (2014) performed a similar calculation that builds on Toomre (1981) and Athanassoula (1984). They added the shear rate of the galactic disk as a new parameter, and showed that spiral arms develop at smaller X when the shear rate is smaller (see Figure 5 of Dobbs & Baba 2014).

Both simulated and observed galaxies show the tendency for the pattern of spiral arms to depend on galaxy disk properties, disk mass, and shear rate. Carlberg & Freedman (1985) found that the wavenumber of the most developed spiral arms is inversely proportional to the disk-to-total-mass ratio. They discussed that this relation results from the existence of a characteristic wavelength with X ∼ 2. D'Onghia (2015) showed that there is an agreement between the analytic prediction and the simulations in terms of the number of spiral arms according to swing amplification. Seigar et al. (2005, 2006) derived the relation of pitch angle and shear rate from observed galaxies, while Grand et al. (2013) and Michikoshi & Kokubo (2014) found a similar relation from simulated galaxies, with stronger shear corresponding to more tightly wound arms.

Density wave theory tries to explain spiral arms as longstanding quasi-static structures, with their amplitude and pattern speed nearly constant with time. However, recently performed N-body simulations suggest a fundamentally different view. These works showed that spiral arms are transient structures and alternate between formation and decay (Carlberg & Freedman 1985; Bottema 2003; Fujii et al. 2011; Sellwood 2011; Grand et al. 2012a, 2012b; Baba et al. 2013; D'Onghia et al. 2013; Roca-Fàbrega et al. 2013). This picture is partly similar to swing amplification theory, but some features of the recently simulated arms cannot be explained in terms of swing amplification. For example, some work proposed that spiral arms co-rotate with stars at each radius (Wada et al. 2011; Grand et al. 2012b; Baba et al. 2013; D'Onghia et al. 2013; Roca-Fàbrega et al. 2013). So swing amplification cannot grasp all the aspects of spiral arm formation in simulated galaxies.

Another approach was taken by D'Onghia et al. (2013), which adds further complexity to the formation mechanism of spiral arms. Their simulation contains many perturbers in low-mass N-body stellar disks. Perturbers are modeled as point masses, each of which is heavier than individual N-body particles. Their gravity induces local density enhancement around each perturber, which they call "wakelets." Their proposition is that global spiral arms are formed by connections of wakelets. Their simulation, however, hampers a deeper understanding of the spiral arm formation. Each simulation contains 1000 perturbers distributed with the same profile as the disk and assumed to be co-rotating on circular orbits. It is difficult to isolate interaction of a certain pair of perturbers and analyze the connection process of wakelets because of this complexity.

The purpose of our study is to overcome this difficulty and clarify the fundamental mechanism of the connection of wakelets. We performed N-body simulations with perturbers used in a similar way to D'Onghia et al. (2013). However our simulation is better controlled. We introduce a smaller number of perturbers, and arrange them regularly in a pair of concentric rings in the disk plane. This setup enables isolating each connection process of wakelets, leading to a more straightforward interpretation of numerical results.

We found that two wakelets, which are orbiting at different galactocentric radii, interact non-linearly with each other when the inner wakelet overtakes the outer one rotating more slowly around the galactic center. Density enhancement caused by this non-linear interaction connects the two wakelets, thereby forming a longer density enhancement. Successive operation of such an interaction is considered to create a global spiral arm extending over the entire disk. The wavenumber of spiral arms developed by this mechanism is found to be consistent with the prediction of the swing amplification. We also investigated how the strength of wakelet connection depends on disk mass and shear rate. The results can be naturally understood by invoking swing amplification.

The rest of this paper is organized as follows. In Section 2, we describe simulation models and the numerical methods. We show our results and analyze the effects of non-linear interaction in Section 3. In Section 4, we show the dependence on disk mass and shear rate. Finally, we discuss the role of non-linear interaction in the formation and evolution of spiral arms in Section 5.

2.1. Galaxy Models

In our simulation, each model galaxy consists of a static dark matter halo and a three-dimensional N-body exponential disk. Several works show that the characteristic properties of spiral arms depend on disk mass and shear rate, as mentioned in the introduction. Therefore, we introduce these two parameters in our disk models.

Shear rate is defined as

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{1}{2}\left(1-\displaystyle \frac{R}{{V}_{{\rm{c}}}}\displaystyle \frac{{{dV}}_{{\rm{c}}}}{{dR}}\right),\end{eqnarray} \tag{ 3 }$$

where V_c is the rotation velocity at radius R. In our model, a rotation curve is prescribed so that the shear rate is constant for all radii. In this case, the rotation velocity becomes ${V}_{{\rm{c}}}\propto {R}^{-{C}_{\mathrm{shear}}}$, where

$$\begin{eqnarray}&&{C}_{\mathrm{shear}}\equiv 2{\rm{\Gamma }}-1.\end{eqnarray} \tag{ 4 }$$

According to Seigar et al. (2005), shear rate of real galaxies is 0.2 < Γ < 0.8, so −0.6 < C_shear < 0.6. For the fiducial case that V_c = 200 km s⁻¹ at R = 8 kpc,

$$\begin{eqnarray}&&{V}_{{\rm{c}}}=200{\left[\displaystyle \frac{R}{8{\rm{kpc}}}\right]}^{-{C}_{\mathrm{shear}}}\;\mathrm{km}\;{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 5 }$$

We modified this rotation curve to

$$\begin{eqnarray}&&{V}_{{\rm{c}}}=\displaystyle \frac{200}{1+{({R}_{{\rm{d}}}/R)}^{2}}{\left[\displaystyle \frac{R}{8{\rm{kpc}}}\right]}^{-{C}_{\mathrm{shear}}}\;\mathrm{km}\;{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 6 }$$

where R_d is the factor which prevents rotation velocity from diverging at the center when C_shear > 0.

We take an exponential profile as the density profile of the stellar disk component. The rotation velocity arising from disk gravity is then

$$\begin{eqnarray}&&{V}_{{\rm{c}},\mathrm{disk}}^{2}=\displaystyle \frac{2{{GM}}_{\mathrm{disk}}{y}^{2}}{{R}_{{\rm{s}}}}[{I}_{0}(y){K}_{0}(y)-{I}_{1}(y){K}_{1}(y)],\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&y\equiv \displaystyle \frac{R}{2{R}_{{\rm{s}}}},\end{eqnarray} \tag{ 8 }$$

where R_s is the scale radius and I and K are modified Bessel functions (Binney & Tremaine 2008). Here we assumed that R_s = 3 kpc. Strictly speaking, this velocity profile is for an infinitesimally thin 2D exponential disk. We use this profile for our 3D simulation, because it is not necessary to setup a exact equilibrium initially. Actually we evolve an isolated model for some time before we introduce perturbers, and make the disk relax into equilibrium in a practical sense.

We use five models that have different disk masses or shear rates. The parameters of each model are listed in Table 1. First, we show results for the standard model in Section 3. In Section 4, we describe other models and discuss parameter dependence.

Table 1.  Galaxy Models

Model	Cshear	Mdisk
Standard model	0.0	1.5 × 1010M⊙
Low mass model	0.0	0.8 × 1010M⊙
High mass model	0.0	3.0 × 1010M⊙
Weak shear model	−0.3	1.5 × 1010M⊙
Strong shear model	0.3	1.5 × 1010M⊙

Download table as:  ASCIITypeset image

V_c and V_c,disk are calculated once we specify the C_shear and M_disk. Then the rotation velocity caused by the dark matter halo is given by

$$\begin{eqnarray}&&{V}_{{\rm{c}},\mathrm{halo}}^{2}={V}_{{\rm{c}}}^{2}-{V}_{{\rm{c}},\mathrm{disk}}^{2}.\end{eqnarray} \tag{ 9 }$$

These velocities for the standard model are shown in Figure 1.

Zoom In Zoom Out Reset image size

The disk-to-total-mass ratio f_d is the important parameter because this parameter is strongly correlated with the number of spiral arms (Carlberg & Freedman 1985). m(R), which is the number of spiral arms as a function of radius, is expected to be given by

$$\begin{eqnarray}&&m(R)\sim \displaystyle \frac{1}{{f}_{{\rm{d}}}(R)}\sim {\left(\displaystyle \frac{{V}_{{\rm{c}}}(R)}{{V}_{{\rm{c}},\mathrm{disk}}(R)}\right)}^{2},\end{eqnarray} \tag{ 10 }$$

and shown in Figure 1.

We model only the disk component as an N-body system (N = 3 × 10⁵), and treat the halo as a static gravitational field. Initial velocities of disk stellar particles are determined by solving Jeans' equation following Hernquist (1993).

2.2. Numerical Method

We use the GRAPE system of the National Astronomical Observatory of Japan for numerical computation.

First, we make up a stable disk by using a high Q value, Q_min = 1.7. We evolved the system for about 8 Gyr from this state to remove initial fluctuations that arise because the disk is not initially in a rigorous dynamical equilibrium. Figure 2 shows the development of  the Q value for this phase. Eventually, the minimum value of Q increases to about 2. The final disk does not have any significant spiral structures, and we take it as a stable disk.

Zoom In Zoom Out Reset image size

We added perturbers to this stable disk in the next step. Perturbers are expected to play a role as the seed to form the wakelets in our numerical experiments. The mass of each perturber is 5.0 × 10⁷M_⊙ (about 0.3% of the disk mass), and this choice does NOT have any astrophysical grounding (e.g., giant molecular clouds). We discuss  possible astronomical origins of wakelets in Section 5.2. In order to isolate the dynamical behavior of perturbers, the arrangement of perturbers is especially important.

Perturbers are placed equally spaced in a circle or two circles around the disk center. We calculated four models as follows. Models A and B are single ring models. The former has perturbers located at 6 kpc, while the latter has perturbers at 9 kpc. Model C is a double-ring model in which perturbers are located at 6 kpc and 9 kpc (see Figure 3). The number of perturbers at each circle was varied from 3 to 10. In these models, each perturber is made to move on a circular orbit with the same velocity as the initial stellar rotational velocity at the same radius. Namely, perturbers co-rotate with the equilibrium stellar disk. Additionally, we calculated Model D, which has no perturber. First, we show the case in which six perturbers are placed in each ring.

Zoom In Zoom Out Reset image size

Figure 3 shows snapshots of the surface density distribution for the standard model. Models A, B, C, and D are displayed from left to right. In models with perturbers, six perturbers are placed in each ring. For each model, density distribution is shown at three different times from top to bottom. We added perturbers at t = 0.

Large dots indicate the location of the perturbers at each time, with green points showing the perturbers that were located at x > 0 and y = 0 at T = 0 Myr.

Note that a density enhancement appears around each perturber. These structures can be considered wakelets. As the disk evolves, further density enhancements are created by the interaction of two wakelets. For example, we can see that global spiral arms are formed by the connection of wakelets in model C.

We used a Fourier transform,

$$\begin{eqnarray}&&A=\displaystyle \frac{1}{{N}_{f}}\displaystyle \sum _{k}\displaystyle \sum _{j}\mathrm{exp}[{im}({\phi }_{j,{t}_{k}}-{{\rm{\Omega }}}_{p}{t}_{k})],\end{eqnarray} \tag{ 11 }$$

to analyze the pattern speed of spiral arms. A, N_f, ${\phi }_{j,{t}_{k}}$, and Ω_p are Fourier amplitude, the particle number used for Fourier transform, the azimuth angle of j-th particle at t = t_k and the pattern speed, respectively. We divided the disk into a series of concentric annuli, and the particles located in each annulus were used for Fourier transformation. We thus find the Fourier amplitude and the pattern speed to be functions of the radius. The width of each annulus is 0.25 kpc, and the time span of  the window function is 0 Myr < t < 300 Myr.

Figure 4 shows the results of this Fourier transform for Models A, B, and C. It is shown that density enhancement is made around each perturber indicated by green dots. Namely, pertubers also form wakelets around themselves in Fourier space. It is also clear that each wakelet has a radial extension and has a roughly constant pattern speed. Figure 4 shows that the radial extent of each wakelet is limited by inner and outer Lindbrad resonances. This feature is similar to that of local modes described by Sellwood & Carlberg (2014). The relation of our wakelets and those local modes is discussed in Section 5.

Zoom In Zoom Out Reset image size

We are interested in a possible "non-linear interaction of wakelets." Here, what "non-linear" means is the effect other than linear superposition of effects of inner and outer perturbers. We use the phrase "non-linear interaction of wakelets" in order to emphasize that interaction of wakelets causes the density development more than linear superposition of wakelets. We devised the following method for picking up the effects of non-linear interaction.

We combine the distribution of particles for the four models (Models A, B, C, or D) as follows. Superimposition of Models A and B would have simply caused the effects of inner and outer perturbers to be linearly combined. On the other hand, superimposition of Models C and D would have an additional effect on the non-linear interaction between the inner and outer wakelets. Therefore, the "enhanced surface density" profile,

$$\begin{eqnarray}&&\delta {\rm{\Sigma }}=({{\rm{\Sigma }}}_{{\rm{C}}}+{{\rm{\Sigma }}}_{{\rm{D}}})-({{\rm{\Sigma }}}_{{\rm{A}}}+{{\rm{\Sigma }}}_{{\rm{B}}}),\end{eqnarray} \tag{ 12 }$$

gives the enhancement by non-linear interaction. In other words, if wakelets do not interact with each other, Equation (12) shows that δΣ ∼ 0 because Model C shows the linear superimposition of the effects of inner and outer perturbers in that case.

We calculate the enhanced surface density at each time step. Figure 5 show the results for the case when the number of perturbers at each circle is six. Note that non-linear interaction makes a density enhancement between the locations of inner and outer perturbers, in particular as shown at t = 300 Myr. This enhancement connects the inner and outer wakelets temporarily.

Zoom In Zoom Out Reset image size

Fourier transform of this "enhanced surface density" gives the pattern speed of non-linear enhanced structure. In order to see the effect of varying the number of perturbers, we carried out the simulations by placing 3–10 perturbers at each ring. The aim of this numerical experiment is to investigate the most developed wavenumber by non-linear interaction.

The number of spiral arms formed by non-linear interaction is equal to the number of perturbers at each circle. Figure 6 shows the Fourier amplitude for non-linear enhancement as seen in Figure 5 when the number of perturbers at each circle is varied from 3 to 10. The width of the time window function for Fourier transform is 100 Myr centered around t = 300 Myr. These results are obtained by calculating Equation (11) for δ Σ. We use the number of perturbers at each circle as the azimuthal wavenumber m for Fourier transform; Equation (11).

Zoom In Zoom Out Reset image size

When the number of perturbers at each circle is 6 (m = 6), the Figure 6 shows a strong amplification between 6 and 9 kpc. This amplification suggests the presence of a non-linearly developed structure. More specifically, global spiral arms of Model C as shown in Figure 3 include not only simple linear superposition of inner and outer wakelets but also non-linear interaction of those wakelets. Note that this structure has a pattern speed that is between those of the two perturbers. Comparing the results for different m suggests that spiral arms are most developed by non-linear interaction when there are 5 or 6 perturbers at each radius for our standard model. Note that these wave numbers (m = 5, 6) are consistent with m(R) between 6 and 9 kpc in Figure 1. It is very interesting that non-linear interaction does not always develop when two wakelets encounter each other. But this result does NOT mean that the number of spiral arms depends on the number of perturbers. Figure 6 only shows that global spiral arms whose wavenumbers are 5 or 6 are selectively developed by non-linear interaction, while spiral arms with other wavenumbers grow inefficiently. Thus the number of global spiral arms developed by non-linear interaction is 5 or 6, independently of the number of perturbers or wakelets. This prediction is consistent with the results of D'Onghia et al. (2013) in which the number of spiral arms does not depend on the number of perturbers.

These non-linear enhancements via the interaction of wakelets are suggestive of the non-linear mode coupling theory suggested by Tagger et al. (1987) and the presence of mode coupling between bar and spiral arms demonstrated by some simulations (e.g., Masset & Tagger 1997; Quillen et al. 2011). Spiral arms may be the product of mode coupling of wakelets.

We also simulated other disk models with different mass or shear rate. We investigated the relationship of the wavenumber and amplitude of non-linear enhancements for these disks in the same way as for the standard model and found the wavenumber of the most developed structure for each disk. Hereafter we refer to this wavenumber as the characteristic wavenumber.

First, we simulated higher-mass and lower-mass disk models than the standard model. Figures 7 and 8 are the results of Fourier transform for "enhanced surface density" of low-mass and high-mass models. Comparing three models, the standard, the high-mass, and the low-mass models, shows that the characteristic wavenumber depends on the disk mass. It is clear that the lower-mass disk develops spiral arms with larger wavenumbers as a result of non-linear interaction of wakelets. For the low-mass model, non-linear enhancement is the strongest at m = 6–8, whereas for the high-mass model non-linear enhancement is the strongest at m = 4. This dependence is consistent with the tendency observed in N-body simulations of disk galaxies (Carlberg & Freedman 1985; Bottema 2003).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, we carried out simulations introducing stronger and weaker shear in the disk. Figures 9 and 10 give the results of Fourier transform for "enhanced surface density" in the weak-shear and strong-shear models. These results show that weaker shear leads to a smaller number of spiral arms. We discuss this dependence on shear rate in Section 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To summarize, we found that the characteristic wavenumber of spiral arms depends on the disk properties, disk mass, or shear rate. A disk with lower mass or weaker shear develops a smaller number of spiral arms through non-linear interaction of wakelets.

Our simulation clearly detected non-linear interaction of wakelets originally caused by perturbers rotating in the disk. It was also demonstrated that this non-linear interaction plays a fundamental role in the formation of global spiral arms. Here we discuss the relevance of these findings to the previous results and try to interpret numerical results reported by other authors that seemingly contradict each other.

5.1. The Dependence on Mass and Shear

5.2. Ubiquity of Non-linear Interaction

5.3. Co-rotation Wave and Local Mode

5.4. Other Spiral Structure

We analyzed the interaction of wakelets by using N-body simulations including perturbing point masses, which are heavier than individual N-body particles and act as the seeds for wakelets. Consequently, we obtained the following results.

• 1.  
  Two adjacent wakelets, which are orbiting at different galactocentric radii, interact non-linearly with each other when the inner wakelet overtakes the outer one rotating more slowly around the galactic center. This non-linear interaction makes density enhancement and connects the two wakelets, thereby creating a global spiral arm extending over the entire disk. (See Sections 3 and 5.2.)
• 2.  
  The wavenumbers of spiral arms developed by this mechanism depend on disk mass and shear rate. This dependence is consistent with the prediction of the swing amplification and suggests that the structure formed by the connection of wakelets experiences swing amplification and develops into global spiral arms. (See Sections 4 and 5.1.)
• 3.  
  In our simulation, non-linearly enhanced structures have pattern speeds between those of the inner and outer perturbers. It is also located spatially between the two perturbers that co-rotate with disk stars. This result provides unification of previous results, namely local wave modes (Mata-Chávez et al. 2014; Sellwood & Carlberg 2014) and spiral arms co-rotating with stellar particles at each radius (Wada et al. 2011; Grand et al. 2012b; Baba et al. 2013; D'Onghia et al. 2013; Roca-Fàbrega et al. 2013). (See Section 5.3.)

Numerical computations in this paper were performed on the GRAPE system at the Center for Computational Astrophysics, National Astronomical Observatory of Japan.