© ESO, 2013

1. Introduction

The morphology of spiral galaxies, as laid out in the Hubble classification (Hubble 1926), can be broadly characterised by the tightness of spiral arm structure and the size of the central region or bulge. In this classification scheme, more tightly wound spiral arms are associated with large central mass concentrations. The strong correlation between central mass concentration and pitch angle predicted by modal density wave theory (e.g. Lin & Shu 1964; Roberts et al. 1975; Seiden & Gerola 1979; Bertin et al. 1989) is in accordance with this. However, there are complications in the Hubble classification scheme insofar as that this relation was derived from optical information of galaxies only. The correlation is not observed in the near-infrared wavelengths (de Jong 1996; Seigar & James 1998a,b), and some observational studies in the infrared waveband highlight a difference in morphology from that seen in the optical (e.g. Block et al. 1994; Thornley 1996; Grosbol & Patsis 1998). Moreover, the correlation between Hubble type and pitch angle has been shown to be weak (Kennicutt 1981) and the model predictions from density wave theory for spiral arm properties have been shown to have systematic offsets to observations (Kennicutt & Hodge 1982).

Despite these uncertainties in the Hubble type-pitch angle relation, more recent observations have shown convincing evidence for a correlation between spiral arm pitch angle and the shear rate of differentially rotating discs of spiral galaxies. Seigar et al. (2005) derived shear rates from the rotation curves of a sample of several barred galaxies and used Fourier analysis to draw the spiral shape. They found evidence for the shear rate dependency of the spiral arm pitch angle. Because the rotation curve shape is determined by the mass distribution, this is essentially a correlation between the central mass concentration and spiral arm pitch angle. This survey was later extended and the conclusion strengthened by Seigar et al. (2006).

The shear rate-pitch angle correlation is also supported by the analytical work based on swing amplification theory (Goldreich & Lynden-Bell 1965; Toomre 1981) by Julian & Toomre (1966; see also Fuchs 2001), which calculated the spatial distribution of the response of the density of the differentially rotating stellar disc to a large perturbing mass. They showed that the density enhancement in this context is predicted to show smaller pitch angles (hence a more tightly wound structure) with increasing amount of shear present.

While theoretical and observational studies provide evidence for the shear rate-pitch angle relation, it has yet to be explored in N-body simulations. In this paper, we aim to study this relation by running a suite of N-body simulations of varying shear rates. For the first time we investigate the pitch angles of individual spiral wave mode patterns in N-body simulations by isolating the spiral wave mode patterns from the system using the conventional spectrogram analysis (e.g. Quillen et al. 2011; Sellwood 2012; Solway et al. 2012; Minchev et al. 2012; Roškar et al. 2012) and calculating the spiral phase of the mth mode. We find that the discs of higher shear rate exhibit systematically smaller pitch angles than their lower shear rate counterparts, as predicted from the theoretical studies mentioned above. We also trace the overall spiral arm feature and measure its pitch angle as a function of time. The motivation for exploring this pitch angle behaviour is that we and other authors have found that the pattern speed of the spiral arms in N-body simulations and observed galaxies decreases with radius in a similar manner to the angular rotation velocity of the disc particles (Merrifield et al. 2005, 2006; Speights & Westpfahl 2011; Wada et al. 2011; Grand et al. 2012a,b; Nelson et al. 2012; Comparetta & Quillen 2012; Baba et al. 2013). Because the pattern speed decreases in this way, the pitch angle decreases with time and leads to transient and recurrent spiral arm features that are seen in many simulations (e.g. Sellwood 2010, 2011, and references therein). The evolving nature of the pitch angle of winding spiral arm features can be compared to the observational work of Seigar et al. (2006), which measures the pitch angle and shear rate of many spiral galaxies and reveals several different observed pitch angles for a given shear rate.

The paper is organised as follows. The simulations are described in Sect. 2, the analysis techniques laid out in Sect. 3, and the results are described in Sects. 4 and 5 in which we also explore some of the other parameter space apart from shear rate. The discussion is presented in Sect. 6, followed by the conclusions in Sect. 7.

2. Simulations

The simulations in this paper are performed with a hierarchical Tree N-body code GCD+ (Kawata & Gibson 2003; Kawata et al. 2013). We run a suite of simulations, each of which consists of a spherical static dark matter halo (and a spherical static stellar bulge component in some cases) and a live stellar disc. The halo and bulge are static rather than live in order to facilitate greater control of the experimental scenarios. A live halo/bulge component will complicate the evolution of the stellar disc with effects such as scattering and heating, and may even act as large perturbing masses that greatly disturb the disc if the mass resolution for the dark matter is too small (D’Onghia et al. 2013). These are unwanted effects, and because the focus of this study is on the stellar disc component only, we have elected to model the external components with static potentials.

