# Classifying and modelling spiral structures in hydrodynamic simulations of astrophysical discs

Classifying and modelling spiral structures in hydrodynamic simulations of astrophysical discs Abstract We demonstrate numerical techniques for automatic identification of individual spiral arms in hydrodynamic simulations of astrophysical discs. Building on our earlier work, which used tensor classification to identify regions that were ‘spiral-like’, we can now obtain fits to spirals for individual arm elements. We show this process can even detect spirals in relatively flocculent spiral patterns, but the resulting fits to logarithmic ‘grand-design’ spirals are less robust. Our methods not only permit the estimation of pitch angles, but also direct measurements of the spiral arm width and pattern speed. In principle, our techniques will allow the tracking of material as it passes through an arm. Our demonstration uses smoothed particle hydrodynamics simulations, but we stress that the method is suitable for any finite-element hydrodynamics system. We anticipate our techniques will be essential to studies of star formation in disc galaxies, and attempts to find the origin of recently observed spiral structure in protostellar discs. methods: numerical, stars: formation, ISM: kinematics and dynamics 1 INTRODUCTION Astrophysical disc structures are found across a wide range of scales – from disc galaxies to discs surrounding active galactic nuclei (AGNs), to discs around protostars and even around protoplanets. Almost as commonly, these discs can be perturbed into producing spiral structures. The means by which spiral structures are produced depends on the object being studied. Broadly, spiral structure is caused by some form of unstable non-axisymmetric perturbation. How this perturbation subsequently evolves into spiral structure depends on its origin, but in most cases is caused by the disc being susceptible to gravitational instability, where (for the gas) the Toomre parameter is   $$Q = \frac{c_{\rm s} \kappa }{\pi {\rm G} \Sigma } \lessapprox 1,$$ (1)where cs is the sound speed of the gas, κ is the epicyclic frequency, G is the gravitational constant, and Σ is the local disc surface density. A similar expression exists for the stability of a stellar disc (where the radial stellar velocity dispersion replaces cs). Spiral structures have been observed in galaxies for some 150 yr, going back to Rosse (1850)’s observations of M51. We now know that spiral galaxies constitute over half of observed massive galaxies (Lintott et al. 2011). This has provided theorists with significant data reserves on which to test spiral generation theories. Spiral structures observed in disc galaxies can be self-generated, either by local instabilities/noise which are swing amplified to generate arms (e.g. Sellwood & Carlberg 1984), or through globally propagating (quasi-stationary) density waves, as predicted by the pre-eminent Lin–Shu density wave theory (Lin & Shu 1964). Spiral structures can also be externally generated through tidal interactions, e.g. during galaxy mergers [Holmberg (1941)’s experiments in this area are particularly illuminating]. Indeed, all three mechanisms (instabilities, density waves, and interactions) can and do collaborate to produce the spiral structures observed both in observations and simulations (see Dobbs & Baba 2014, for a detailed review). Spiral arms drive important physical processes in galaxies. Molecular clouds and recent star formation tend to be concentrated in and near spiral arms (Schinnerer et al. 2013; Heyer & Dame 2015; Ragan et al. 2016; Schinnerer et al. 2017), most likely due to the compression and shocking of gas as it falls into the arm potential (Bonnell et al. 2006; Dobbs & Bonnell 2007; Bonnell, Dobbs & Smith 2013). Spiral morphology is also directly correlated with the surface density of neutral atomic hydrogen in the disc, and the central stellar bulge mass (Davis et al. 2015). The gravitational torques induced by spiral structures will drive material into the inner few kpc of galaxies. The potential asymmetries caused by stellar bars can then deliver this material further inwards, feeding AGNs (see e.g. Alexander & Hickox 2012; Querejeta et al. 2016). The wide-ranging effects of spirals on galaxy evolution have driven a great deal of effort on determining their properties in observations. Attempts to characterize observed spiral structure in galaxies began with visual classification of the total number of arms, and their ‘openness’ (Hubble 1926), which can be quantified by the pitch angle ϕ. The modern Hubble sequence divides spiral galaxies into barred (SB) and unbarred (S) spirals, with a further subclass (a–d) denoting the tightness of the spiral winding. Elmegreen (1990) proposed dividing spirals based on their arm number (grand design, multi-armed, and flocculent). More quantitative measurements of observed spiral galaxies involve Fourier analysis of the image (Rix & Zaritsky 1995; Foyle et al. 2010), which decomposes the azimuthal variations in surface brightness into a Fourier sum over m spiral modes with amplitude Am, e.g.:   $$\frac{\mu (R, \theta )}{\bar{\mu }(R)} = \sum _{m=1}^{\infty } A_m (R) {\rm e}^{{\rm i}m(\theta -\theta _m)}.$$ (2)Pitch angles can be fitted to images assuming a logarithmic spiral of constant ϕ, a prediction of density wave theory (Kennicutt 1981), or via hyperbolic spirals with radially varying pitch angle (e.g. Seiden & Gerola 1979). The pattern speed of an arm can be investigated by studying star formation (and tracers of star formation) inside the arm (Egusa, Sofue & Nakanishi 2004; Egusa et al. 2009). A steady arm of constant pattern speed will set up a sequence of tracers – e.g. HI, CO, 24 $$\mu$$m emission from enshrouded stars, and UV emission from unobscured stars – from upstream to downstream of the arm's corotation. The failure to see this sequence in local spiral galaxies is strong evidence that spiral structures do not persist beyond the local dynamical time (Foyle et al. 2011). Other observational techniques can also provide convincing evidence for the origin of spiral structures. For example, Choi et al. (2015) resolve the stellar populations along the north-east arm of M81, to show it does not possess a constant pattern speed. This points to a kinematic origin, driven by tidal interactions with nearby galaxies M82 and NGC 3077s. Protostellar discs can also self-generate spiral structure through local instabilities. Shortly after the formation of a protostellar system, the protostellar disc has a mass comparable to that of the central star. Such discs can soon assume a marginally unstable state where Q ∼ 1, and the disc mass determines the spiral modes present (Lodato & Rice 2005; Forgan et al. 2011). Interactions with a companion either internal to the disc (such as a protoplanet) or external to it (such as a close stellar encounter) also generate tidally induced spiral arms. Spiral arms in protostellar discs are also efficient outward transporters of angular momentum, driving rapid accretion and assembly of protostars (Laughlin & Rozyczka 1996). Spiral structures with a sufficient surface density contrast can concentrate dust grains (Rice et al. 2004; Clarke & Lodato 2009; Booth & Clarke 2016), promoting grain growth and setting the initial conditions for planet formation via core accretion, as well as altering local chemistry (Ilee et al. 2011; Evans et al. 2015; Ilee et al. 2017). In extremis, spiral arms of sufficiently large density amplitude can induce protostellar discs to fragment into bound objects (Rice, Lodato & Armitage 2005; Forgan & Rice 2011; Tobin et al. 2016), providing an alternate formation channel for gas giants and substellar objects at large orbital semimajor axis (Forgan & Rice 2013; Forgan, Parker & Rice 2015; Vigan et al. 2017). For sufficiently massive protostellar systems, fragmentation can also generate binary star systems (Bonnell & Bate 1994; Bate, Bonnell & Bromm 2003; Tobin et al. 2016). Spiral structure has only recently been observed in protoplanetary discs, initially in near-infrared observations of scattered light, which is typically most sensitive to the upper disc surface. In particular, extended two-armed structures have been detected around SAO206462 (Muto et al. 2012) and MWC 758 (Benisty et al. 2015). These arms are detectable out to relatively large distances (up to 100 au) from the parent star, with pitch angles of the order of 10°. Recent ALMA observations have shown discs with spiral structure that extends down to the disc midplane, for example the recent detection of m = 2 spiral structure around Elias 2-27, with a measured pitch angle of 7.9±0.4° (Pérez et al. 2016) – the origin of this structure is not yet clear (Meru et al. 2017). Characterizing these spirals, and determining their origin, yields crucial information about protostellar accretion and the protostar's approach to the main sequence, as well as the formation of planetary systems. If we are to identify the origin of newly observed protostellar disc spiral structures, or to study how spiral arms govern and drive star formation in galaxies, analysing spiral structure driven in hydrodynamic simulations is crucial. As such, we require tools to identify and characterize spiral arms in these simulations. Regardless of scale, characterizing spiral morphology yields important diagnostics of what is driving the spiral structure. In particular, the number of arms, their amplitude and pitch angles are sensitive to both the driving mechanism and the disc's properties. For example, the spiral wake produced by planets embedded in a protostellar disc adopts a pitch angle which is a function of the disc temperature and rotation profile, as well as the planetary mass (Rafikov 2002; see also Pohl et al. 2015; Zhu et al. 2015). Gravitationally unstable discs also produce spiral structure with pitch angles and arm number that depend in particular on the disc mass (Dong et al. 2015), which may be a good deal larger than the observed disc mass (Forgan et al. 2016b). What is clear is that both theoretical and observational astrophysicists stand to gain a great deal from higher quality characterization of spiral structure in numerical simulations of astrophysical discs. Theorists gain important insights into how spiral structures are generated, and how they affect the thermal and chemical history of gas and dust, and the future evolution of the disc. Observers gain diagnostics for what is driving the spiral structure in their observations, and glean information on disc properties that is generally orthogonal to other methods. Most attempts to characterize spiral structure in simulations rely on Fourier decomposition (e.g. Cossins, Lodato & Clarke 2009; Dobbs et al. 2010; Forgan et al. 2011; Mata-Chavez, Gomez & Puerari 2014; Pettitt, Tasker & Wadsley 2016). This gives important information on the relative strengths of the spiral modes at play, and the pitch angle and pattern speed of the dominant mode. However, it does not give the pitch angle and pattern speed of individual arms. It also does not inform us as to what sections of the disc are currently in the spiral (and which sections are in the inter-arm regions). Recent attempts at spiral arm characterization have moved away from this Fourier analysis. For example, Grand, Kawata & Cropper (2012) attempted to directly identify individual stellar spiral arms in N-body/hydrodynamic simulations of a barred spiral galaxy, by finding the location of a series of density peaks in a range of annuli. However, this does not directly provide data on the arm width, or which fluid elements currently reside in the spiral. In this paper, we demonstrate that judicious use of tensor classification on hydrodynamic simulation data (Forgan et al. 2016a) allows the identification of fluid elements that are inside spiral structures (or in the inter-arm regions). Further analysis allows the isolation of individual arms to obtain their shape parameters, as well as the pattern speed of the wave. Our examples focus on smoothed particle hydrodynamics (SPH) simulations, but we emphasize that our methods only require the ability to compute derivatives, and are therefore applicable to any hydrodynamic simulation. 2 METHODS Our spiral detection algorithm has two distinct parts. In the first part, fluid elements undergo tensor classification (Forgan et al. 2016a) to determine whether their behaviour indicates they are in fact inside a spiral structure. In the second part, the fluid elements identified as spirals are extracted from the main simulation, and a friends-of-friends algorithm is used on this population to identify the spine of each individual spiral. We describe these procedures below. Our code is published on Github at https://github.com/dh4gan/tache. 2.1 Tensor classification We follow the same procedure as described in Forgan et al. (2016a), which itself builds on work originally applied to N-body simulations of the cosmic Web (see e.g. Hahn et al. 2007; Forero-Romero et al. 2009). In our formalism, tensor classification determines the topology of a chosen field at the location of a given SPH particle. We consider the topology of either the gravitational potential Φ, by computing the tidal tensor Tij as follows:   $$T_{ij} = \frac{\mathrm{\partial} ^2 \Phi }{ \mathrm{\partial} x_i \mathrm{\partial} x_j},$$ (3)or the velocity field via the velocity shear tensor σij as follows:   $$\sigma _{ij} = -\frac{1}{2}\left(\frac{\mathrm{\partial} v_i}{\mathrm{\partial} x_j} + \frac{\mathrm{\partial} v_j}{\mathrm{\partial} x_i} \right).$$ (4)The classification proceeds as follows. Once the tensor to be used (T or σ) is selected, its eigenvalues λi and their corresponding eigenvectors $$\mathbf {n_i}$$ are computed for every particle, e.g.:   $$T \mathbf {n}_j = \lambda _j \mathbf {n}_j.$$ (5)The eigenvalues are defined so λ1 ≥ λ2 ≥ λ3. As T and σ are real and symmetric, these eigenvalues are always real. Both tensors assume that the fields being investigated (potential, velocity) are smooth and continuous, so that the derivatives are always defined. We then compute E, which is defined as the number of positive eigenvalues.1 The value of E determines a topological classification for each fluid element: E = 0 → ‘void’ (0-D manifold) E = 1 → ‘sheet’ (1-D manifold) E = 2 → ‘filament’ (2-D manifold) E = 3 → ‘cluster’ (3-D manifold). We can see this by considering the tidal tensor in the context of Zeldovich theory (Zel'dovich 1970). A test particle orbiting a local extremum ∇Φ = 0 has the following (linearized) equation of motion:   $$\ddot{x}_{i} = -T_{ij} x_j.$$ (6)If we operate in a basis where T is diagonal, then the stability of the orbit is determined by the sign of the eigenvalues λi. A single positive eigenvalue ensures a stable orbit along the corresponding axis (or a sheet-like structure, if multiple particles are present). Two positive eigenvalues allow an extra degree of freedom (filaments) and finally all three being positive allow orbits of any degree (clusters). If no eigenvalues are positive, stable orbits cannot be achieved (voids). Tensor classification using the tidal tensor therefore determines the manifold dimension of the local potential, and as such the preferred structure for matter to collapse into if only gravity is present. Classification using the velocity shear tensor instead diagnoses the manifold being sculpted by the flow at any given instant. As a result, the finite elements of any hydrodynamic simulation can be grouped into four components, with each component composed of fluid elements of a matching E. Discs are inherently sheet-like structures, and hence perturbations from axisymmetry are classified as either filaments or clusters. We will see that depending on the strength of spiral structures, fluid elements classified as either filaments or clusters trace spiral structures, and that fluid elements classified as sheets will trace unperturbed disc material in the inter-arm regions. As a general rule, the tidal tensor is better suited to tracing high-amplitude spiral structures that induce a strong surface density perturbation, but the velocity shear tensor is better at revealing low-amplitude spiral structures that produce weak surface density perturbations.2 Tensor classification also yields eigenvectors corresponding to each eigenvalue. In sheets, the eigenvector $$\mathbf {n}_1$$ corresponding to the single positive eigenvalue λ1 defines the normal or symmetry axis of the sheet. In filaments, the eigenvector $$\mathbf {n}_3$$ corresponding to the single negative eigenvalue λ3 provides the flow direction of the filament. We do not use eigenvectors in this analysis, but note that they may also be of interest in, e.g. determining the relative geometry of spiral structures to their host disc. 2.2 Identifying individual spirals Once tensor classification is complete, we extract all particles that correspond to either filament or cluster classification, as these trace the spiral structure in the disc, and discard the other particles. To determine the spines of each spiral arm on this subset, we run a friends-of-friends algorithm as follows. We specify a minimum spherical radius from the centre of the disc, rmin, from which to begin our study. As we are interested in spines, we will also select the top x% percentile in density of the remaining particles. We then select particle i with the largest density, where radius ri > rmin. This particle forms the first component of the spiral's spine. A sphere of radius L is then drawn from ri. The particles inside the sphere are tested, and particle j is selected if it is the densest particle with rj > ri. Particle j then forms the second component of the spiral's spine, j is then set to i, and the algorithm is repeated until a particle can no longer be found that meets the above conditions. Whenever a particle is tested, it is removed from the list of particles and is not tested again, regardless of whether it forms a component of the spiral spine or otherwise. The act of testing assigns the particle to that particular spiral. As a consequence, any particle within a distance L of any spiral spine point belongs to that spiral. If the conditions can no longer be met, a new spiral is begun at the location of the densest particle not yet tested, and the procedure begins again. Note that throughout we are implicitly assuming that the particle density decreases with increasing radius. Three parameters form the identification algorithm: the linking length L, the minimum radius from which to begin, rmin, and the percentile of the population from which to conduct the analysis, x%. The values needed for these three parameters will depend on the resolution and input physics of the simulation being analysed. The linking length L should clearly be smaller than the spiral arm spacing to avoid the algorithm ‘jumping’ between individual arms, while still being larger than the typical smoothing length to obtain a sufficient number of neighbouring particles. The minimum radius should be selected based on prior knowledge of where spiral structure is well resolved (as well as the locations from which the user desires to measure the said structure). The percentile selection depends on how well defined the structures are post-tensor classification. Selecting a low x% will result in easier classification, but will prevent the classification of weak structures. This algorithm delivers Nspiral sets of coordinate points, with each set denoting the spine of a spiral structure. We elect to fit these points via Markov Chain Monte Carlo (MCMC) to a logarithmic spiral   $$r = r_0 + a e^{b\theta },$$ (7)with fit parameters $$(r_0= \sqrt{x_0^2 + y^2_0}, a,b)$$. The origin of the spiral is defined at θ = 0; r = r0 + a. A pure logarithmic spiral will have a constant pitch angle ϕ, which in radians is simply   $$\phi = \arctan b.$$ (8) It is worth noting that if the spiral arms exclusively constitute the densest regions of our simulation, then it is likely that tensor classification will not be necessary to determine the location of the spiral spines, and that the above friends-of-friends algorithm will work perfectly well. However, if there are other dense structures in the simulation (such as condensing clumps), then the algorithm will attempt to erroneously fit spiral spines to them. We have run test calculations on fragmenting discs to show that tensor classification can successfully remove dense structures that are not spiral arms, preventing this problem. Also, without tensor classification, only the spine of the arm can be determined – the width of the spiral arms or the properties of the gas contained within the arm cannot be measured. 3 TESTS We now test our algorithm on three examples of spiral structure in astrophysical discs. The origins of the structure in each case are subtly different. In the first, we consider a self-gravitating protostellar disc, whose spiral structure is entirely due to the self-gravity of the disc gas. In the second case, we consider a spiral galaxy with an imposed fixed potential that induces well-defined structures. In the last, we consider a galaxy with a ‘live potential’ of star particles, with a resulting spiral structure that is highly flocculent. We use SPH for all three tests. SPH is a Lagrangian hydrodynamics method, where the fluid is discretized into individual particles. Each SPH particle i possesses position and velocity vectors $$\mathbf {r}_i$$, $$\mathbf {v}_i$$, internal energy ui, and smoothing length hi, which is adjusted such that a sphere of radius 2hi encompasses a suitable number N of neighbour SPH particles. The local fluid density is established by performing a sum over the neighbour particles’ smoothing kernels   $$\rho (\mathbf {r}) = \sum ^{N}_{i=1} m_i W(\left|\mathbf {r} - \mathbf {r}_i\right|, h).$$ (9)Writing down the system's Lagrangian and applying variational methods (in combination with equation 9) produce an algorithm for solving the equations of hydrodynamics on the particle distribution (with appropriate compensation for shocks and fluid mixing; see Monaghan 1992, 2005; Price 2012, for reviews). Table 1 shows the values of L, rmin, and x% used in the three tests run in this paper. Table 1. The fit parameters used in the spiral detection algorithm for the test cases evaluated in this paper. Section  Simulation  L  rmin  x%  3.1  Protostellar disc  6 au  10 au  20  3.2  Galaxy, analytic potential  0.2 kpc  1 kpc  30  3.3  Galaxy, live potential  0.2 kpc  1 kpc  30  Section  Simulation  L  rmin  x%  3.1  Protostellar disc  6 au  10 au  20  3.2  Galaxy, analytic potential  0.2 kpc  1 kpc  30  3.3  Galaxy, live potential  0.2 kpc  1 kpc  30  View Large 3.1 A self-gravitating protostellar disc For our first test, we return to a simulation previously analysed in Forgan et al. (2016a). The top left-hand panel of Fig. 1 shows an SPH simulation of a self-gravitating disc of mass 0.25 M⊙, orbiting a sink particle representing a star of 1 M⊙. The initial disc surface density profile follows Σ ∝ r−1, and has a maximum radius of 50 au. The sound speed profile fixes the initial Toomre parameter Q = 2 at all radii. Figure 1. View largeDownload slide Spiral arm detection in a self-gravitating protostellar disc. Top left: a snapshot from the full SPH simulation. Top right: the same snapshot, with only the SPH particles identified as exhibiting arm-like behaviour (through classification of the tidal tensor). Bottom left: the spine points of each arm (green crosses) accompanied by the best-fitting logarithmic spirals in red. Bottom right: The posterior distribution of the logarithmic spiral parameters for the uppermost spiral. Figure 1. View largeDownload slide Spiral arm detection in a self-gravitating protostellar disc. Top left: a snapshot from the full SPH simulation. Top right: the same snapshot, with only the SPH particles identified as exhibiting arm-like behaviour (through classification of the tidal tensor). Bottom left: the spine points of each arm (green crosses) accompanied by the best-fitting logarithmic spirals in red. Bottom right: The posterior distribution of the logarithmic spiral parameters for the uppermost spiral. The radiative transfer formalism of Forgan et al. (2009) is used, implementing realistic local cooling alongside typical compressive and shock heating to settle into a self-regulated, marginally stable state, where the Toomre parameter Q ∼ 1.5. Marginally stable self-gravitating discs are able to sustain quasi-steady spiral structure for as long as the Toomre parameter can be held at this value, which depends critically on the disc's total mass, and its ability to cool. The top right-hand panel of Fig. 1 shows the result of tensor classification using the tidal tensor. Given the deep gravitational perturbations caused by this spiral structure, we elect to use E = 3 (cluster) particles to define the minimum potential of the structure. The bottom row of Fig. 1 shows the result of spiral arm identification. Bottom left shows the identified spine points of each spiral (green points), with accompanying logarithmic fits via MCMC (red curves). An example of the derived posteriors from the MCMC fits is shown in the bottom right-hand panel (for the uppermost spiral). We see that logarithmic spirals can produce good-quality fits to the data, especially in the outer disc. The uppermost spiral is well fitted by a logarithmic spiral with (a, b) = (3.634 ± 0.300, 0.2468 ± 0.0123), giving a pitch angle of 13.8° ± 0.6°. The posterior distributions for a and b are well behaved and unimodal, with only a weak degeneracy between a and b. The inner radii of the spiral shows some deviation from the fit, suggesting that pitch angles may not be constant in this region. This is reflected in the best-fitting centre of the spiral deviating by around 1–2 au from the star's actual position at the origin. Repeating this analysis over subsequent 40 time-steps of the simulation (around 1.5 outer rotation periods, Fig. 2) shows that the matched spirals tend to share relatively similar parameters over this interval. The sample mean pitch angle derived over this duration is 13.02°, with a sample standard deviation of around 3.4°. Interestingly, the scale parameter a tends to assume slightly larger values than our example spiral, with a mean of 4.996 au, and a standard deviation of 1.02 au. Figure 2. View largeDownload slide The distribution of pitch angle ϕ (left) and a over the course of 40 snapshots of the self-gravitating disc simulation (approximately 500 yr of evolution). Figure 2. View largeDownload slide The distribution of pitch angle ϕ (left) and a over the course of 40 snapshots of the self-gravitating disc simulation (approximately 500 yr of evolution). Being able to measure the spiral structure over more than one time-step allows us to calculate the pattern speed of a given arm. Given that for a logarithmic spiral   $$\theta = \frac{1}{b} \ln (r/a),$$ (10)we can estimate the pattern speed $$\Omega _{\rm p} = \dot{\theta }$$ given $$\dot{b}$$ and $$\dot{a}$$:   $$\dot{\theta } = \frac{-\dot{b}}{b} \theta -\frac{\dot{a}}{ab}.$$ (11)Note that $$\dot{\theta }$$ has a dependence on θ only if the pitch angle is time-dependent. Having computed $$\dot{a}$$ and $$\dot{b}$$ from the two fits, we can evaluate the pattern speed over the range θ = [0, 2π] and take an average. For our example spiral, we evaluate (a, b) for two snapshots and compute $$\Omega _{\rm p} = 3.5 \times 10^{-9}\, \rm {rad \, s^{-1}}$$. The corotation radius, where Ω = Ωp, is approximately 14.6 au (Fig. 3). Figure 3. View largeDownload slide The ratio of the pattern speed to the bulk angular velocity Ωp/Ω as a function of radius. Corotation occurs when this ratio is unity, which occurs at 14.6 au, interior to which the disc is dominated by artificial viscosity. Figure 3. View largeDownload slide The ratio of the pattern speed to the bulk angular velocity Ωp/Ω as a function of radius. Corotation occurs when this ratio is unity, which occurs at 14.6 au, interior to which the disc is dominated by artificial viscosity. Artificial viscosity dominates the disc interior to this radius (see e.g. Artymowicz & Lubow 1994; Murray 1996; Lodato & Rice 2004; Lodato & Price 2010; Forgan et al. 2011). The Toomre Q parameter also reaches the instability regime just beyond this radius (∼20 au). As a result, the total disc stress has an artificial minimum at corotation, rather than continuing to decrease with decreasing radius (Rice & Armitage 2009). This shows how the inner disc resolution affects the launching point of spiral structures. Indeed, unstable non-axisymmetric normal modes in self-gravitating gaseous discs are expected to trail the flow at all radii (Papaloizou & Savonije 1991). 3.2 A disc galaxy with fixed arm potential In our second test, we use a galaxy model with analytically defined spiral structure imposed by an external gravitational potential. The galactic potential is represented by a combination of an axisymmetric term plus a perturbation of the spiral arms   $$\Phi = \Phi _{0} + \Phi _{{\rm pert}}.$$ (12) The first component is given by a logarithmic potential (e.g. Binney & Tremaine 2008)   $$\Phi _0 = \frac{1}{2} v^2_0 \left(R^2 + R^2_{\rm c} + (z/z_{\rm q})^2\right),$$ (13)which has a rotation curve given by   $$v_{\rm c}(R) = v_0 \frac{R}{\sqrt{R_{\rm c}^2 + R^2}},$$ (14)where v0 is a velocity parameter, Rc is a characteristic radius, and zq is a vertical scale factor. This produces a flat rotation curve at values larger than Rc. For our simulations, $$v_0 = 120 \, \rm {km \,s^{-1}}$$, Rc = 2.5 kpc, and zq = 0.7. The velocity parameter corresponds to that of the rotation curve of the local spiral galaxy M33 (Corbelli & Salucci 2000; Sie Kam et al. 2015). The spiral arm perturbation uses the scheme of Cox & Gómez (2002)   $$\Phi _{{\rm pert}}(R, \theta ,z, t) = \sum _n A_n(R,z) \cos \left(n \Gamma (R, \theta ,t) \right),$$ (15)where An describes the perturbation amplitude which decays with both R and z (see Cox & Gómez 2002, for details). Γ determines the form of spiral generated by the perturbation   $$\Gamma (R, \theta ,t) = N \left(\theta + \Omega _{\rm p} t - \frac{\ln (R/R_0)}{\tan \phi } - \theta _{\rm p} \right).$$ (16)The number of spiral arms N = 4, and their pattern speed is $$\Omega _{\rm p} = 23 \, \rm {km\,s^{-1}\,kpc^{-1}}$$ (the pitch angle ϕ is 15°). The constant θp defines the location of the arm at a fiducial radius R0. The gaseous component is assumed to have an exponential surface density profile: $$\Sigma (R) = \Sigma _0 e^{-R/R_{\rm d}}$$, where Σ0 is the central surface density and Rd is the scale radius. The gas self-gravity is not included. This analytic potential model is advantageous for its reduced computational expense, and also allows us to have a well-defined spiral structure for testing the spiral detection algorithm. Fig. 4 shows the full SPH simulation (top left), the particles identified as ‘spiral-like’ using the velocity shear tensor (top right), and the particles identified as belonging to the inter-arm regions (E = 1, ‘sheet’ classification, bottom left). The spines of the four main spirals are easily identified (bottom right of Fig. 4), and are well fitted by logarithmic spirals each with b ≈ 0.25 ± 5 per cent (i.e. ϕ ≈ 14° ± 0.7°), representing the input perturbation model well. Figure 4. View largeDownload slide The evolution of a galactic disc under an analytic spiral potential. Top left: the full SPH simulation containing 2 m gas particles. Top right: the particles identified as ‘spiral-like’ (filament class). Bottom left: particles identified as moving in a disc configuration (sheet class). Bottom right: Fits to the spiral structure. Figure 4. View largeDownload slide The evolution of a galactic disc under an analytic spiral potential. Top left: the full SPH simulation containing 2 m gas particles. Top right: the particles identified as ‘spiral-like’ (filament class). Bottom left: particles identified as moving in a disc configuration (sheet class). Bottom right: Fits to the spiral structure. The separation of arm and inter-arm particles in azimuth is clear, if we bin particles by class in an annulus of 0.1 kpc (Fig. 5). We can in fact apply Gaussian fits to the four arms to determine their angular width, resulting in a full width half-maximum (FWHM) of 0.2825 rad. This corresponds to a rather large distance of around 0.5 kpc at this radius, but this is principally due to the winding of the arms aligning them with the azimuthal vector. Figure 5. View largeDownload slide The number of particles classified as ‘arm’ versus ‘inter-arm’ in the galactic disc under the analytic potential, as a function of θ, in an annulus with inner and outer radii R = 1.95, 2.05 kpc. Also plotted are Gaussian fits to the arm histograms, with each Gaussian possessing a full width half-maximum (FWHM) of 0.2825 rad. The galaxy rotates in the positive θ direction. Figure 5. View largeDownload slide The number of particles classified as ‘arm’ versus ‘inter-arm’ in the galactic disc under the analytic potential, as a function of θ, in an annulus with inner and outer radii R = 1.95, 2.05 kpc. Also plotted are Gaussian fits to the arm histograms, with each Gaussian possessing a full width half-maximum (FWHM) of 0.2825 rad. The galaxy rotates in the positive θ direction. We can compute the width of a spiral arm along its entire extent by interrogating the particles that have been identified as belonging to it. We do this for the rightmost spiral in the bottom right-hand panel of Fig. 4, and compute each particle's minimum distance to the spiral spine. The left-hand panel of Fig. 6 shows a 2D histogram of the data (we overplot percentiles of the bin counts for convenience). We also extract 1D histograms by binning in radius (the right-hand panel of Fig. 6) to show that the spiral arm width is essentially constant with radius (with the FWHM of the 1D histograms being approximately 0.1 kpc). Figure 6. View largeDownload slide The width of an individual spiral arm in the galaxy model with analytic potential, as a function of radius. Left: For a population of particles identified by our algorithms as belonging to a specific arm, we compute their minimum distance from the arm and bin them in the above 2D histogram. We overplot percentiles of the bin counts to highlight how the typical separation from an arm is constant with radius. Right: 1D histograms of distance from the arm for the same population, for all radii, and inner and outer radii. Figure 6. View largeDownload slide The width of an individual spiral arm in the galaxy model with analytic potential, as a function of radius. Left: For a population of particles identified by our algorithms as belonging to a specific arm, we compute their minimum distance from the arm and bin them in the above 2D histogram. We overplot percentiles of the bin counts to highlight how the typical separation from an arm is constant with radius. Right: 1D histograms of distance from the arm for the same population, for all radii, and inner and outer radii. As for the protostellar disc, we can compute pattern speeds by comparing two snapshots of the simulation. We find pattern speeds of around $$20.4\, \rm {km\,s^{-1}\,kpc^{-1}}$$ (6.6 × 10−16rad s−1), with a 1σ uncertainty of around 19 per cent. Given the rotation curve, these parameters define the corotation radius as   $$R_{\rm co} = \sqrt{\left(\frac{v_0}{\Omega _{\rm p}}\right)^2 - R^2_{\rm c}} = 5.33 \, \rm {kpc}.$$ (17)The analytic spiral arm potential has a pattern speed of $$23\, \rm {km\,s^{-1}\,kpc^{-1}}$$ (7.45 × 10−16rad s−1), and corotation at 4.9 kpc. Hence, the analytic value resides within the 1σ uncertainties of our calculation. 3.3 A disc galaxy with a live stellar potential We also develop a N-body version of the model galaxy, which consists of a live stellar disc, a bulge, and a gaseous component. The dark matter halo is represented by a static potential in order to focus the computational efforts in the disc dynamics. The stellar disc is represented by an exponential-isothermal density profile   $$\rho _{\rm d}(R,z) = \frac{M_{\rm d}}{4 \pi R_{\rm d}^2 z_{\rm d}} \exp \left(-\frac{R}{R_{\rm d}} \right) \mathrm{sech}^2 \left( \frac{z}{z_{\rm d}} \right)\,,$$ (18)where Rd and zd are the radial and vertical scale lengths, respectively. The stellar bulge follows a Hernquist (1990) profile given by   $$\rho _{\rm b}(R) = \frac{M_{\rm b} }{2 \pi R_{\rm b}^3} \frac{1}{(R/R_{\rm b}) (1 + R/R_{\rm b})^3}\,,$$ (19)where Mb is the mass of the bulge and Rb is the scale radius. For the dark matter halo, a Navarro et al. (1997) profile is used, which has a density profile of the form   $$\rho _{\rm h}(R) = \frac{\rho _0}{(R/R_{\rm h}) (1 + R/R_{\rm h})^2}\,,$$ (20)where Rh is the scale radius and ρ0 is a central density parameter. The physical parameters chosen for this model are shown in Table 2, which are based on results from observed data and works modelling the rotation curve of M33 (e.g. Regan & Vogel 1994; Corbelli & Salucci 2000; Hague & Wilkinson 2015; Sie Kam et al. 2015) Table 2. Galaxy model parameters representative of M33. The parameters are: Md – the total baryonic mass of the galactic disc; Rd, zd – the scale length and scale height of the stellar disc; fg – the gas fraction of the disc; f⋆ – the stellar fraction of the disc; Q – the Toomre parameter; m – expected number of spiral modes (at a given radius); ηb – the Efstathiou parameter for bar instability, which in this case is less than the critical value; T – the orbital period of the disc at a given radius; Mb – the mass of the bulge; Rb - the scale length of the bulge, Mh - the mass of the dark matter halo; Rh – the scale length of the halo; c – the halo concentration parameter. Disc    Md(109 M⊙)  9.0  Rd(kpc)  2.5  zd(kpc)  0.2  fg  0.15  f⋆  0.85  Q(1.5Rd)  1.5  m(R = 7kpc)  5.6  ηb  0.9  T(R = 2Rd)(Myr)  284  Bulge    Mb(108 M⊙)  3.0  Rb(kpc)  0.4  Halo    Mh(1011 M⊙)  5.7  Rh(kpc)  33.8  c  4.0  Disc    Md(109 M⊙)  9.0  Rd(kpc)  2.5  zd(kpc)  0.2  fg  0.15  f⋆  0.85  Q(1.5Rd)  1.5  m(R = 7kpc)  5.6  ηb  0.9  T(R = 2Rd)(Myr)  284  Bulge    Mb(108 M⊙)  3.0  Rb(kpc)  0.4  Halo    Mh(1011 M⊙)  5.7  Rh(kpc)  33.8  c  4.0  View Large The initial conditions are obtained by using the method of McMillan & Dehnen (2007), which is publicly available through the code mkgalaxy.3 It generates a self-consistent model of a galaxy composed by a halo, disc, and bulge. It has the advantage of creating a stellar disc with initial velocities sampled from a distribution function representative of a disc rather than assuming a local Maxwellian distribution in velocity space (Dehnen 1999). The disc particles have a velocity distribution consistent with the rotation curve of the galaxy. However, this code only produces a collisionless model of the galaxy and the gas initialization has to be treated separately. To build the full galaxy model, a disc component containing the total number of particles and the total combined stellar and gas mass is first generated. Then, to initialize the gaseous component, the original disc is divided into the two groups of particles and are assigned the following particle masses: mg = Mg/Ng, m⋆ = M⋆/N⋆ for the gas and stellar components, respectively. At this point, the model has a gas component and a stellar component, with the positions and velocities obtained from the initial conditions generator. This is not an equilibrium condition for the gas itself, so the system has to be evolved for some time in order to allow it to settle in the galaxy. This scheme has been previously tested in Ramón-Fox & Aceves (2014) and a similar method also in Dobbs et al. (2010). Fig. 7 shows the full simulation, its arm and inter-arm components, and fits to the spiral structure. The flocculent spiral structures are clearly identified by classification using the velocity shear tensor (E = 2), but the spine fitting algorithm only outlines the larger arms with any clarity. Figure 7. View largeDownload slide The evolution of a galactic disc under a live stellar potential. Top left: the full SPH simulation containing 2 m gas particles, and 2 m star particles. Top right: the gas particles identified as ‘spiral-like’ (filament class). Bottom left: gas particles identified as moving in a disc configuration (sheet class). Bottom right: fits to the spiral structure. Arrows indicate four spiral arms selected for further analysis (see Fig. 8). Figure 7. View largeDownload slide The evolution of a galactic disc under a live stellar potential. Top left: the full SPH simulation containing 2 m gas particles, and 2 m star particles. Top right: the gas particles identified as ‘spiral-like’ (filament class). Bottom left: gas particles identified as moving in a disc configuration (sheet class). Bottom right: fits to the spiral structure. Arrows indicate four spiral arms selected for further analysis (see Fig. 8). In Table 3, we presented calculated pitch angles and pattern speeds for the longest wave in the upper left and lower right quadrants of the galaxy, as well as two other shorter waves in the lower left quadrant (identified by arrows in the bottom right-hand panel of Fig. 7). The long waves yield corotation radii of approximately 1.8 and 4.2 kpc, respectively (Fig. 8). These values bracket the scale length of the stellar disc ($$R_{\rm d} = 2.5 \, \rm {kpc}$$), and denote the region where the rotation curves of both the stellar and gaseous components begin to flatten, and the stellar bulge and gaseous disc enter corotation. Figure 8. View largeDownload slide The rotation curve of the galactic disc under a live stellar potential, with the pattern speeds of four spiral waves marked by horizontal lines. Regions with full lines indicate the radial extent of the spiral. Dashed lines are added to illustrate the corotation radius, denoted by where the horizontal lines meet the rotation curve. Shaded regions indicate the 1σ uncertainties in the pattern speed. The colours correspond with the arrow indicators in the bottom right-hand panel of Fig. 7. Figure 8. View largeDownload slide The rotation curve of the galactic disc under a live stellar potential, with the pattern speeds of four spiral waves marked by horizontal lines. Regions with full lines indicate the radial extent of the spiral. Dashed lines are added to illustrate the corotation radius, denoted by where the horizontal lines meet the rotation curve. Shaded regions indicate the 1σ uncertainties in the pattern speed. The colours correspond with the arrow indicators in the bottom right-hand panel of Fig. 7. Table 3. Fits to the four spiral arms selected from the galaxy simulation with a live stellar potential. The arms 1, 2, 3, and 4 correspond to the blue, green, orange, and purple arrows in Fig. 7 (and similarly for the lines in Fig. 8). We compute uncertainties for Rcorot by computing its value at the upper and lower 1σ limits for Ωp. Arm  ϕ (°)  Ωp (km s−1 kpc−1)  Rcorot (kpc)  1  17.59 ± 0.87  42.57 ± 10.59  1.82$$^{+1.01}_{-0.61}$$  2  16.81 ± 0.85  24.34 ± 7.54  4.25$$^{+2.62}_{-1.41}$$  3  13.01 ± 0.46  7.27 ± 0.28  17.17$$^{+0.81}_{-0.61}$$  4  16.88 ± 1.52  14.86 ± 6.14  7.88$$^{+6.05}_{-2.62}$$  Arm  ϕ (°)  Ωp (km s−1 kpc−1)  Rcorot (kpc)  1  17.59 ± 0.87  42.57 ± 10.59  1.82$$^{+1.01}_{-0.61}$$  2  16.81 ± 0.85  24.34 ± 7.54  4.25$$^{+2.62}_{-1.41}$$  3  13.01 ± 0.46  7.27 ± 0.28  17.17$$^{+0.81}_{-0.61}$$  4  16.88 ± 1.52  14.86 ± 6.14  7.88$$^{+6.05}_{-2.62}$$  View Large The smaller spiral waves reside inside this radius, and appear to have much longer corotation radii (approximately 7 and 17 kpc, respectively). These appear to be less robust, transient features. 4 DISCUSSION By identifying simulation regions that are ‘spiral-like’, we can use the complement of these regions to identify inter-arm regions. This is a valuable distinction for studies of star formation in spiral galaxies. We can use this technique to trace the passage of a spiral shock through gas, recording the density, temperature, and chemical evolution of the gas as it forms molecular clouds, and eventually stars. For example, as gas transitions from the warm diffuse phase into the dense cool phase, we will be able to time this transition relative to when the gas entered and left a spiral structure. This will be of great use in investigating the deviation of low-surface density clouds from the Kennicutt–Schmidt star formation relation (Bonnell et al. 2013), and the relationship in general between spiral structure and molecular clouds (Duarte-Cabral & Dobbs 2016). Spiral structure also drives amplification and reversal of the galactic magnetic field, particularly at large velocity jumps over the spiral shock and at corotation (Dobbs et al. 2016). This process (of relevance to the boundness of molecular clouds) can be studied in finer detail if the fluid elements driving these magnetic amplifications and reversals can be more rigorously identified. Being able to separate the inter-arm component from the arm component has its own benefits. When computing ‘unperturbed’ disc properties, it is common to simply take azimuthal averages of the disc, and assuming the perturbed component is insignificant compared to the wider disc. If the spiral perturbation amplitude is large, this assumption fails. Being able to compute averages in the knowledge that the spiral structure has been extracted will allow a more accurate computation of the spiral perturbation itself (ΔΣ/Σ), as well as a decontaminated study of the unperturbed disc. Foyle et al. (2011) note that star formation rate tracers can be prominent in the inter-arm regions of local spiral galaxies – effective characterization of inter-arm regions in simulations is needed to investigate this phenomenon, and link it to the triggering of star formation by spiral arms, as well as in situ star formation processes. Our spiral identification technique can also be applied to the stellar component (provided that the velocity shear tensor can be appropriately calculated for individual star particles). This will allow us to measure offsets between stellar and gaseous spiral structure (cf. Pettitt et al. 2016), as well as track the motion of stars relative to both the gaseous and stellar arms. This has extra relevance to attempts to elucidate the Milky Way's spiral structure. For example, using molecular clouds to trace spiral arms has demonstrable systematic offsets due to use of the kinematic distance (Ramón-Fox & Bonnell 2018). The feeding of AGNs begins with the outward transport of angular momentum at large scales, allowing matter to flow into the inner 2–3 kpc (where stellar bars can take over). Lagrangian methods like SPH will be able to trace a fluid element's journey from large distances towards the nucleus, recording its encounters with spiral structure along the way. Spiral classification has important benefits for the study of self-gravitating protostellar discs, and discs which are perturbed by internal and external companions. In the first instance, full characterization of a spiral arm can identify whether it is being self-generated by the disc or is due to a companion (Meru et al. 2017). Strong spiral arms in protostellar discs are typically induced by self-gravity, but non-axisymmetric structure can also be induced by the magneto-rotational instability (e.g. Nelson 2005). Spiral classification will provide important discriminants between the two. Spiral arms in protostellar discs drive chemical evolution (Evans et al. 2015) and potentially grain growth (Rice et al. 2004; Booth & Clarke 2016). Timing the perturbation of gas properties (and dust particles) to its passage through spiral arm passage will be critical to understanding their effectiveness. Spiral classification is also valuable for the study of fragmentation in self-gravitating discs. The interaction of fragments with spiral arms has a significant effect on their survival rate (Hall, Forgan & Rice 2017) and their chemical inventory (Ilee et al. 2017). Being able to observe how material is transferred between the spiral arm and the fragment will provide important insights into how fragments accrete, and how their orbital elements, physical properties, and chemistry are governed by interactions with spiral structure. Spirals driven by companions in protostellar discs offer key diagnostics of the disc and the companion, as mentioned in the ‘Introduction’ section. Automatic characterization of spirals allows simulators to run increasingly large banks of simulations over a wider parameter space. This opens up the possibility of potentially identifying new diagnostics from spiral arm data, as well as improved fitting of observations from grids of model runs. Throughout this paper, we have focused on classification via a single tensor, and the choice of tensor has depended largely on the physics at play. In the majority of cases, where the spiral perturbation amplitude is relatively weak, spiral detection proceeds best via the velocity shear tensor (and identifying E = 2 particles). In the limit where surface density perturbations are large, tidal tensor analysis is usually preferred. Our example of the self-gravitating protostellar disc is relatively extreme, in that we find best results for E = 3 particles. Our choices have reflected some knowledge of the physics of the system – the converse is equally true, that analysing using different tensors yields insight into the governing physics. In Forgan et al. (2016a), we showed the benefits of classifying a system using multiple tensors, and then correlating the data. In an example showing the classification of structure around a supernova explosion, the velocity shear tensor identifies regions at the inner edge of the cavity driven by the supernova; the tidal tensor identifies regions being swept into collapsing structures; correlating data identifies regions under collapse as a direct result of the supernova explosion. Spiral detection under multiple tensors can yield important information about the relative amplitudes of individual arms. Finally, it is worth noting that our algorithm is made of two largely separate components: tensor classification and spiral spine identification. The latter component can be swapped for another algorithm depending on the user's aims. Our spine identification approach works best on the most dense regions of the simulation, and is therefore not well suited to the outer regions (as can be seen by its failure to find some structures visible ‘by eye’). Ideally, a spine identification algorithm should use the tensor adopted for classification. For example, arms detected using the velocity shear tensor should be detected using percentiles of some measure of velocity shear (say the maximum eigenvalue) rather than percentiles of density. We leave this as a route to follow in future work. 5 CONCLUSIONS We have demonstrated analysis techniques that isolate and fit individual spiral arms in astrophysical disc simulations. For simulations of both protostellar and galactic discs, we are able to fit logarithmic spirals to each arm, and compute pitch angles. If multiple snapshots are available, we are also able to compute the pattern speed of spiral density waves, and identify where the wave corotates with the bulk flow. A key innovation of these techniques is the ability to identify when fluid elements are inside a spiral structure, or in the inter-arm region. Our techniques rely purely on computing derivatives of the flow at a spatial location, and the subsequent application of friends-of-friends algorithms. Therefore, our approach is applicable to data from any finite spatial element hydrodynamic solver, such as grid-based solvers using adaptive mesh refinement, such as ENZO (Bryan et al. 2014), Voronoi mesh solvers such as AREPO (Springel 2010), or indeed meshless hydrodynamic systems such as GIZMO (Hopkins 2015). We expect our methods will prove extremely useful to simulators studying the role of spirals in the evolution of disc structures throughout the Universe. ACKNOWLEDGEMENTS DHF, FGR-F, and IAB gratefully acknowledge support from the ECOGAL project, grant agreement 291227, funded by the European Research Council under ERC-2011-ADG. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. This work relied on the compute resources of the St Andrews MHD Cluster, and the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC facility (www.dirac.ac.uk). Surface density plots were created using SPLASH (Price 2007). Corner plots were produced using the corner.py module (Foreman-Mackey 2016). Our code is published on Github at https://github.com/dh4gan/tache. The authors warmly thank the anonymous reviewers for their careful reading and astute suggestions. Footnotes 1 In practice, E is defined as the number of particles whose eigenvalues exceed a small, non-zero threshold, see Forgan et al. (2016a) for details. 2 Tensor classification is designed to operate on smooth continuous fields. This is not a concern as SPH fields are guaranteed to be smooth, but it is unclear how strong discontinuities will affect tensor analysis. If reliable derivatives of a field cannot be computed for a simulation region, that region cannot be reliably classified. 3 The mkgalaxy code can be obtained from the NEMO Stellar Dynamics Toolbox package (Teuben 1995) (https://bima.astro.umd.edu/nemo/). REFERENCES Alexander D., Hickox R., 2012, New Astron. Rev. , 56, 93 CrossRef Search ADS   Artymowicz P., Lubow S. H., 1994, ApJ , 421, 651 CrossRef Search ADS   Bate M. R., Bonnell I. A., Bromm V., 2003, MNRAS , 339, 577 CrossRef Search ADS   Benisty M. et al.  , 2015, A&A , 578, L6 CrossRef Search ADS   Binney J., Tremaine S., 2008, Galactic Dynamics . Princeton Univ. Press, University of Princeton Bonnell I. A., Bate M. R., 1994, MNRAS , 271, 999 CrossRef Search ADS   Bonnell I. A., Dobbs C. L., Robitaille T. P., Pringle J. E., 2006, MNRAS , 365, 37 CrossRef Search ADS   Bonnell I. A., Dobbs C. L., Smith R. J., 2013, MNRAS , 430, 1790 CrossRef Search ADS   Booth R. A., Clarke C. J., 2016, MNRAS , 458, 2676 CrossRef Search ADS   Bryan G. L. et al.  , 2014, ApJ , 211, 19 CrossRef Search ADS   Choi Y., Dalcanton J. J., Williams B. F., Weisz D. R., Skillman E. D., Fouesneau M., Dolphin A. E., 2015, ApJ , 810, 9 CrossRef Search ADS   Clarke C. J., Lodato G., 2009, MNRAS , 398, L6 CrossRef Search ADS   Corbelli E., Salucci P., 2000, MNRAS , 311, 441 CrossRef Search ADS   Cossins P., Lodato G., Clarke C. J., 2009, MNRAS , 393, 1157 CrossRef Search ADS   Cox D., Gómez G., 2002, ApJS , 142, 261 CrossRef Search ADS   Davis B. L. et al.  , 2015, ApJ , 802, L13 CrossRef Search ADS   Dehnen W., 1999, AJ , 118, 1201 CrossRef Search ADS   Dobbs C., Baba J., 2014, Publ. Astron. Soc. Aust. , 31, e035 CrossRef Search ADS   Dobbs C. L., Bonnell I. A., 2007, MNRAS , 374, 1115 CrossRef Search ADS   Dobbs C. L., Theis C., Pringle J. E., Bate M. R., 2010, MNRAS , 403, 625 CrossRef Search ADS   Dobbs C. L., Price D. J., Pettitt A. R., Bate M. R., Tricco T. S., 2016, MNRAS , 461, 4482 CrossRef Search ADS   Dong R., Hall C., Rice K., Chiang E., 2015, ApJ , 812, L32 CrossRef Search ADS   Duarte-Cabral A., Dobbs C. L., 2016, MNRAS , 458, 3667 CrossRef Search ADS   Egusa F., Sofue Y., Nakanishi H., 2004, PASJ , 56, L45 CrossRef Search ADS   Egusa F., Kohno K., Sofue Y., Nakanishi H., Komugi S., 2009, ApJ , 697, 1870 CrossRef Search ADS   Elmegreen B. G., 1990, Ann. New York Acad. Sci. , 596, 40 CrossRef Search ADS   Evans M. G., Ilee J. D., Boley A. C., Caselli P., Durisen R. H., Hartquist T. W., Rawlings J. M. C., 2015, MNRAS , 453, 1147 CrossRef Search ADS   Foreman-Mackey D., 2016, J. Open Source Softw. , 24 Forero-Romero J. E., Hoffman Y., Gottlöber S., Klypin a., Yepes G., 2009, MNRAS , 396, 1815 CrossRef Search ADS   Forgan D., Rice K., 2011, MNRAS , 417, 1928 CrossRef Search ADS   Forgan D., Rice K., 2013, MNRAS , 432, 3168 CrossRef Search ADS   Forgan D. H., Rice K., Stamatellos D., Whitworth A. P., 2009, MNRAS , 394, 882 CrossRef Search ADS   Forgan D., Rice K., Cossins P., Lodato G., 2011, MNRAS , 410, 994 CrossRef Search ADS   Forgan D., Parker R. J., Rice K., 2015, MNRAS , 447, 836 CrossRef Search ADS   Forgan D., Bonnell I., Lucas W., Rice K., 2016a, MNRAS , 457, 2501 CrossRef Search ADS   Forgan D. H., Ilee J. D., Cyganowski C. J., Brogan C. L., Hunter T. R., 2016b, MNRAS , 463, 957 CrossRef Search ADS   Foyle K., Rix H.-W., Walter F., Leroy A. K., 2010, ApJ , 725, 534 CrossRef Search ADS   Foyle K., Rix H.-W., Dobbs C. L., Leroy A. K., Walter F., 2011, ApJ , 735, 101 CrossRef Search ADS   Grand R. J. J., Kawata D., Cropper M., 2012, MNRAS , 426, 167 CrossRef Search ADS   Hague P., Wilkinson M., 2015, ApJ , 800, 15 CrossRef Search ADS   Hahn O., Porciani C., Carollo C. M., Dekel A., 2007, MNRAS , 375, 489 CrossRef Search ADS   Hall C., Forgan D., Rice K., 2017, MNRAS , 470, 2517 CrossRef Search ADS   Hernquist L., 1990, ApJ , 356, 359 CrossRef Search ADS   Heyer M., Dame T., 2015, ARA&A , 53, 583 CrossRef Search ADS   Holmberg E., 1941, ApJ , 94, 385 CrossRef Search ADS   Hopkins P. F., 2015, MNRAS , 450, 53 CrossRef Search ADS   Hubble E. P., 1926, ApJ , 64, 321 CrossRef Search ADS   Ilee J. D., Boley A. C., Caselli P., Durisen R. H., Hartquist T. W., Rawlings J. M. C., 2011, MNRAS , 417, 2950 CrossRef Search ADS   Ilee J. D. et al.  , 2017, MNRAS , 472, 189 CrossRef Search ADS   Kennicutt R. C. J., 1981, AJ , 86, 1847 CrossRef Search ADS   Laughlin G., Rozyczka M., 1996, ApJ , 456, 279 CrossRef Search ADS   Lin C. C., Shu F. H., 1964, ApJ , 140, 646 CrossRef Search ADS   Lintott C. et al.  , 2011, MNRAS , 410, 166 CrossRef Search ADS   Lodato G., Price D. J., 2010, MNRAS , 405, 1212 Lodato G., Rice W. K. M., 2004, MNRAS , 351, 630 CrossRef Search ADS   Lodato G., Rice W. K. M., 2005, MNRAS , 358, 1489 CrossRef Search ADS   Mata-Chavez M. D., Gomez G. C., Puerari I., 2014, MNRAS , 444, 3756 CrossRef Search ADS   McMillan P. J., Dehnen W., 2007, MNRAS , 378, 541 CrossRef Search ADS   Meru F., Juhász A., Ilee J. D., Clarke C. J., Rosotti G. P., Booth R. A., 2017, ApJ , 839, L24 CrossRef Search ADS   Monaghan J. J., 1992, ARA&A , 30, 543 CrossRef Search ADS   Monaghan J. J., 2005, Rep. Prog. Phys. , 68, 1703 CrossRef Search ADS   Murray J. R., 1996, MNRAS , 279, 402 CrossRef Search ADS   Muto T. et al.  , 2012, ApJ , 748, L22 CrossRef Search ADS   Navarro J., Frenk C., White S., 1997, ApJ , 490, 493 CrossRef Search ADS   Nelson R. P., 2005, A&A , 443, 1067 CrossRef Search ADS   Papaloizou J. C., Savonije G. J., 1991, MNRAS , 248, 353 CrossRef Search ADS   Pérez L. M. et al.  , 2016, Science , 353, 1519 CrossRef Search ADS PubMed  Pettitt A. R., Tasker E. J., Wadsley J. W., 2016, MNRAS , 458, 3990 CrossRef Search ADS   Pohl A., Pinilla P., Benisty M., Ataiee S., Juhász A., Dullemond C. P., Van Boekel R., Henning T., 2015, MNRAS , 453, 1768 CrossRef Search ADS   Price D. J., 2007, Publ. Astron. Soc. Aust. , 24, 159 CrossRef Search ADS   Price D. J., 2012, J. Comput. Phys. , 231, 759 CrossRef Search ADS   Querejeta M. et al.  , 2016, A&A , 588, A33 CrossRef Search ADS   Rafikov R. R., 2002, ApJ , 569, 997 CrossRef Search ADS   Ragan S. E., Moore T. J. T., Eden D. J., Hoare M. G., Elia D., Molinari S., 2016, MNRAS , 462, 3123 CrossRef Search ADS   Ramón-Fox F., Aceves H., 2014, in Seigar M., Treuthardt P., eds, ASP Conf. Ser., Vol. 480, Structure and Dynamics of Disk Galaxies . Astron. Soc. Pac., San Francisco, p. 229 Ramón-Fox F. G., Bonnell I. A., 2018, MNRAS , 474, 2028 CrossRef Search ADS   Regan M., Vogel S., 1994, ApJ , 434, 536 CrossRef Search ADS   Rice W. K. M., Armitage P. J., 2009, MNRAS , 396, 2228 CrossRef Search ADS   Rice W. K. M., Lodato G., Pringle J. E., Armitage P. J., Bonnell I. A., 2004, MNRAS , 355, 543 CrossRef Search ADS   Rice W. K. M., Lodato G., Armitage P. J., 2005, MNRAS , 364, L56 CrossRef Search ADS   Rix H.-W., Zaritsky D., 1995, ApJ , 447, 82 CrossRef Search ADS   Rosse E. O., 1850, Phil. Trans. R. Soc. London , 140, 499 CrossRef Search ADS   Schinnerer E. et al.  , 2013, ApJ , 779, 42 CrossRef Search ADS   Schinnerer E. et al.  , 2017, ApJ , 836, 62 CrossRef Search ADS   Seiden P. E., Gerola H., 1979, ApJ , 233, 56 CrossRef Search ADS   Sellwood J. A., Carlberg R. G., 1984, ApJ , 282, 61 CrossRef Search ADS   Sie Kam Z., Carignan C., Chemin L., Amram P., Epinat B., 2015, MNRAS , 449, 4048 CrossRef Search ADS   Springel V., 2010, MNRAS , 401, 791 CrossRef Search ADS   Teuben P., 1995, in Shaw R., Panye H., Hayes J., eds, ASP Conf. Ser., Vol. 77, Astronomical Data Analysis Software and Systems IV . Astron. Soc. Pac., San Francisco, p. 398 Tobin J. J. et al.  , 2016, Nature , 538, 483 CrossRef Search ADS PubMed  Vigan A. et al.  , 2017, A&A , 603, A3 CrossRef Search ADS   Zel'dovich Y. B., 1970, A&A , 5, 84 Zhu Z., Dong R., Stone J. M., Rafikov R. R., 2015, ApJ , 813, 88 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# Classifying and modelling spiral structures in hydrodynamic simulations of astrophysical discs