The dark matter halo density profile follows that of Navarro et al. (1997) with the addition of an exponential truncation term (Rodionov & Athanassoula 2011): (1)where ρ_c is the characteristic density described by Navarro et al. (1997), the concentration parameter, c = r₂₀₀/r_s, and x = r/r₂₀₀. The truncation term, exp   (−x²), is introduced in our initial condition generator for a live halo simulation. Although we use a static dark matter halo in this paper, we retain the profile of Eq. (1) because this term does not change the dark matter density profile in the inner region, which is the focus of this paper. The scale length is r_s, and r₂₀₀ is the radius inside which the mean density of the dark matter sphere is equal to 200ρ_crit (where ; the critical density for closure): (2)where M_vir is the virial mass of the galaxy.

We assume Ω₀ = 0.266, Ω_b = 0.0044, and H₀ = 71   km   s^-1   Mpc^-1.

The spherical static stellar bulge component is modelled by the Hernquist profile (Hernquist 1990), which is described by: (3)where M_b is the total bulge mass and a is the scale length. The scale length is set to the effective radius, R_e = 1.8153a. We apply a compacting factor, b, to scale from the empirical relation of the bulge effective radius (Shen et al. 2003): (4)Hence, the resultant scale-length is defined by a = bR_e/1.8153.

The stellar disc is assumed to follow an exponential surface density profile: (5)The fiducial number of disc particles used is N = 1 × 10⁶. Numbers of this order are reported to be sufficient to minimise numerical heating (Fujii et al. 2011). Although larger particle numbers reduce numerical heating further, we note that the effect is always present (i.e. it does not disappear at a particular resolution), and that a compromise between parameter space and resolution must be made for suites of simulations such as the one presented in this study.

We apply a fixed softening length, ϵ, for star particles with the spline softening suggested by Price & Monaghan (2007). The softening length¹ is dependent on the particle mass, therefore the base value of ϵ = 340 pc for the particle mass, m_p = 5 × 10⁴   M_⊙ varies between simulations that have different particle masses. The model parameters for the simulations are summarised in Table 1, and the rotation curves are shown in Fig. 1.

There are three groups of rotation curves. Simulation group R (R, R2, R3, R4) has a rising rotation curve. Simulation R is an extreme case, where we set a large halo mass with a low concentration parameter, c, in order to extend mass to the outer regions of the disc. Because of such a low concentration of dark matter mass in the central region, the disc mass must be lowered in order to prevent a bar from forming (Ostriker & Peebles 1973). In this way, we avoid the added complication of the bar component and restrict the study to spiral galaxies only. Simulations R2, R3, and R4 are less extreme cases, which explore intermediate shear rates and different disc to halo mass ratios. To produce the flat (simulations F, Fa, Fb, Fc, F2, and F3) and Keplerian-like (simulation K) rotation curves, a bulge component is included. For simulation K, this is a very compact and massive bulge. Although this case is unrealistic, we include it in order to emphasise the effect of galactic shear on spiral morphology.

Fig. 1 Circular velocity at t = 0 for simulation R (very thick dashed red), R2 (thick dashed blue), R3 (medium dashed green), R4 (thin dashed cyan), F (thick dot-dashed green), F2 (thin dot-dashed red), F3 (medium dot-dashed blue), and K (solid blue).

Fig. 2 Galactic shear rate, Γ, for all simulations. Colours are the same as Fig. 1. Note the reduced radial range compared to Fig. 1.

The radial profile of the galactic shear rate at t = 0, given by: (6)for each simulation is shown in Fig. 2. This suite of simulations represents a range of shear rates, which is the principal variable we want to investigate. However, there are other parameters that may affect the pitch angle, such as the disc-halo mass ratio, ζ, softening length, ϵ, and resolution. We also explore these parameters, mainly with simulation group F.

We set the initial Toomre stability parameter, Q, for all our simulations to approximately 1 over the radial range 4 < R < 10 kpc, which allows the spiral structure to grow². The radial dependence in shown in Fig. 3.

3. Method of analysis

Here we present the analysis method of our two techniques for measuring pitch angles: mode pattern analysis and direct spiral arm peak trace method. An important difference between these techniques is that the mode pattern analysis assumes that the spiral arms are constructed by one or multiple density waves of mode, m, which describe patterns of m spiral arms with a constant pitch angle. The direct spiral arm peak trace method does not assume any theory, but simply analyses the pitch angle of the overall spiral arm feature. The distinction between these two methods is that while both characterise the spiral arm as a logarithmic spiral of fixed pitch angle at all radii of interest, in the direct method the pitch angle and amplitude of the spiral arm changes with time. However, in the mode analysis, changes in the spiral arm (in particular the winding) may only occur through the changing superposition of the various mode patterns present.

Fig. 3 Toomre stability parameter, Q, at t = 0 for all simulations. Colours are the same as Fig. 1.

Before we describe these two analysis techniques, we define the pitch angle which we will use with both. Given the positional information (R,θ) of a density enhancement, we can fit logarithmic spiral arms, described by: (7)where θ is the azimuth coordinate, R is the radial coordinate, and B and C are constants. Logarithmic spirals have pitch angles, φ, given by (Binney & Tremaine 2008): (8)where the distance, d_θ, is the spatial distance of the density enhancement in the azimuthal direction defined as d_θ = RΔθ. The pitch angle of a logarithmic spiral is constant with radius. The next step is to recover the positional information (R,θ) required to apply the logarithmic chi-squared fitting using Eq. (7) and calculate the pitch angle of the fit using Eq. (8).

3.1. Mode pattern analysis method

By construction, a wave mode pattern has a constant pattern speed, . Therefore the shape of a wave mode pattern is time independent i.e. the pitch angle is constant over time. In this analysis, we focus on strong patterns because their behaviour is most evident. In order to find patterns of significant amplitude, we first search for dominant modes i.e. wave modes of m spiral arms that exhibit large amplitudes. The amplitude of a given wave mode, m, is calculated from the quantities: (9)where θ_i is the azimuthal angle between the radial vector of the particle and a common reference vector. The amplitude is then calculated as: (10)The mean amplitude in a radial range 4−10 kpc is calculated using Eq. (10) for modes m = 1−7 over the entire 2 Gyr of the evolution for each simulation. This is shown in the top row of Fig. 4. In each simulation, prominent modes are identified for analysis. We aim to extract the positional information of the patterns. The adopted procedure is to compute their power spectra by taking the Fourier transform of the time sequence of each component in Eq. (9) (Quillen et al. 2011): (11)where h(t) denotes the Hanning function used to reduce the aliasing. T₁ and T₂ denote the beginning and end of the time window of the Fourier transform. This is chosen to be at around a relatively late epoch of the simulation (when the system is more stable) and is centred around a peak of the most dominant mode present in each case. It spans Δt = 256 Myr, which is a typical life time of a spiral arm as shown in the next section.

The amplitude in each frequency as a function of radius is then calculated via: (12)Because simulations generally possess several patterns for a given mode that can overlap in radius (e.g. see Fig. 4 of Roškar et al. 2012), care must be taken when computing the spiral phase of a pattern. In this technique each wave mode pattern is characterised by a pattern speed given by , which is constant over radius. Individual patterns should be selected by isolating a horizontal ridge (a single pattern speed) over a radial range where the signal significantly stands out from the noise. In each of the galaxies, we focus on the most dominant patterns and look at the three quantities, , , and A^m(R,ω) on the real and imaginary axis for each radial pixel in a ridge. We then calculate the real spiral arm phase position within the domain 0 to 2π as: (13)where is the spiral phase of the pattern at each radial bin, which is retrieved by considering only the Fourier coefficients of a single ω. Because this quantity spans a domain of 2πm, the spiral phase, , is divided by m in order to yield the real phase position of the wave mode pattern as a function of radius. This provides the azimuthal and radial values required for the calculation of the pitch angle using Eqs. (7) and (8).

Fig. 4 Top row: amplitudes calculated from Eq. (10) and averaged over a radial range of 4−10 kpc of spiral modes m = 1 (thick red); 2 (thick green); 3 (thick blue); 4 (thick yellow); 5 (thin red); 6 (thin green); 7 (thin blue), and 8 (thin yellow) normalised to the axisymmetric m = 0 mode, as a function of time for simulations R (left), F (middle), and K (right). Vertical dashed lines represent the time window of a Fourier transform applied in Sect. 3.1. Second row: power spectra calculated from Eq. (12) of simulation R for the m = 8 mode (left), F for the m = 3 mode (middle), and K for the m = 2 mode (right). Prominent ridges (dark pixels) span between 4−10 kpc in most cases. Third row: in polar coordinates, the density map of the dominant density wave mode pattern selected from rows of , 30, and 24 km s-1 kpc-1 for simulations R, F, and K respectively. White regions indicate areas of low density and black regions indicate areas of high density. Contours emphasis the highest density regions. Bottom panels: dominant mode pattern positions (black points) calculated from Eq. (13) in the azimuth-radius plane for the corresponding patterns in the row above. The red lines show the lines of best fit for each pattern. The right side of each panel shows the radial amplitude profile, which is used to weight the fitting.

3.2. Direct spiral arm peak trace method

The method we use to trace the spiral arm peak position directly is a particle density weighting method, in which we select a point near the spiral arm of interest at some start radius (~5 kpc), define an azimuth range that encapsulates the width of the spiral arm and weight by particle density to find the peak position (see Grand et al. 2012a, for more details). This is iterated over a radial range until the spiral arm peak position is drawn out. Several spiral arms are traced over a range of snapshots between 1 and 2 Gyr of the simulation evolution. Spiral arms are only traced when they show a single density peak over azimuth for each radius in the radial range chosen for fitting. The pitch angles are then calculated using Eqs. (7) and (8).