, Volume 476 (2) – May 1, 2018
12 pages

/lp/ou_press/classifying-and-modelling-spiral-structures-in-hydrodynamic-Ixyvgyc3xF
Publisher
The Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty331
Publisher site
See Article on Publisher Site

### Abstract

Abstract We demonstrate numerical techniques for automatic identification of individual spiral arms in hydrodynamic simulations of astrophysical discs. Building on our earlier work, which used tensor classification to identify regions that were ‘spiral-like’, we can now obtain fits to spirals for individual arm elements. We show this process can even detect spirals in relatively flocculent spiral patterns, but the resulting fits to logarithmic ‘grand-design’ spirals are less robust. Our methods not only permit the estimation of pitch angles, but also direct measurements of the spiral arm width and pattern speed. In principle, our techniques will allow the tracking of material as it passes through an arm. Our demonstration uses smoothed particle hydrodynamics simulations, but we stress that the method is suitable for any finite-element hydrodynamics system. We anticipate our techniques will be essential to studies of star formation in disc galaxies, and attempts to find the origin of recently observed spiral structure in protostellar discs. methods: numerical, stars: formation, ISM: kinematics and dynamics 1 INTRODUCTION Astrophysical disc structures are found across a wide range of scales – from disc galaxies to discs surrounding active galactic nuclei (AGNs), to discs around protostars and even around protoplanets. Almost as commonly, these discs can be perturbed into producing spiral structures. The means by which spiral structures are produced depends on the object being studied. Broadly, spiral structure is caused by some form of unstable non-axisymmetric perturbation. How this perturbation subsequently evolves into spiral structure depends on its origin, but in most cases is caused by the disc being susceptible to gravitational instability, where (for the gas) the Toomre parameter is   $$Q = \frac{c_{\rm s} \kappa }{\pi {\rm G} \Sigma } \lessapprox 1,$$ (1)where cs is the sound speed of the gas, κ is the epicyclic frequency, G is the gravitational constant, and Σ is the local disc surface density. A similar expression exists for the stability of a stellar disc (where the radial stellar velocity dispersion replaces cs). Spiral structures have been observed in galaxies for some 150 yr, going back to Rosse (1850)’s observations of M51. We now know that spiral galaxies constitute over half of observed massive galaxies (Lintott et al. 2011). This has provided theorists with significant data reserves on which to test spiral generation theories. Spiral structures observed in disc galaxies can be self-generated, either by local instabilities/noise which are swing amplified to generate arms (e.g. Sellwood & Carlberg 1984), or through globally propagating (quasi-stationary) density waves, as predicted by the pre-eminent Lin–Shu density wave theory (Lin & Shu 1964). Spiral structures can also be externally generated through tidal interactions, e.g. during galaxy mergers [Holmberg (1941)’s experiments in this area are particularly illuminating]. Indeed, all three mechanisms (instabilities, density waves, and interactions) can and do collaborate to produce the spiral structures observed both in observations and simulations (see Dobbs & Baba 2014, for a detailed review). Spiral arms drive important physical processes in galaxies. Molecular clouds and recent star formation tend to be concentrated in and near spiral arms (Schinnerer et al. 2013; Heyer & Dame 2015; Ragan et al. 2016; Schinnerer et al. 2017), most likely due to the compression and shocking of gas as it falls into the arm potential (Bonnell et al. 2006; Dobbs & Bonnell 2007; Bonnell, Dobbs & Smith 2013). Spiral morphology is also directly correlated with the surface density of neutral atomic hydrogen in the disc, and the central stellar bulge mass (Davis et al. 2015). The gravitational torques induced by spiral structures will drive material into the inner few kpc of galaxies. The potential asymmetries caused by stellar bars can then deliver this material further inwards, feeding AGNs (see e.g. Alexander & Hickox 2012; Querejeta et al. 2016). The wide-ranging effects of spirals on galaxy evolution have driven a great deal of effort on determining their properties in observations. Attempts to characterize observed spiral structure in galaxies began with visual classification of the total number of arms, and their ‘openness’ (Hubble 1926), which can be quantified by the pitch angle ϕ. The modern Hubble sequence divides spiral galaxies into barred (SB) and unbarred (S) spirals, with a further subclass (a–d) denoting the tightness of the spiral winding. Elmegreen (1990) proposed dividing spirals based on their arm number (grand design, multi-armed, and flocculent). More quantitative measurements of observed spiral galaxies involve Fourier analysis of the image (Rix & Zaritsky 1995; Foyle et al. 2010), which decomposes the azimuthal variations in surface brightness into a Fourier sum over m spiral modes with amplitude Am, e.g.:   $$\frac{\mu (R, \theta )}{\bar{\mu }(R)} = \sum _{m=1}^{\infty } A_m (R) {\rm e}^{{\rm i}m(\theta -\theta _m)}.$$ (2)Pitch angles can be fitted to images assuming a logarithmic spiral of constant ϕ, a prediction of density wave theory (Kennicutt 1981), or via hyperbolic spirals with radially varying pitch angle (e.g. Seiden & Gerola 1979). The pattern speed of an arm can be investigated by studying star formation (and tracers of star formation) inside the arm (Egusa, Sofue & Nakanishi 2004; Egusa et al. 2009). A steady arm of constant pattern speed will set up a sequence of tracers – e.g. HI, CO, 24 $$\mu$$m emission from enshrouded stars, and UV emission from unobscured stars – from upstream to downstream of the arm's corotation. The failure to see this sequence in local spiral galaxies is strong evidence that spiral structures do not persist beyond the local dynamical time (Foyle et al. 2011). Other observational techniques can also provide convincing evidence for the origin of spiral structures. For example, Choi et al. (2015) resolve the stellar populations along the north-east arm of M81, to show it does not possess a constant pattern speed. This points to a kinematic origin, driven by tidal interactions with nearby galaxies M82 and NGC 3077s. Protostellar discs can also self-generate spiral structure through local instabilities. Shortly after the formation of a protostellar system, the protostellar disc has a mass comparable to that of the central star. Such discs can soon assume a marginally unstable state where Q ∼ 1, and the disc mass determines the spiral modes present (Lodato & Rice 2005; Forgan et al. 2011). Interactions with a companion either internal to the disc (such as a protoplanet) or external to it (such as a close stellar encounter) also generate tidally induced spiral arms. Spiral arms in protostellar discs are also efficient outward transporters of angular momentum, driving rapid accretion and assembly of protostars (Laughlin & Rozyczka 1996). Spiral structures with a sufficient surface density contrast can concentrate dust grains (Rice et al. 2004; Clarke & Lodato 2009; Booth & Clarke 2016), promoting grain growth and setting the initial conditions for planet formation via core accretion, as well as altering local chemistry (Ilee et al. 2011; Evans et al. 2015; Ilee et al. 2017). In extremis, spiral arms of sufficiently large density amplitude can induce protostellar discs to fragment into bound objects (Rice, Lodato & Armitage 2005; Forgan & Rice 2011; Tobin et al. 2016), providing an alternate formation channel for gas giants and substellar objects at large orbital semimajor axis (Forgan & Rice 2013; Forgan, Parker & Rice 2015; Vigan et al. 2017). For sufficiently massive protostellar systems, fragmentation can also generate binary star systems (Bonnell & Bate 1994; Bate, Bonnell & Bromm 2003; Tobin et al. 2016). Spiral structure has only recently been observed in protoplanetary discs, initially in near-infrared observations of scattered light, which is typically most sensitive to the upper disc surface. In particular, extended two-armed structures have been detected around SAO206462 (Muto et al. 2012) and MWC 758 (Benisty et al. 2015). These arms are detectable out to relatively large distances (up to 100 au) from the parent star, with pitch angles of the order of 10°. Recent ALMA observations have shown discs with spiral structure that extends down to the disc midplane, for example the recent detection of m = 2 spiral structure around Elias 2-27, with a measured pitch angle of 7.9±0.4° (Pérez et al. 2016) – the origin of this structure is not yet clear (Meru et al. 2017). Characterizing these spirals, and determining their origin, yields crucial information about protostellar accretion and the protostar's approach to the main sequence, as well as the formation of planetary systems. If we are to identify the origin of newly observed protostellar disc spiral structures, or to study how spiral arms govern and drive star formation in galaxies, analysing spiral structure driven in hydrodynamic simulations is crucial. As such, we require tools to identify and characterize spiral arms in these simulations. Regardless of scale, characterizing spiral morphology yields important diagnostics of what is driving the spiral structure. In particular, the number of arms, their amplitude and pitch angles are sensitive to both the driving mechanism and the disc's properties. For example, the spiral wake produced by planets embedded in a protostellar disc adopts a pitch angle which is a function of the disc temperature and rotation profile, as well as the planetary mass (Rafikov 2002; see also Pohl et al. 2015; Zhu et al. 2015). Gravitationally unstable discs also produce spiral structure with pitch angles and arm number that depend in particular on the disc mass (Dong et al. 2015), which may be a good deal larger than the observed disc mass (Forgan et al. 2016b). What is clear is that both theoretical and observational astrophysicists stand to gain a great deal from higher quality characterization of spiral structure in numerical simulations of astrophysical discs. Theorists gain important insights into how spiral structures are generated, and how they affect the thermal and chemical history of gas and dust, and the future evolution of the disc. Observers gain diagnostics for what is driving the spiral structure in their observations, and glean information on disc properties that is generally orthogonal to other methods. Most attempts to characterize spiral structure in simulations rely on Fourier decomposition (e.g. Cossins, Lodato & Clarke 2009; Dobbs et al. 2010; Forgan et al. 2011; Mata-Chavez, Gomez & Puerari 2014; Pettitt, Tasker & Wadsley 2016). This gives important information on the relative strengths of the spiral modes at play, and the pitch angle and pattern speed of the dominant mode. However, it does not give the pitch angle and pattern speed of individual arms. It also does not inform us as to what sections of the disc are currently in the spiral (and which sections are in the inter-arm regions). Recent attempts at spiral arm characterization have moved away from this Fourier analysis. For example, Grand, Kawata & Cropper (2012) attempted to directly identify individual stellar spiral arms in N-body/hydrodynamic simulations of a barred spiral galaxy, by finding the location of a series of density peaks in a range of annuli. However, this does not directly provide data on the arm width, or which fluid elements currently reside in the spiral. In this paper, we demonstrate that judicious use of tensor classification on hydrodynamic simulation data (Forgan et al. 2016a) allows the identification of fluid elements that are inside spiral structures (or in the inter-arm regions). Further analysis allows the isolation of individual arms to obtain their shape parameters, as well as the pattern speed of the wave. Our examples focus on smoothed particle hydrodynamics (SPH) simulations, but we emphasize that our methods only require the ability to compute derivatives, and are therefore applicable to any hydrodynamic simulation. 2 METHODS Our spiral detection algorithm has two distinct parts. In the first part, fluid elements undergo tensor classification (Forgan et al. 2016a) to determine whether their behaviour indicates they are in fact inside a spiral structure. In the second part, the fluid elements identified as spirals are extracted from the main simulation, and a friends-of-friends algorithm is used on this population to identify the spine of each individual spiral. We describe these procedures below. Our code is published on Github at https://github.com/dh4gan/tache. 2.1 Tensor classification We follow the same procedure as described in Forgan et al. (2016a), which itself builds on work originally applied to N-body simulations of the cosmic Web (see e.g. Hahn et al. 2007; Forero-Romero et al. 2009). In our formalism, tensor classification determines the topology of a chosen field at the location of a given SPH particle. We consider the topology of either the gravitational potential Φ, by computing the tidal tensor Tij as follows:   $$T_{ij} = \frac{\mathrm{\partial} ^2 \Phi }{ \mathrm{\partial} x_i \mathrm{\partial} x_j},$$ (3)or the velocity field via the velocity shear tensor σij as follows:   $$\sigma _{ij} = -\frac{1}{2}\left(\frac{\mathrm{\partial} v_i}{\mathrm{\partial} x_j} + \frac{\mathrm{\partial} v_j}{\mathrm{\partial} x_i} \right).$$ (4)The classification proceeds as follows. Once the tensor to be used (T or σ) is selected, its eigenvalues λi and their corresponding eigenvectors $$\mathbf {n_i}$$ are computed for every particle, e.g.:   $$T \mathbf {n}_j = \lambda _j \mathbf {n}_j.$$ (5)The eigenvalues are defined so λ1 ≥ λ2 ≥ λ3. As T and σ are real and symmetric, these eigenvalues are always real. Both tensors assume that the fields being investigated (potential, velocity) are smooth and continuous, so that the derivatives are always defined. We then compute E, which is defined as the number of positive eigenvalues.1 The value of E determines a topological classification for each fluid element: E = 0 → ‘void’ (0-D manifold) E = 1 → ‘sheet’ (1-D manifold) E = 2 → ‘filament’ (2-D manifold) E = 3 → ‘cluster’ (3-D manifold). We can see this by considering the tidal tensor in the context of Zeldovich theory (Zel'dovich 1970). A test particle orbiting a local extremum ∇Φ = 0 has the following (linearized) equation of motion:   $$\ddot{x}_{i} = -T_{ij} x_j.$$ (6)If we operate in a basis where T is diagonal, then the stability of the orbit is determined by the sign of the eigenvalues λi. A single positive eigenvalue ensures a stable orbit along the corresponding axis (or a sheet-like structure, if multiple particles are present). Two positive eigenvalues allow an extra degree of freedom (filaments) and finally all three being positive allow orbits of any degree (clusters). If no eigenvalues are positive, stable orbits cannot be achieved (voids). Tensor classification using the tidal tensor therefore determines the manifold dimension of the local potential, and as such the preferred structure for matter to collapse into if only gravity is present. Classification using the velocity shear tensor instead diagnoses the manifold being sculpted by the flow at any given instant. As a result, the finite elements of any hydrodynamic simulation can be grouped into four components, with each component composed of fluid elements of a matching E. Discs are inherently sheet-like structures, and hence perturbations from axisymmetry are classified as either filaments or clusters. We will see that depending on the strength of spiral structures, fluid elements classified as either filaments or clusters trace spiral structures, and that fluid elements classified as sheets will trace unperturbed disc material in the inter-arm regions. As a general rule, the tidal tensor is better suited to tracing high-amplitude spiral structures that induce a strong surface density perturbation, but the velocity shear tensor is better at revealing low-amplitude spiral structures that produce weak surface density perturbations.2 Tensor classification also yields eigenvectors corresponding to each eigenvalue. In sheets, the eigenvector $$\mathbf {n}_1$$ corresponding to the single positive eigenvalue λ1 defines the normal or symmetry axis of the sheet. In filaments, the eigenvector $$\mathbf {n}_3$$ corresponding to the single negative eigenvalue λ3 provides the flow direction of the filament. We do not use eigenvectors in this analysis, but note that they may also be of interest in, e.g. determining the relative geometry of spiral structures to their host disc. 2.2 Identifying individual spirals Once tensor classification is complete, we extract all particles that correspond to either filament or cluster classification, as these trace the spiral structure in the disc, and discard the other particles. To determine the spines of each spiral arm on this subset, we run a friends-of-friends algorithm as follows. We specify a minimum spherical radius from the centre of the disc, rmin, from which to begin our study. As we are interested in spines, we will also select the top x% percentile in density of the remaining particles. We then select particle i with the largest density, where radius ri > rmin. This particle forms the first component of the spiral's spine. A sphere of radius L is then drawn from ri. The particles inside the sphere are tested, and particle j is selected if it is the densest particle with rj > ri. Particle j then forms the second component of the spiral's spine, j is then set to i, and the algorithm is repeated until a particle can no longer be found that meets the above conditions. Whenever a particle is tested, it is removed from the list of particles and is not tested again, regardless of whether it forms a component of the spiral spine or otherwise. The act of testing assigns the particle to that particular spiral. As a consequence, any particle within a distance L of any spiral spine point belongs to that spiral. If the conditions can no longer be met, a new spiral is begun at the location of the densest particle not yet tested, and the procedure begins again. Note that throughout we are implicitly assuming that the particle density decreases with increasing radius. Three parameters form the identification algorithm: the linking length L, the minimum radius from which to begin, rmin, and the percentile of the population from which to conduct the analysis, x%. The values needed for these three parameters will depend on the resolution and input physics of the simulation being analysed. The linking length L should clearly be smaller than the spiral arm spacing to avoid the algorithm ‘jumping’ between individual arms, while still being larger than the typical smoothing length to obtain a sufficient number of neighbouring particles. The minimum radius should be selected based on prior knowledge of where spiral structure is well resolved (as well as the locations from which the user desires to measure the said structure). The percentile selection depends on how well defined the structures are post-tensor classification. Selecting a low x% will result in easier classification, but will prevent the classification of weak structures. This algorithm delivers Nspiral sets of coordinate points, with each set denoting the spine of a spiral structure. We elect to fit these points via Markov Chain Monte Carlo (MCMC) to a logarithmic spiral   $$r = r_0 + a e^{b\theta },$$ (7)with fit parameters $$(r_0= \sqrt{x_0^2 + y^2_0}, a,b)$$. The origin of the spiral is defined at θ = 0; r = r0 + a. A pure logarithmic spiral will have a constant pitch angle ϕ, which in radians is simply   $$\phi = \arctan b.$$ (8) It is worth noting that if the spiral arms exclusively constitute the densest regions of our simulation, then it is likely that tensor classification will not be necessary to determine the location of the spiral spines, and that the above friends-of-friends algorithm will work perfectly well. However, if there are other dense structures in the simulation (such as condensing clumps), then the algorithm will attempt to erroneously fit spiral spines to them. We have run test calculations on fragmenting discs to show that tensor classification can successfully remove dense structures that are not spiral arms, preventing this problem. Also, without tensor classification, only the spine of the arm can be determined – the width of the spiral arms or the properties of the gas contained within the arm cannot be measured. 3 TESTS We now test our algorithm on three examples of spiral structure in astrophysical discs. The origins of the structure in each case are subtly different. In the first, we consider a self-gravitating protostellar disc, whose spiral structure is entirely due to the self-gravity of the disc gas. In the second case, we consider a spiral galaxy with an imposed fixed potential that induces well-defined structures. In the last, we consider a galaxy with a ‘live potential’ of star particles, with a resulting spiral structure that is highly flocculent. We use SPH for all three tests. SPH is a Lagrangian hydrodynamics method, where the fluid is discretized into individual particles. Each SPH particle i possesses position and velocity vectors $$\mathbf {r}_i$$, $$\mathbf {v}_i$$, internal energy ui, and smoothing length hi, which is adjusted such that a sphere of radius 2hi encompasses a suitable number N of neighbour SPH particles. The local fluid density is established by performing a sum over the neighbour particles’ smoothing kernels   $$\rho (\mathbf {r}) = \sum ^{N}_{i=1} m_i W(\left|\mathbf {r} - \mathbf {r}_i\right|, h).$$ (9)Writing down the system's Lagrangian and applying variational methods (in combination with equation 9) produce an algorithm for solving the equations of hydrodynamics on the particle distribution (with appropriate compensation for shocks and fluid mixing; see Monaghan 1992, 2005; Price 2012, for reviews). Table 1 shows the values of L, rmin, and x% used in the three tests run in this paper. Table 1. The fit parameters used in the spiral detection algorithm for the test cases evaluated in this paper. Section  Simulation  L  rmin  x%  3.1  Protostellar disc  6 au  10 au  20  3.2  Galaxy, analytic potential  0.2 kpc  1 kpc  30  3.3  Galaxy, live potential  0.2 kpc  1 kpc  30  Section  Simulation  L  rmin  x%  3.1  Protostellar disc  6 au  10 au  20  3.2  Galaxy, analytic potential  0.2 kpc  1 kpc  30  3.3  Galaxy, live potential  0.2 kpc  1 kpc  30  View Large 3.1 A self-gravitating protostellar disc For our first test, we return to a simulation previously analysed in Forgan et al. (2016a). The top left-hand panel of Fig. 1 shows an SPH simulation of a self-gravitating disc of mass 0.25 M⊙, orbiting a sink particle representing a star of 1 M⊙. The initial disc surface density profile follows Σ ∝ r−1, and has a maximum radius of 50 au. The sound speed profile fixes the initial Toomre parameter Q = 2 at all radii. Figure 1. View largeDownload slide Spiral arm detection in a self-gravitating protostellar disc. Top left: a snapshot from the full SPH simulation. Top right: the same snapshot, with only the SPH particles identified as exhibiting arm-like behaviour (through classification of the tidal tensor). Bottom left: the spine points of each arm (green crosses) accompanied by the best-fitting logarithmic spirals in red. Bottom right: The posterior distribution of the logarithmic spiral parameters for the uppermost spiral. Figure 1. View largeDownload slide Spiral arm detection in a self-gravitating protostellar disc. Top left: a snapshot from the full SPH simulation. Top right: the same snapshot, with only the SPH particles identified as exhibiting arm-like behaviour (through classification of the tidal tensor). Bottom left: the spine points of each arm (green crosses) accompanied by the best-fitting logarithmic spirals in red. Bottom right: The posterior distribution of the logarithmic spiral parameters for the uppermost spiral. The radiative transfer formalism of Forgan et al. (2009) is used, implementing realistic local cooling alongside typical compressive and shock heating to settle into a self-regulated, marginally stable state, where the Toomre parameter Q ∼ 1.5. Marginally stable self-gravitating discs are able to sustain quasi-steady spiral structure for as long as the Toomre parameter can be held at this value, which depends critically on the disc's total mass, and its ability to cool. The top right-hand panel of Fig. 1 shows the result of tensor classification using the tidal tensor. Given the deep gravitational perturbations caused by this spiral structure, we elect to use E = 3 (cluster) particles to define the minimum potential of the structure. The bottom row of Fig. 1 shows the result of spiral arm identification. Bottom left shows the identified spine points of each spiral (green points), with accompanying logarithmic fits via MCMC (red curves). An example of the derived posteriors from the MCMC fits is shown in the bottom right-hand panel (for the uppermost spiral). We see that logarithmic spirals can produce good-quality fits to the data, especially in the outer disc. The uppermost spiral is well fitted by a logarithmic spiral with (a, b) = (3.634 ± 0.300, 0.2468 ± 0.0123), giving a pitch angle of 13.8° ± 0.6°. The posterior distributions for a and b are well behaved and unimodal, with only a weak degeneracy between a and b. The inner radii of the spiral shows some deviation from the fit, suggesting that pitch angles may not be constant in this region. This is reflected in the best-fitting centre of the spiral deviating by around 1–2 au from the star's actual position at the origin. Repeating this analysis over subsequent 40 time-steps of the simulation (around 1.5 outer rotation periods, Fig. 2) shows that the matched spirals tend to share relatively similar parameters over this interval. The sample mean pitch angle derived over this duration is 13.02°, with a sample standard deviation of around 3.4°. Interestingly, the scale parameter a tends to assume slightly larger values than our example spiral, with a mean of 4.996 au, and a standard deviation of 1.02 au. Figure 2. View largeDownload slide The distribution of pitch angle ϕ (left) and a over the course of 40 snapshots of the self-gravitating disc simulation (approximately 500 yr of evolution). Figure 2. View largeDownload slide The distribution of pitch angle ϕ (left) and a over the course of 40 snapshots of the self-gravitating disc simulation (approximately 500 yr of evolution). Being able to measure the spiral structure over more than one time-step allows us to calculate the pattern speed of a given arm. Given that for a logarithmic spiral   $$\theta = \frac{1}{b} \ln (r/a),$$ (10)we can estimate the pattern speed $$\Omega _{\rm p} = \dot{\theta }$$ given $$\dot{b}$$ and $$\dot{a}$$:   $$\dot{\theta } = \frac{-\dot{b}}{b} \theta -\frac{\dot{a}}{ab}.$$ (11)Note that $$\dot{\theta }$$ has a dependence on θ only if the pitch angle is time-dependent. Having computed $$\dot{a}$$ and $$\dot{b}$$ from the two fits, we can evaluate the pattern speed over the range θ = [0, 2π] and take an average. For our example spiral, we evaluate (a, b) for two snapshots and compute $$\Omega _{\rm p} = 3.5 \times 10^{-9}\, \rm {rad \, s^{-1}}$$. The corotation radius, where Ω = Ωp, is approximately 14.6 au (Fig. 3). Figure 3. View largeDownload slide The ratio of the pattern speed to the bulk angular velocity Ωp/Ω as a function of radius. Corotation occurs when this ratio is unity, which occurs at 14.6 au, interior to which the disc is dominated by artificial viscosity. Figure 3. View largeDownload slide The ratio of the pattern speed to the bulk angular velocity Ωp/Ω as a function of radius. Corotation occurs when this ratio is unity, which occurs at 14.6 au, interior to which the disc is dominated by artificial viscosity. Artificial viscosity dominates the disc interior to this radius (see e.g. Artymowicz & Lubow 1994; Murray 1996; Lodato & Rice 2004; Lodato & Price 2010; Forgan et al. 2011). The Toomre Q parameter also reaches the instability regime just beyond this radius (∼20 au). As a result, the total disc stress has an artificial minimum at corotation, rather than continuing to decrease with decreasing radius (Rice & Armitage 2009). This shows how the inner disc resolution affects the launching point of spiral structures. Indeed, unstable non-axisymmetric normal modes in self-gravitating gaseous discs are expected to trail the flow at all radii (Papaloizou & Savonije 1991). 3.2 A disc galaxy with fixed arm potential In our second test, we use a galaxy model with analytically defined spiral structure imposed by an external gravitational potential. The galactic potential is represented by a combination of an axisymmetric term plus a perturbation of the spiral arms   $$\Phi = \Phi _{0} + \Phi _{{\rm pert}}.$$ (12) The first component is given by a logarithmic potential (e.g. Binney & Tremaine 2008)   $$\Phi _0 = \frac{1}{2} v^2_0 \left(R^2 + R^2_{\rm c} + (z/z_{\rm q})^2\right),$$ (13)which has a rotation curve given by   $$v_{\rm c}(R) = v_0 \frac{R}{\sqrt{R_{\rm c}^2 + R^2}},$$ (14)where v0 is a velocity parameter, Rc is a characteristic radius, and zq is a vertical scale factor. This produces a flat rotation curve at values larger than Rc. For our simulations, $$v_0 = 120 \, \rm {km \,s^{-1}}$$, Rc = 2.5 kpc, and zq = 0.7. The velocity parameter corresponds to that of the rotation curve of the local spiral galaxy M33 (Corbelli & Salucci 2000; Sie Kam et al. 2015). The spiral arm perturbation uses the scheme of Cox & Gómez (2002)   $$\Phi _{{\rm pert}}(R, \theta ,z, t) = \sum _n A_n(R,z) \cos \left(n \Gamma (R, \theta ,t) \right),$$ (15)where An describes the perturbation amplitude which decays with both R and z (see Cox & Gómez 2002, for details). Γ determines the form of spiral generated by the perturbation   $$\Gamma (R, \theta ,t) = N \left(\theta + \Omega _{\rm p} t - \frac{\ln (R/R_0)}{\tan \phi } - \theta _{\rm p} \right).$$ (16)The number of spiral arms N = 4, and their pattern speed is $$\Omega _{\rm p} = 23 \, \rm {km\,s^{-1}\,kpc^{-1}}$$ (the pitch angle ϕ is 15°). The constant θp defines the location of the arm at a fiducial radius R0. The gaseous component is assumed to have an exponential surface density profile: $$\Sigma (R) = \Sigma _0 e^{-R/R_{\rm d}}$$, where Σ0 is the central surface density and Rd is the scale radius. The gas self-gravity is not included. This analytic potential model is advantageous for its reduced computational expense, and also allows us to have a well-defined spiral structure for testing the spiral detection algorithm. Fig. 4 shows the full SPH simulation (top left), the particles identified as ‘spiral-like’ using the velocity shear tensor (top right), and the particles identified as belonging to the inter-arm regions (E = 1, ‘sheet’ classification, bottom left). The spines of the four main spirals are easily identified (bottom right of Fig. 4), and are well fitted by logarithmic spirals each with b ≈ 0.25 ± 5 per cent (i.e. ϕ ≈ 14° ± 0.7°), representing the input perturbation model well. Figure 4. View largeDownload slide The evolution of a galactic disc under an analytic spiral potential. Top left: the full SPH simulation containing 2 m gas particles. Top right: the particles identified as ‘spiral-like’ (filament class). Bottom left: particles identified as moving in a disc configuration (sheet class). Bottom right: Fits to the spiral structure. Figure 4. View largeDownload slide The evolution of a galactic disc under an analytic spiral potential. Top left: the full SPH simulation containing 2 m gas particles. Top right: the particles identified as ‘spiral-like’ (filament class). Bottom left: particles identified as moving in a disc configuration (sheet class). Bottom right: Fits to the spiral structure. The separation of arm and inter-arm particles in azimuth is clear, if we bin particles by class in an annulus of 0.1 kpc (Fig. 5). We can in fact apply Gaussian fits to the four arms to determine their angular width, resulting in a full width half-maximum (FWHM) of 0.2825 rad. This corresponds to a rather large distance of around 0.5 kpc at this radius, but this is principally due to the winding of the arms aligning them with the azimuthal vector. Figure 5. View largeDownload slide The number of particles classified as ‘arm’ versus ‘inter-arm’ in the galactic disc under the analytic potential, as a function of θ, in an annulus with inner and outer radii R = 1.95, 2.05 kpc. Also plotted are Gaussian fits to the arm histograms, with each Gaussian possessing a full width half-maximum (FWHM) of 0.2825 rad. The galaxy rotates in the positive θ direction. Figure 5. View largeDownload slide The number of particles classified as ‘arm’ versus ‘inter-arm’ in the galactic disc under the analytic potential, as a function of θ, in an annulus with inner and outer radii R = 1.95, 2.05 kpc. Also plotted are Gaussian fits to the arm histograms, with each Gaussian possessing a full width half-maximum (FWHM) of 0.2825 rad. The galaxy rotates in the positive θ direction. We can compute the width of a spiral arm along its entire extent by interrogating the particles that have been identified as belonging to it. We do this for the rightmost spiral in the bottom right-hand panel of Fig. 4, and compute each particle's minimum distance to the spiral spine. The left-hand panel of Fig. 6 shows a 2D histogram of the data (we overplot percentiles of the bin counts for convenience). We also extract 1D histograms by binning in radius (the right-hand panel of Fig. 6) to show that the spiral arm width is essentially constant with radius (with the FWHM of the 1D histograms being approximately 0.1 kpc). Figure 6. View largeDownload slide The width of an individual spiral arm in the galaxy model with analytic potential, as a function of radius. Left: For a population of particles identified by our algorithms as belonging to a specific arm, we compute their minimum distance from the arm and bin them in the above 2D histogram. We overplot percentiles of the bin counts to highlight how the typical separation from an arm is constant with radius. Right: 1D histograms of distance from the arm for the same population, for all radii, and inner and outer radii. Figure 6. View largeDownload slide The width of an individual spiral arm in the galaxy model with analytic potential, as a function of radius. Left: For a population of particles identified by our algorithms as belonging to a specific arm, we compute their minimum distance from the arm and bin them in the above 2D histogram. We overplot percentiles of the bin counts to highlight how the typical separation from an arm is constant with radius. Right: 1D histograms of distance from the arm for the same population, for all radii, and inner and outer radii. As for the protostellar disc, we can compute pattern speeds by comparing two snapshots of the simulation. We find pattern speeds of around $$20.4\, \rm {km\,s^{-1}\,kpc^{-1}}$$ (6.6 × 10−16rad s−1), with a 1σ uncertainty of around 19 per cent. Given the rotation curve, these parameters define the corotation radius as   $$R_{\rm co} = \sqrt{\left(\frac{v_0}{\Omega _{\rm p}}\right)^2 - R^2_{\rm c}} = 5.33 \, \rm {kpc}.$$ (17)The analytic spiral arm potential has a pattern speed of $$23\, \rm {km\,s^{-1}\,kpc^{-1}}$$ (7.45 × 10−16rad s−1), and corotation at 4.9 kpc. Hence, the analytic value resides within the 1σ uncertainties of our calculation. 3.3 A disc galaxy with a live stellar potential We also develop a N-body version of the model galaxy, which consists of a live stellar disc, a bulge, and a gaseous component. The dark matter halo is represented by a static potential in order to focus the computational efforts in the disc dynamics. The stellar disc is represented by an exponential-isothermal density profile   $$\rho _{\rm d}(R,z) = \frac{M_{\rm d}}{4 \pi R_{\rm d}^2 z_{\rm d}} \exp \left(-\frac{R}{R_{\rm d}} \right) \mathrm{sech}^2 \left( \frac{z}{z_{\rm d}} \right)\,,$$ (18)where Rd and zd are the radial and vertical scale lengths, respectively. The stellar bulge follows a Hernquist (1990) profile given by   $$\rho _{\rm b}(R) = \frac{M_{\rm b} }{2 \pi R_{\rm b}^3} \frac{1}{(R/R_{\rm b}) (1 + R/R_{\rm b})^3}\,,$$ (19)where Mb is the mass of the bulge and Rb is the scale radius. For the dark matter halo, a Navarro et al. (1997) profile is used, which has a density profile of the form   $$\rho _{\rm h}(R) = \frac{\rho _0}{(R/R_{\rm h}) (1 + R/R_{\rm h})^2}\,,$$ (20)where Rh is the scale radius and ρ0 is a central density parameter. The physical parameters chosen for this model are shown in Table 2, which are based on results from observed data and works modelling the rotation curve of M33 (e.g. Regan & Vogel 1994; Corbelli & Salucci 2000; Hague & Wilkinson 2015; Sie Kam et al. 2015) Table 2. Galaxy model parameters representative of M33. The parameters are: Md – the total baryonic mass of the galactic disc; Rd, zd – the scale length and scale height of the stellar disc; fg – the gas fraction of the disc; f⋆ – the stellar fraction of the disc; Q – the Toomre parameter; m – expected number of spiral modes (at a given radius); ηb – the Efstathiou parameter for bar instability, which in this case is less than the critical value; T – the orbital period of the disc at a given radius; Mb – the mass of the bulge; Rb - the scale length of the bulge, Mh - the mass of the dark matter halo; Rh – the scale length of the halo; c – the halo concentration parameter. Disc    Md(109 M⊙)  9.0  Rd(kpc)  2.5  zd(kpc)  0.2  fg  0.15  f⋆  0.85  Q(1.5Rd)  1.5  m(R = 7kpc)  5.6  ηb  0.9  T(R = 2Rd)(Myr)  284  Bulge    Mb(108 M⊙)  3.0  Rb(kpc)  0.4  Halo    Mh(1011 M⊙)  5.7  Rh(kpc)  33.8  c  4.0  Disc    Md(109 M⊙)  9.0  Rd(kpc)  2.5  zd(kpc)  0.2  fg  0.15  f⋆  0.85  Q(1.5Rd)  1.5  m(R = 7kpc)  5.6  ηb  0.9  T(R = 2Rd)(Myr)  284  Bulge    Mb(108 M⊙)  3.0  Rb(kpc)  0.4  Halo    Mh(1011 M⊙)  5.7  Rh(kpc)  33.8  c  4.0  View Large The initial conditions are obtained by using the method of McMillan & Dehnen (2007), which is publicly available through the code mkgalaxy.3 It generates a self-consistent model of a galaxy composed by a halo, disc, and bulge. It has the advantage of creating a stellar disc with initial velocities sampled from a distribution function representative of a disc rather than assuming a local Maxwellian distribution in velocity space (Dehnen 1999). The disc particles have a velocity distribution consistent with the rotation curve of the galaxy. However, this code only produces a collisionless model of the galaxy and the gas initialization has to be treated separately. To build the full galaxy model, a disc component containing the total number of particles and the total combined stellar and gas mass is first generated. Then, to initialize the gaseous component, the original disc is divided into the two groups of particles and are assigned the following particle masses: mg = Mg/Ng, m⋆ = M⋆/N⋆ for the gas and stellar components, respectively. At this point, the model has a gas component and a stellar component, with the positions and velocities obtained from the initial conditions generator. This is not an equilibrium condition for the gas itself, so the system has to be evolved for some time in order to allow it to settle in the galaxy. This scheme has been previously tested in Ramón-Fox & Aceves (2014) and a similar method also in Dobbs et al. (2010). Fig. 7 shows the full simulation, its arm and inter-arm components, and fits to the spiral structure. The flocculent spiral structures are clearly identified by classification using the velocity shear tensor (E = 2), but the spine fitting algorithm only outlines the larger arms with any clarity. Figure 7. View largeDownload slide The evolution of a galactic disc under a live stellar potential. Top left: the full SPH simulation containing 2 m gas particles, and 2 m star particles. Top right: the gas particles identified as ‘spiral-like’ (filament class). Bottom left: gas particles identified as moving in a disc configuration (sheet class). Bottom right: fits to the spiral structure. Arrows indicate four spiral arms selected for further analysis (see Fig. 8). Figure 7. View largeDownload slide The evolution of a galactic disc under a live stellar potential. Top left: the full SPH simulation containing 2 m gas particles, and 2 m star particles. Top right: the gas particles identified as ‘spiral-like’ (filament class). Bottom left: gas particles identified as moving in a disc configuration (sheet class). Bottom right: fits to the spiral structure. Arrows indicate four spiral arms selected for further analysis (see Fig. 8). In Table 3, we presented calculated pitch angles and pattern speeds for the longest wave in the upper left and lower right quadrants of the galaxy, as well as two other shorter waves in the lower left quadrant (identified by arrows in the bottom right-hand panel of Fig. 7). The long waves yield corotation radii of approximately 1.8 and 4.2 kpc, respectively (Fig. 8). These values bracket the scale length of the stellar disc ($$R_{\rm d} = 2.5 \, \rm {kpc}$$), and denote the region where the rotation curves of both the stellar and gaseous components begin to flatten, and the stellar bulge and gaseous disc enter corotation. Figure 8. View largeDownload slide The rotation curve of the galactic disc under a live stellar potential, with the pattern speeds of four spiral waves marked by horizontal lines. Regions with full lines indicate the radial extent of the spiral. Dashed lines are added to illustrate the corotation radius, denoted by where the horizontal lines meet the rotation curve. Shaded regions indicate the 1σ uncertainties in the pattern speed. The colours correspond with the arrow indicators in the bottom right-hand panel of Fig. 7. Figure 8. View largeDownload slide The rotation curve of the galactic disc under a live stellar potential, with the pattern speeds of four spiral waves marked by horizontal lines. Regions with full lines indicate the radial extent of the spiral. Dashed lines are added to illustrate the corotation radius, denoted by where the horizontal lines meet the rotation curve. Shaded regions indicate the 1σ uncertainties in the pattern speed. The colours correspond with the arrow indicators in the bottom right-hand panel of Fig. 7. Table 3. Fits to the four spiral arms selected from the galaxy simulation with a live stellar potential. The arms 1, 2, 3, and 4 correspond to the blue, green, orange, and purple arrows in Fig. 7 (and similarly for the lines in Fig. 8). We compute uncertainties for Rcorot by computing its value at the upper and lower 1σ limits for Ωp. Arm  ϕ (°)  Ωp (km s−1 kpc−1)  Rcorot (kpc)  1  17.59 ± 0.87  42.57 ± 10.59  1.82$$^{+1.01}_{-0.61}$$  2  16.81 ± 0.85  24.34 ± 7.54  4.25$$^{+2.62}_{-1.41}$$  3  13.01 ± 0.46  7.27 ± 0.28  17.17$$^{+0.81}_{-0.61}$$  4  16.88 ± 1.52  14.86 ± 6.14  7.88$$^{+6.05}_{-2.62}$$  Arm  ϕ (°)  Ωp (km s−1 kpc−1)  Rcorot (kpc)  1  17.59 ± 0.87  42.57 ± 10.59  1.82$$^{+1.01}_{-0.61}$$  2  16.81 ± 0.85  24.34 ± 7.54  4.25$$^{+2.62}_{-1.41}$$  3  13.01 ± 0.46  7.27 ± 0.28  17.17$$^{+0.81}_{-0.61}$$  4  16.88 ± 1.52  14.86 ± 6.14  7.88$$^{+6.05}_{-2.62}$$  View Large The smaller spiral waves reside inside this radius, and appear to have much longer corotation radii (approximately 7 and 17 kpc, respectively). These appear to be less robust, transient features. 4 DISCUSSION By identifying simulation regions that are ‘spiral-like’, we can use the complement of these regions to identify inter-arm regions. This is a valuable distinction for studies of star formation in spiral galaxies. We can use this technique to trace the passage of a spiral shock through gas, recording the density, temperature, and chemical evolution of the gas as it forms molecular clouds, and eventually stars. For example, as gas transitions from the warm diffuse phase into the dense cool phase, we will be able to time this transition relative to when the gas entered and left a spiral structure. This will be of great use in investigating the deviation of low-surface density clouds from the Kennicutt–Schmidt star formation relation (Bonnell et al. 2013), and the relationship in general between spiral structure and molecular clouds (Duarte-Cabral & Dobbs 2016). Spiral structure also drives amplification and reversal of the galactic magnetic field, particularly at large velocity jumps over the spiral shock and at corotation (Dobbs et al. 2016). This process (of relevance to the boundness of molecular clouds) can be studied in finer detail if the fluid elements driving these magnetic amplifications and reversals can be more rigorously identified. Being able to separate the inter-arm component from the arm component has its own benefits. When computing ‘unperturbed’ disc properties, it is common to simply take azimuthal averages of the disc, and assuming the perturbed component is insignificant compared to the wider disc. If the spiral perturbation amplitude is large, this assumption fails. Being able to compute averages in the knowledge that the spiral structure has been extracted will allow a more accurate computation of the spiral perturbation itself (ΔΣ/Σ), as well as a decontaminated study of the unperturbed disc. Foyle et al. (2011) note that star formation rate tracers can be prominent in the inter-arm regions of local spiral galaxies – effective characterization of inter-arm regions in simulations is needed to investigate this phenomenon, and link it to the triggering of star formation by spiral arms, as well as in situ star formation processes. Our spiral identification technique can also be applied to the stellar component (provided that the velocity shear tensor can be appropriately calculated for individual star particles). This will allow us to measure offsets between stellar and gaseous spiral structure (cf. Pettitt et al. 2016), as well as track the motion of stars relative to both the gaseous and stellar arms. This has extra relevance to attempts to elucidate the Milky Way's spiral structure. For example, using molecular clouds to trace spiral arms has demonstrable systematic offsets due to use of the kinematic distance (Ramón-Fox & Bonnell 2018). The feeding of AGNs begins with the outward transport of angular momentum at large scales, allowing matter to flow into the inner 2–3 kpc (where stellar bars can take over). Lagrangian methods like SPH will be able to trace a fluid element's journey from large distances towards the nucleus, recording its encounters with spiral structure along the way. Spiral classification has important benefits for the study of self-gravitating protostellar discs, and discs which are perturbed by internal and external companions. In the first instance, full characterization of a spiral arm can identify whether it is being self-generated by the disc or is due to a companion (Meru et al. 2017). Strong spiral arms in protostellar discs are typically induced by self-gravity, but non-axisymmetric structure can also be induced by the magneto-rotational instability (e.g. Nelson 2005). Spiral classification will provide important discriminants between the two. Spiral arms in protostellar discs drive chemical evolution (Evans et al. 2015) and potentially grain growth (Rice et al. 2004; Booth & Clarke 2016). Timing the perturbation of gas properties (and dust particles) to its passage through spiral arm passage will be critical to understanding their effectiveness. Spiral classification is also valuable for the study of fragmentation in self-gravitating discs. The interaction of fragments with spiral arms has a significant effect on their survival rate (Hall, Forgan & Rice 2017) and their chemical inventory (Ilee et al. 2017). Being able to observe how material is transferred between the spiral arm and the fragment will provide important insights into how fragments accrete, and how their orbital elements, physical properties, and chemistry are governed by interactions with spiral structure. Spirals driven by companions in protostellar discs offer key diagnostics of the disc and the companion, as mentioned in the ‘Introduction’ section. Automatic characterization of spirals allows simulators to run increasingly large banks of simulations over a wider parameter space. This opens up the possibility of potentially identifying new diagnostics from spiral arm data, as well as improved fitting of observations from grids of model runs. Throughout this paper, we have focused on classification via a single tensor, and the choice of tensor has depended largely on the physics at play. In the majority of cases, where the spiral perturbation amplitude is relatively weak, spiral detection proceeds best via the velocity shear tensor (and identifying E = 2 particles). In the limit where surface density perturbations are large, tidal tensor analysis is usually preferred. Our example of the self-gravitating protostellar disc is relatively extreme, in that we find best results for E = 3 particles. Our choices have reflected some knowledge of the physics of the system – the converse is equally true, that analysing using different tensors yields insight into the governing physics. In Forgan et al. (2016a), we showed the benefits of classifying a system using multiple tensors, and then correlating the data. In an example showing the classification of structure around a supernova explosion, the velocity shear tensor identifies regions at the inner edge of the cavity driven by the supernova; the tidal tensor identifies regions being swept into collapsing structures; correlating data identifies regions under collapse as a direct result of the supernova explosion. Spiral detection under multiple tensors can yield important information about the relative amplitudes of individual arms. Finally, it is worth noting that our algorithm is made of two largely separate components: tensor classification and spiral spine identification. The latter component can be swapped for another algorithm depending on the user's aims. Our spine identification approach works best on the most dense regions of the simulation, and is therefore not well suited to the outer regions (as can be seen by its failure to find some structures visible ‘by eye’). Ideally, a spine identification algorithm should use the tensor adopted for classification. For example, arms detected using the velocity shear tensor should be detected using percentiles of some measure of velocity shear (say the maximum eigenvalue) rather than percentiles of density. We leave this as a route to follow in future work. 5 CONCLUSIONS We have demonstrated analysis techniques that isolate and fit individual spiral arms in astrophysical disc simulations. For simulations of both protostellar and galactic discs, we are able to fit logarithmic spirals to each arm, and compute pitch angles. If multiple snapshots are available, we are also able to compute the pattern speed of spiral density waves, and identify where the wave corotates with the bulk flow. A key innovation of these techniques is the ability to identify when fluid elements are inside a spiral structure, or in the inter-arm region. Our techniques rely purely on computing derivatives of the flow at a spatial location, and the subsequent application of friends-of-friends algorithms. Therefore, our approach is applicable to data from any finite spatial element hydrodynamic solver, such as grid-based solvers using adaptive mesh refinement, such as ENZO (Bryan et al. 2014), Voronoi mesh solvers such as AREPO (Springel 2010), or indeed meshless hydrodynamic systems such as GIZMO (Hopkins 2015). We expect our methods will prove extremely useful to simulators studying the role of spirals in the evolution of disc structures throughout the Universe. ACKNOWLEDGEMENTS DHF, FGR-F, and IAB gratefully acknowledge support from the ECOGAL project, grant agreement 291227, funded by the European Research Council under ERC-2011-ADG. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. This work relied on the compute resources of the St Andrews MHD Cluster, and the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC facility (www.dirac.ac.uk). Surface density plots were created using SPLASH (Price 2007). Corner plots were produced using the corner.py module (Foreman-Mackey 2016). Our code is published on Github at https://github.com/dh4gan/tache. The authors warmly thank the anonymous reviewers for their careful reading and astute suggestions. Footnotes 1 In practice, E is defined as the number of particles whose eigenvalues exceed a small, non-zero threshold, see Forgan et al. (2016a) for details. 2 Tensor classification is designed to operate on smooth continuous fields. This is not a concern as SPH fields are guaranteed to be smooth, but it is unclear how strong discontinuities will affect tensor analysis. If reliable derivatives of a field cannot be computed for a simulation region, that region cannot be reliably classified. 3 The mkgalaxy code can be obtained from the NEMO Stellar Dynamics Toolbox package (Teuben 1995) (https://bima.astro.umd.edu/nemo/). REFERENCES Alexander D., Hickox R., 2012, New Astron. Rev. , 56, 93 CrossRef Search ADS   Artymowicz P., Lubow S. H., 1994, ApJ , 421, 651 CrossRef Search ADS   Bate M. R., Bonnell I. A., Bromm V., 2003, MNRAS , 339, 577 CrossRef Search ADS   Benisty M. et al.  , 2015, A&A , 578, L6 CrossRef Search ADS   Binney J., Tremaine S., 2008, Galactic Dynamics . Princeton Univ. Press, University of Princeton Bonnell I. A., Bate M. R., 1994, MNRAS , 271, 999 CrossRef Search ADS   Bonnell I. A., Dobbs C. L., Robitaille T. P., Pringle J. E., 2006, MNRAS , 365, 37 CrossRef Search ADS   Bonnell I. A., Dobbs C. L., Smith R. J., 2013, MNRAS , 430, 1790 CrossRef Search ADS   Booth R. A., Clarke C. J., 2016, MNRAS , 458, 2676 CrossRef Search ADS   Bryan G. L. et al.  , 2014, ApJ , 211, 19 CrossRef Search ADS   Choi Y., Dalcanton J. J., Williams B. F., Weisz D. R., Skillman E. D., Fouesneau M., Dolphin A. E., 2015, ApJ , 810, 9 CrossRef Search ADS   Clarke C. J., Lodato G., 2009, MNRAS , 398, L6 CrossRef Search ADS   Corbelli E., Salucci P., 2000, MNRAS , 311, 441 CrossRef Search ADS   Cossins P., Lodato G., Clarke C. J., 2009, MNRAS , 393, 1157 CrossRef Search ADS   Cox D., Gómez G., 2002, ApJS , 142, 261 CrossRef Search ADS   Davis B. L. et al.  , 2015, ApJ , 802, L13 CrossRef Search ADS   Dehnen W., 1999, AJ , 118, 1201 CrossRef Search ADS   Dobbs C., Baba J., 2014, Publ. Astron. Soc. Aust. , 31, e035 CrossRef Search ADS   Dobbs C. L., Bonnell I. A., 2007, MNRAS , 374, 1115 CrossRef Search ADS   Dobbs C. L., Theis C., Pringle J. E., Bate M. R., 2010, MNRAS , 403, 625 CrossRef Search ADS   Dobbs C. L., Price D. J., Pettitt A. R., Bate M. R., Tricco T. S., 2016, MNRAS , 461, 4482 CrossRef Search ADS   Dong R., Hall C., Rice K., Chiang E., 2015, ApJ , 812, L32 CrossRef Search ADS   Duarte-Cabral A., Dobbs C. L., 2016, MNRAS , 458, 3667 CrossRef Search ADS   Egusa F., Sofue Y., Nakanishi H., 2004, PASJ , 56, L45 CrossRef Search ADS   Egusa F., Kohno K., Sofue Y., Nakanishi H., Komugi S., 2009, ApJ , 697, 1870 CrossRef Search ADS   Elmegreen B. G., 1990, Ann. New York Acad. Sci. , 596, 40 CrossRef Search ADS   Evans M. G., Ilee J. D., Boley A. C., Caselli P., Durisen R. H., Hartquist T. W., Rawlings J. M. C., 2015, MNRAS , 453, 1147 CrossRef Search ADS   Foreman-Mackey D., 2016, J. Open Source Softw. , 24 Forero-Romero J. E., Hoffman Y., Gottlöber S., Klypin a., Yepes G., 2009, MNRAS , 396, 1815 CrossRef Search ADS   Forgan D., Rice K., 2011, MNRAS , 417, 1928 CrossRef Search ADS   Forgan D., Rice K., 2013, MNRAS , 432, 3168 CrossRef Search ADS   Forgan D. H., Rice K., Stamatellos D., Whitworth A. P., 2009, MNRAS , 394, 882 CrossRef Search ADS   Forgan D., Rice K., Cossins P., Lodato G., 2011, MNRAS , 410, 994 CrossRef Search ADS   Forgan D., Parker R. J., Rice K., 2015, MNRAS , 447, 836 CrossRef Search ADS   Forgan D., Bonnell I., Lucas W., Rice K., 2016a, MNRAS , 457, 2501 CrossRef Search ADS   Forgan D. H., Ilee J. D., Cyganowski C. J., Brogan C. L., Hunter T. R., 2016b, MNRAS , 463, 957 CrossRef Search ADS   Foyle K., Rix H.-W., Walter F., Leroy A. K., 2010, ApJ , 725, 534 CrossRef Search ADS   Foyle K., Rix H.-W., Dobbs C. L., Leroy A. K., Walter F., 2011, ApJ , 735, 101 CrossRef Search ADS   Grand R. J. J., Kawata D., Cropper M., 2012, MNRAS , 426, 167 CrossRef Search ADS   Hague P., Wilkinson M., 2015, ApJ , 800, 15 CrossRef Search ADS   Hahn O., Porciani C., Carollo C. M., Dekel A., 2007, MNRAS , 375, 489 CrossRef Search ADS   Hall C., Forgan D., Rice K., 2017, MNRAS , 470, 2517 CrossRef Search ADS   Hernquist L., 1990, ApJ , 356, 359 CrossRef Search ADS   Heyer M., Dame T., 2015, ARA&A , 53, 583 CrossRef Search ADS   Holmberg E., 1941, ApJ , 94, 385 CrossRef Search ADS   Hopkins P. F., 2015, MNRAS , 450, 53 CrossRef Search ADS   Hubble E. P., 1926, ApJ , 64, 321 CrossRef Search ADS   Ilee J. D., Boley A. C., Caselli P., Durisen R. H., Hartquist T. W., Rawlings J. M. C., 2011, MNRAS , 417, 2950 CrossRef Search ADS   Ilee J. D. et al.  , 2017, MNRAS , 472, 189 CrossRef Search ADS   Kennicutt R. C. J., 1981, AJ , 86, 1847 CrossRef Search ADS   Laughlin G., Rozyczka M., 1996, ApJ , 456, 279 CrossRef Search ADS   Lin C. C., Shu F. H., 1964, ApJ , 140, 646 CrossRef Search ADS   Lintott C. et al.  , 2011, MNRAS , 410, 166 CrossRef Search ADS   Lodato G., Price D. J., 2010, MNRAS , 405, 1212 Lodato G., Rice W. K. M., 2004, MNRAS , 351, 630 CrossRef Search ADS   Lodato G., Rice W. K. M., 2005, MNRAS , 358, 1489 CrossRef Search ADS   Mata-Chavez M. D., Gomez G. C., Puerari I., 2014, MNRAS , 444, 3756 CrossRef Search ADS   McMillan P. J., Dehnen W., 2007, MNRAS , 378, 541 CrossRef Search ADS   Meru F., Juhász A., Ilee J. D., Clarke C. J., Rosotti G. P., Booth R. A., 2017, ApJ , 839, L24 CrossRef Search ADS   Monaghan J. J., 1992, ARA&A , 30, 543 CrossRef Search ADS   Monaghan J. J., 2005, Rep. Prog. Phys. , 68, 1703 CrossRef Search ADS   Murray J. R., 1996, MNRAS , 279, 402 CrossRef Search ADS   Muto T. et al.  , 2012, ApJ , 748, L22 CrossRef Search ADS   Navarro J., Frenk C., White S., 1997, ApJ , 490, 493 CrossRef Search ADS   Nelson R. P., 2005, A&A , 443, 1067 CrossRef Search ADS   Papaloizou J. C., Savonije G. J., 1991, MNRAS , 248, 353 CrossRef Search ADS   Pérez L. M. et al.  , 2016, Science , 353, 1519 CrossRef Search ADS PubMed  Pettitt A. R., Tasker E. J., Wadsley J. W., 2016, MNRAS , 458, 3990 CrossRef Search ADS   Pohl A., Pinilla P., Benisty M., Ataiee S., Juhász A., Dullemond C. P., Van Boekel R., Henning T., 2015, MNRAS , 453, 1768 CrossRef Search ADS   Price D. J., 2007, Publ. Astron. Soc. Aust. , 24, 159 CrossRef Search ADS   Price D. J., 2012, J. Comput. Phys. , 231, 759 CrossRef Search ADS   Querejeta M. et al.  , 2016, A&A , 588, A33 CrossRef Search ADS   Rafikov R. R., 2002, ApJ , 569, 997 CrossRef Search ADS   Ragan S. E., Moore T. J. T., Eden D. J., Hoare M. G., Elia D., Molinari S., 2016, MNRAS , 462, 3123 CrossRef Search ADS   Ramón-Fox F., Aceves H., 2014, in Seigar M., Treuthardt P., eds, ASP Conf. Ser., Vol. 480, Structure and Dynamics of Disk Galaxies . Astron. Soc. Pac., San Francisco, p. 229 Ramón-Fox F. G., Bonnell I. A., 2018, MNRAS , 474, 2028 CrossRef Search ADS   Regan M., Vogel S., 1994, ApJ , 434, 536 CrossRef Search ADS   Rice W. K. M., Armitage P. J., 2009, MNRAS , 396, 2228 CrossRef Search ADS   Rice W. K. M., Lodato G., Pringle J. E., Armitage P. J., Bonnell I. A., 2004, MNRAS , 355, 543 CrossRef Search ADS   Rice W. K. M., Lodato G., Armitage P. J., 2005, MNRAS , 364, L56 CrossRef Search ADS   Rix H.-W., Zaritsky D., 1995, ApJ , 447, 82 CrossRef Search ADS   Rosse E. O., 1850, Phil. Trans. R. Soc. London , 140, 499 CrossRef Search ADS   Schinnerer E. et al.  , 2013, ApJ , 779, 42 CrossRef Search ADS   Schinnerer E. et al.  , 2017, ApJ , 836, 62 CrossRef Search ADS   Seiden P. E., Gerola H., 1979, ApJ , 233, 56 CrossRef Search ADS   Sellwood J. A., Carlberg R. G., 1984, ApJ , 282, 61 CrossRef Search ADS   Sie Kam Z., Carignan C., Chemin L., Amram P., Epinat B., 2015, MNRAS , 449, 4048 CrossRef Search ADS   Springel V., 2010, MNRAS , 401, 791 CrossRef Search ADS   Teuben P., 1995, in Shaw R., Panye H., Hayes J., eds, ASP Conf. Ser., Vol. 77, Astronomical Data Analysis Software and Systems IV . Astron. Soc. Pac., San Francisco, p. 398 Tobin J. J. et al.  , 2016, Nature , 538, 483 CrossRef Search ADS PubMed  Vigan A. et al.  , 2017, A&A , 603, A3 CrossRef Search ADS   Zel'dovich Y. B., 1970, A&A , 5, 84 Zhu Z., Dong R., Stone J. M., Rafikov R. R., 2015, ApJ , 813, 88 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### Monthly Plan • Read unlimited articles • Personalized recommendations • No expiration • Print 20 pages per month • 20% off on PDF purchases • Organize your research • Get updates on your journals and topic searches$49/month

14-day Free Trial

Best Deal — 39% off

### Annual Plan

• All the features of the Professional Plan, but for 39% off!
• Billed annually
• No expiration
• For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588$360/year

billed annually

14-day Free Trial