We remind the reader that this pitch angle is derived from the spiral arm line that traces out the overall density enhancement directly, which varies with time. This is different from the time independent pitch angle calculated from the positional information of the wave mode patterns derived from the power spectra (see Sect. 3.1). The latter bears the assumption of a density wave of constant pattern speed and fixed pitch angle, whereas the former bears no assumptions at all.

4. Results of fiducial simulations

First, we show the results of three fiducial simulations, R, F, and K in Table 1, which represent rising, flat, and decreasing rotation curves respectively. In the next section, we will show results of the other simulations in Table 1 to examine the robustness of the relation between pitch angle and the shear rate shown in this section.

4.1. Pitch angle of the mode patterns

The amplitude for several wave modes is shown for each of the fiducial simulations R, F, and K as a function of time in the top row of Fig. 4. A_m is normalised to the axisymmetric amplitude, A₀, and averaged over the radial range 4−10 kpc, which defines the region of spiral structure. The strong mode patterns are isolated by the vertical dashed lines in the top row of Fig. 4, which define the time window for the Fourier transform. The time window used is ΔT = 256 Myr. Because the top row of Fig. 4 shows that wave mode patterns appear to grow and fade on this time scale, this time window length enables the isolation of individual wave mode patterns. Although this results in limited frequency resolution, the positional information will be more reliable than that calculated from longer time windows, which may convolve multiple patterns in the Fourier analysis. However, we have confirmed that the use of longer time windows has a negligible effect on the pitch angle values.

For each of our fiducial simulations, the power spectrum of the dominant mode highlighted in the top row of Fig. 4 is calculated from the square of the amplitude given in Eq. (12), and shown as a function of radius and the pattern speed, , in the second row of Fig. 4. A wave mode pattern is eligible to be analysed if its maximum power, , is greater than 50% of the maximum power of the strongest pattern, (i.e. ): all other patterns are considered subsidiary. There are typically several patterns in each simulation that fulfil this criterion.

To demonstrate the fitting process, we focus on the most dominant patterns in each of the simulations R, F, and K. The density maps of these dominant wave mode patterns in real space polar coordinates are shown in the third row of Fig. 4. This is calculated from a sinusoidal wave of the form: ρ = A^m(R) [cos(m(θ − θ_p(R))) + sin(m(θ − θ_p(R)))] . The amplitudes and phases of each radial bin are calculated from the power spectrum in the second row of Fig. 4 using Eqs. (12) and (13) respectively. Grey scale images highlight positive (black) and negative (white) normalised density, and contours emphasise the high density regions. The plots show coherent spiral structure with well defined pitch angles where the density contrast is high.

Fig. 5 Mode pitch angles for the fiducial set of simulations, R (red circles), F (green crosses), and K (blue diamonds) as a function of shear rate.

Fig. 6 Snapshots of the disc density in polar coordinates. Density contours are overlaid in white. The traced spiral arm position is highlighted with a black line. The double peak structure at R ~ 5.5 and  ~9 kpc at snapshots t = 1.152 and t = 1.2 Gyr prevents an unambiguous fitting to a single peak, and this defines the time range in which the spiral arm can be traced.

The bottom row of Fig. 4 shows the logarithmic chi-squared fitting of the most dominant patterns in each simulation. The right side of each panel shows the normalised pattern amplitude as a function of radius, which reflects the relative strength of a pattern at a given radius. The logarithmic fitting is weighted by the amplitude shown in the right panel, and is represented by the red line (left panel). The fits are satisfactory for the radial ranges where the patterns are strong, and produce reliable pitch angles. The fitting of all other selected patterns for these simulations are very similar to those shown in the bottom row of Fig. 4. The derived pitch angles are given in Table 2.

Figure 5 shows the pitch angle dependence with shear rate (Eq. (6)). All the pitch angle values clearly show a dependence on shear rate. Simulations with higher shear rate show smaller pitch angles. This is in accordance with the qualitative trend expected of the pitch angle-shear relation from theoretical studies (e.g. Lin & Shu 1964; Julian & Toomre 1966). It is interesting to note that modes of different m and different pattern speeds in the same simulation (e.g. m = 3 and 4 in simulation F) show similar pitch angles.

4.2. Direct pitch angles of overall spiral arm features

As described in Sect. 3.2, we trace the evolution of the overall spiral arm feature directly by use of the particle density weighting method. Figure 6 demonstrates an example of the application of the arm tracing criteria to one of the spiral arms in simulation K. Because it is possible to reliably trace spiral arms which show only single peak structure for the radial range considered for fitting, we reject those snapshots that show the spiral arm with indistinct or double peak structure, which typically occurs during spiral arm formation (t = 1.152 Gyr in Fig. 6) and after the arm shows bifurcation or breaking (t = 1.2 Gyr in Fig. 6).

The results for several spiral arms in each fiducial simulation are shown in Fig. 7. It is clear that every spiral arm pitch angle decreases with time, which is consistent with winding, co-rotating spiral arms which have been reported in Wada et al. (2011), Grand et al. (2012a,b), and Baba et al. (2013). Note that this winding is also seen in the previous formalism with mode analysis, but only through a superposition of the different mode patterns: the individual mode patterns of course are defined as being formed of fixed pattern speed, , at all radii of interest. The mean of the mode pattern pitch angles calculated in the previous section is highlighted by the horizontal lines in Fig. 7. The direct pitch angle values follow the same trend with shear rate as the mode pattern pitch angles presented in Sect. 4.1, but simulations of different shear rate can overlap in direct pitch angle owing to the spread in pitch angle values produced by the winding mechanism of the spiral arm features. A snapshot of a time when direct and mode pattern pitch angles are approximately the same is shown in Fig. 8 for simulation R, F, and K. This shows the pitch angle − shear trend clearly³.

Fig. 7 Pitch angle evolution of the overall spiral arm feature for simulations R (red circles), F (green triangles), and K (blue diamonds). In all cases the pitch angle decreases with time, which indicates the winding nature of the overall density peak. The horizontal lines represent the mean mode pattern pitch angle, determined from the patterns in Fig. 4 and shown in Table 2 for simulations R (dot-dashed red), F (dashed green), and K (solid blue). Note that the range of directly measured spiral arm pitch angles clearly map out separate domains about the mode pattern pitch angles of their respective galaxies.

The winding nature of the spiral arms means that each spiral arm can exhibit several pitch angles over the spiral arm lifetime. Figure 9 shows these pitch angles plotted against galactic shear, which clearly shows that the pitch angle decreases for increasing shear rate. The range of pitch angles becomes smaller with increasing shear rate as well. This trend and scatter shown in Fig. 9 are both consistent with the pitch angle-shear rate correlation and scatter seen in real observations (e.g. Fig. 3 of Seigar et al. 2006). This may indicate that observers are seeing spiral arms at varying stages of their evolution, and therefore detect a range of pitch angles at a given rate of shear. To test the validity of these results, we explore the effect of other parameters on pitch angle in the next section.

5. Parameter survey

Up to this point, we have presented results only from the fiducial simulations R, F, and K, which clearly show the relationship between pitch angle and shear rate owing to their very different rates of shear. We now explore the effects on the pitch angle of the other parameters that vary between them.

Fig. 8 Face on view of each simulation (from left to right: simulations R, F, and K) when the directly measured spiral arm pitch angle coincides with the calculated mode pattern pitch angle. The spirals become increasingly tight going from left to right.

5.1. Resolution and softening length

We investigate the numerical robustness of the simulations by examining the effect of the number of particles and the choice of softening length. We start with simulations Fa, Fb, and Fc, which use N = 5 × 10⁶ particles with different softening lengths (see Table 1) together with the fiducial F. They are identical in every other parameter to the fiducial F simulation. The top row of Fig. 10 shows their wave mode amplitudes and dominant mode pattern phase positions. There are some differences between the higher resolution simulations, Fa, Fb, and Fc. For example, the m = 5 mode shows significant amplitude in Fc.

Because the softening length relates to the particle mass as , a direct comparison to explore the effect of resolution is between simulation F and Fb. The spiral structure grows slightly more slowly in simulation Fb (as well as the other higher resolution simulations) than in simulation F, but modes of m = 3 and 4 remain strong in all of these simulations. The difference in level of spiral structure growth for the different particle number is as expected (Fujii et al. 2011).

Fig. 9 All directly calculated spiral arm feature pitch angles plotted as a function of galactic shear for simulations R (red circles), F (green triangles), and K (blue diamonds).

The chi-squared fitting of the most dominant patterns in each simulation is shown in the bottom row of Fig. 10. The mode pattern pitch angles for all three higher resolution simulations are given in Table 2, and are all very similar to the fiducial F mode pattern pitch angles.

Figure 11 shows the pitch angles of several spiral arms that we analysed using the direct trace of the spiral arm features. Again, the arms are winding with time, and the range of pitch angles are consistent with simulation F in Fig. 7. In Fig. 11, at around t = 1.6 Gyr, simulation Fa shows a spiral arm that forms with an initial pitch angle of φ = 41 degrees, and is quickly wound. Although this initial pitch angle is high compared to that of the other arms, the later pitch angle measurements for this spiral arm overlap the range of pitch angles of all the other arms in simulations F, Fa, Fb, and Fc.

The general agreement between the mode pattern pitch angles and the range of direct pitch angles over the simulations F, Fa, Fb, and Fc indicates that the fiducial resolution of N = 1 million particles is sufficient to capture robust pitch angles. Moreover, the variation of the softening length in the assumed range does not appear to be a significant factor either, owing to the very similar mode pattern pitch angles given in Table 2 and directly measured pitch angles shown in Fig. 11.

5.2. Disc-halo mass ratio

Another variable in our simulations is the disc mass to halo mass ratio. To see whether or not this parameter affects the pitch angle, we perform the same analysis on simulations F2 and F3, which display shear rates within  ~2% of the fiducial simulation F, with lower and higher disc-halo mass ratios respectively (see Table 1). This ratio, ζ, is calculated as the ratio of the disc mass to the external mass within two radial scale lengths (as performed in D’Onghia et al. 2013). The amplitudes and density mode pattern phase positions are shown in Fig. 12. The mode pattern pitch angles calculated from the fitting in the bottom rows in Fig. 12 is presented in Table 2. The pitch angle values of F2 and F3 are similar to that of F. The directly measured pitch angles from the spiral arm feature shown in Fig. 13 also show little difference between the simulations, with perhaps the exception of the F3 spiral arm beginning t = 1 Gyr at φ ~ 40°. Overall, these results indicate that the disc to halo mass ratio does not affect the pitch angle of the spiral features, but instead the number of spiral arms, m. For example, in Fig. 12 the higher disc-mass ratio simulation, F3, displays more power in lower wave mode numbers (m = 2,3) whereas the lowest disc-halo mass ratio simulation, F2, shows the m = 7 mode to be most prominent. This is consistent with previous studies (Julian & Toomre 1966; Toomre 1981; Efstathiou et al. 1982; Carlberg & Freedman 1985; D’Onghia et al. 2013).

Fig. 10 Top row: amplitudes of the m = 1−7 wave mode numbers (colours as in top row of Fig. 4). Bottom row: phase positions of the strong mode patterns identified in top row. From left to right: simulations Fa (m = 4), Fb (m = 4), and Fc (m = 4) respectively.

We also performed simulations of intermediate shear rate values between simulations R and F with a slight alteration of disc-halo mass ratio. These simulations, labelled R2, R3, and R4 (in order from higher to lower shear), have no bulge. Figure 14 shows the direct pitch angle of several spiral arms in these simulations. While they are similar to each other, the range of pitch angles covers a slightly higher range than that of simulation F but slightly lower than that of simulation R. This agrees with the intermediate shear values shown in Fig. 2. Table 2 shows the measured pitch angle of the wave modes, which also indicates the intermediate mode pattern pitch angles between simulations R and F.

To examine the trends together, we plot the mode pattern pitch angles of simulations F, F2, F3, K, R, R2, R3, and R4 as a function of shear rate in Fig. 15. This figure shows a clear correlation between pitch angle and shear rate, which is the main finding of this paper.

The lack of effect of disc-halo mass ratio on pitch angle in combination with the difference in pitch angle between simulations F and R, which both have the same mass ratio, are convincing evidence that the shear rate is the dominant driver of pitch angle in N-body simulations of spiral galaxies.

Fig. 11 As for Fig. 7 but for simulations Fa (blue squares), Fb (red circles), and Fc (green triangles).

6. Discussion

We have shown that in N-body simulations, the measured pitch angles (measured both through the wave mode patterns and directly tracing the spiral arm features) correlates with shear rate. The range of direct pitch angles produced is in agreement with observation. We explored other simulation parameters, and show that the pitch angle is not significantly affected by the disc-halo mass ratio, resolution or softening length. One other parameter whose effect we could not explore is the stability parameter, Q, owing to the fact that it cannot be directly specified and it evolves over time (Fujii et al. 2011). Although we could not test this parameter directly, we note that the Q parameter is reported from analytical studies (e.g. Julian & Toomre 1966; Athanassoula 1984; Fuchs 2001) to have negligible effect on the pitch angle of swing-amplified patches. Also, the density wave theory of Lin & Shu (1964) does not show an explicit correlation between the pitch angle and the Q parameter. Therefore, we expect the major driver of the pitch angle value of spiral arms in N-body simulations to be the shear rate. However, this aspect still needs further study.

The observed correlation between the pitch angle of the density wave mode and galactic shear rate is qualitatively consistent with the prediction of the classic theories of both density wave theory (Lin & Shu 1964) and swing amplification theory (Julian & Toomre 1966; Toomre 1981).

In the context of swing amplification theory, spiral structure grows from density perturbations as the stellar material swings from an open to a tightly wound structure, so as to exhibit a range of inclination angles. Therefore the pitch angle may correspond to the inclination angle when each density perturbation is most amplified, around a specific inclination angle, which is correlated to shear rate (Julian & Toomre 1966).

Fig. 12 As in Fig. 10, but for simulations F (m = 3), F2 (m = 7), and F3 (m = 3).

In the context of the Lin-Shu density wave theory, each wave mode can be interpreted as a standing wave mode of constant pitch angle and fixed pattern speed. Lin & Shu (1964) demonstrate that the pitch angle of such waves is lower for higher central mass concentrations, i.e. a higher shear rate. However, there must be more than one wave mode to manifest the winding of the spiral arm, which must then be interpreted in terms of a superposition of multiple mode patterns, which changes with time (e.g. Comparetta & Quillen 2012). In this interpretation, the wave mode patterns in the inner disc region must have a faster pattern speed than that in the outer region, and must overlap at some intermediate radii. Therefore, the pitch angle begins larger than that measured for the wave mode, and then approaches the mode pitch angle while the density grows (constructive interfering). The waves then pass and move away from one another, which decreases the pitch angle further. This leads to a stretch in the azimuthal direction of the overall spiral arm density.

Fig. 13 As for Fig. 7 but for simulations F (green triangles), F2 (blue squares), and F3 (red circles).

Fig. 14 As for Fig. 7 but for simulations R2 (blue squares), R3 (green triangles), and R4 (red circles).

Fig. 15 The mode pitch angles as a function of shear for simulations, R (red circles), R2, R3, R4 (magenta plusses), F, F2, F5 (green crosses) and K (blue diamonds).

If multiple wave modes are the driving mechanism of spiral arms, the N-body simulations suggest that there are many patterns of various multiplicity, m, that are short-lived (as seen from the transient nature of the mode amplitudes in the top row of Fig. 4 for example) and recurrent. However, it is worth noting that such waves are some distance from the large scale, long timescale structures that classic spiral density wave theory was developed to produce. The formation and evolution of such wave modes should be non-linear and complicated (D’Onghia et al. 2013; Baba et al. 2013), which deserves further study, and is beyond the scope of this paper.

7. Conclusions

For the first time, to our knowledge, we have analysed the pitch angle of the spiral arm features directly and the pitch angle of the wave mode pattern in N-body simulations of disc galaxies with different galactic shear rate. The former pitch angle is derived from tracing the physical movement of the actual surface density of the spiral arms, and the latter is calculated from Fourier analysis that aims to isolate density wave mode patterns from the system that may contribute to the overall movement of the spiral arms. We presented and compared the results of both techniques, and come to the following conclusions.

• 1.

  We find that the pitch angle measured both through the wavemode analysis and direct analysis is correlated with the rate ofgalactic shear: the pitch angle is smaller for higher galactic shearrate and vice versa. This is consistent qualitatively with theanalytical predictions based on density wave theory (Lin &Shu 1964) and swing amplification theory inJulian & Toomre (1966), which we demonstrate inN-body simulations for the first time.

• 2.

  The direct pitch angles of the overall spiral arm density enhancement decrease with time, as the spiral arms grow from a relatively open arm morphology, then wind over time to become more tightly wound until they disrupt. This is consistent with previous simulations that reported winding and co-rotating spiral arms (Wada et al. 2011; Grand et al. 2012a,b; Baba et al. 2013).

• 3.

  The range of the direct pitch angles resulting from the winding spiral arm features is correlated with their shear rate: the direct pitch angle range tends to be smaller for the system with higher galactic shear and vice versa. The range of direct pitch angles at a given shear rate is similar to the scatter seen from the observed relation between the pitch angle and the shear rate in spiral galaxies reported in Seigar et al. (2006). This is consistent with the view that real galaxies exhibit transient and winding spiral arms.

Our N-body simulations demonstrate the relation between the pitch angle and the galactic shear rate. Although we explored several parameters, such as disc-total mass ratio and simulation resolution, this area of study is far from completion. We also used a fixed dark matter halo for simplicity, and left out the gas component. In real galaxies, there are also constant minor mergers and tidal interactions with satellite galaxies, which we have not explored. However, we suggest that this study highlights the relation between pitch angle and the galactic shear rate, and encourages further studies with more realistic and complicated models. If this relation is a dominant mechanism to determine the pitch angle of the spiral arms, because the late type spiral galaxies tend to have rising rotation curves, this relation will become key to explain the correlation between the pitch angle and the Hubble type (Hubble 1926; Kennicutt 1981).

Acknowledgments

The authors thank Jerry Sellwood for suggestions regarding the mode pattern analysis. The authors acknowledge the support of the UK’s Science & Technology Facilities Council (STFC Grant ST/H00260X/1). The calculations for this paper were performed on Cray XT4 at Centre for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan and the DiRAC Facilities, Legion and COSMOS, jointly funded by STFC and the Large Facilities Capital Fund of BIS. The authors acknowledge support of the STFC-funded Miracle and COSMOS Consortium (part of the DiRAC facility) in providing access to the UCL Legion and Cambridge COSMOS High Performance Computing Facilities. The authors additionally acknowledge the support of UCL’s Research Computing team with the use of the

Legion facility. This work was carried out, in part, through the Gaia Research for European Astronomy Training (GREAT-ITN) network. The research leading to these results has received funding from the European Union Seventh Framework Programme ([FP7/2007-2013] under grant agreement number 264895.

References

All Tables

All Figures

	Fig. 1 Circular velocity at t = 0 for simulation R (very thick dashed red), R2 (thick dashed blue), R3 (medium dashed green), R4 (thin dashed cyan), F (thick dot-dashed green), F2 (thin dot-dashed red), F3 (medium dot-dashed blue), and K (solid blue).
In the text

	Fig. 2 Galactic shear rate, Γ, for all simulations. Colours are the same as Fig. 1. Note the reduced radial range compared to Fig. 1.
In the text

	Fig. 3 Toomre stability parameter, Q, at t = 0 for all simulations. Colours are the same as Fig. 1.
In the text

Fig. 4 Top row: amplitudes calculated from Eq. (10) and averaged over a radial range of 4−10 kpc of spiral modes m = 1 (thick red); 2 (thick green); 3 (thick blue); 4 (thick yellow); 5 (thin red); 6 (thin green); 7 (thin blue), and 8 (thin yellow) normalised to the axisymmetric m = 0 mode, as a function of time for simulations R (left), F (middle), and K (right). Vertical dashed lines represent the time window of a Fourier transform applied in Sect. 3.1. Second row: power spectra calculated from Eq. (12) of simulation R for the m = 8 mode (left), F for the m = 3 mode (middle), and K for the m = 2 mode (right). Prominent ridges (dark pixels) span between 4−10 kpc in most cases. Third row: in polar coordinates, the density map of the dominant density wave mode pattern selected from rows of , 30, and 24 km s-1 kpc-1 for simulations R, F, and K respectively. White regions indicate areas of low density and black regions indicate areas of high density. Contours emphasis the highest density regions. Bottom panels: dominant mode pattern positions (black points) calculated from Eq. (13) in the azimuth-radius plane for the corresponding patterns in the row above. The red lines show the lines of best fit for each pattern. The right side of each panel shows the radial amplitude profile, which is used to weight the fitting.

In the text

	Fig. 5 Mode pitch angles for the fiducial set of simulations, R (red circles), F (green crosses), and K (blue diamonds) as a function of shear rate.
In the text

	Fig. 6 Snapshots of the disc density in polar coordinates. Density contours are overlaid in white. The traced spiral arm position is highlighted with a black line. The double peak structure at R ~ 5.5 and  ~9 kpc at snapshots t = 1.152 and t = 1.2 Gyr prevents an unambiguous fitting to a single peak, and this defines the time range in which the spiral arm can be traced.
In the text

Fig. 7 Pitch angle evolution of the overall spiral arm feature for simulations R (red circles), F (green triangles), and K (blue diamonds). In all cases the pitch angle decreases with time, which indicates the winding nature of the overall density peak. The horizontal lines represent the mean mode pattern pitch angle, determined from the patterns in Fig. 4 and shown in Table 2 for simulations R (dot-dashed red), F (dashed green), and K (solid blue). Note that the range of directly measured spiral arm pitch angles clearly map out separate domains about the mode pattern pitch angles of their respective galaxies.

In the text

	Fig. 8 Face on view of each simulation (from left to right: simulations R, F, and K) when the directly measured spiral arm pitch angle coincides with the calculated mode pattern pitch angle. The spirals become increasingly tight going from left to right.
In the text

	Fig. 9 All directly calculated spiral arm feature pitch angles plotted as a function of galactic shear for simulations R (red circles), F (green triangles), and K (blue diamonds).
In the text

	Fig. 10 Top row: amplitudes of the m = 1−7 wave mode numbers (colours as in top row of Fig. 4). Bottom row: phase positions of the strong mode patterns identified in top row. From left to right: simulations Fa (m = 4), Fb (m = 4), and Fc (m = 4) respectively.
In the text

	Fig. 11 As for Fig. 7 but for simulations Fa (blue squares), Fb (red circles), and Fc (green triangles).
In the text

	Fig. 12 As in Fig. 10, but for simulations F (m = 3), F2 (m = 7), and F3 (m = 3).
In the text

	Fig. 13 As for Fig. 7 but for simulations F (green triangles), F2 (blue squares), and F3 (red circles).
In the text

	Fig. 14 As for Fig. 7 but for simulations R2 (blue squares), R3 (green triangles), and R4 (red circles).
In the text

	Fig. 15 The mode pitch angles as a function of shear for simulations, R (red circles), R2, R3, R4 (magenta plusses), F, F2, F5 (green crosses) and K (blue diamonds).
In the text