© The Authors 2023

1. Introduction

The rotation frequency Ω_p, also called the pattern speed, of morphological features such as bars or spiral arms is a crucial parameter in most dynamical studies of disk galaxies because it allows describing the local dynamics in a rotating frame where time dependence is reduced. Physical and numerical reality is more complicated, however, because most self-consistent galaxy models, even when they do not include a dissipative gas component, show that several flexible patterns with distinct speeds coexist (Sparke & Sellwood 1987; Wozniak 2015, 2020; Wu et al. 2016, 2018). This implies that the system is not strictly time-invariant at meso- and macro-scales in any rotating frame. Therefore, periodic or secular angular speed variations must be expected, and the amplitude of the time fluctuations is an aspect of galactic dynamics that needs to be better examined. The constancy of spiral arm pattern speeds has been a fundamental assumption in the density wave theory of Lin & Shu (1964) all along, but this theory has been challenged (Sellwood 2011): Spiral arms might be winding up, or be transient (Toomre 1991) in view of the efficient swing amplification of small perturbations in cold disks (Julian & Toomre 1966). Long-lived (∼5 Gyr) spirals are not excluded (Saha & Elmegreen 2016) in particular conditions, however. More quantitative numerical investigations of these questions in fully consistent models are therefore required. In the meanwhile, evaluations of the instantaneous pattern speeds of spatially limited features are required in order to be able to verify whether long-time pattern speed averages, as commonly calculated, make sense for disk dynamics, and to quantitatively determine the errors when assuming a single global pattern speed.

For solid bodies, there is no internal streaming, thus the instantaneous rotation speed and instantaneous pattern speed are identical. More realistically, in stellar and fluid systems, the mean rotation of their matter and the rotation of one or several patterns may all be different, and differ spatially, as was clearly demonstrated in Wu et al. (2018), because bars and spirals may be flexible structures.

The precise definition of “pattern speed” is rarely specified, while a “pattern” may mean different things: it can be just a special morphological feature, or the shape of a mass distribution in space, kinematical space, or phase-space, or in more abstract spaces, such as frequency (Fourier) space. To characterize a shape, the gradient of observable quantities, such as the mean mass distribution, must be available. The motion of this gradient allows us to follow the shape in time, and, if it rotates uniformly around a center, it can describe the pattern rotation speed. Assuming a single pattern rotation speed therefore assumes that all gradients of the observable quantity over the examined region rotate rigidly, at the same angular speed, around a common center.

That a stellar or particle system displays a single pattern speed everywhere is very particular. Typical galaxies or N-body models of galaxies may display an approximate constant pattern speed over only a subset of the whole. When a pattern speed is estimated, it therefore is important to understand the possible bias introduced by selecting a specific particle subset. Unless a continuous differentiable function is involved, in which case the motion of the function gradient at any point yields a speed, a pattern in galaxies or in simulations is a statistical property over many particles. Choosing a specific subset may therefore introduce a selection bias. For example, the pattern speed of a single particle can only be its speed. The pattern speed of a well-grouped set of particles (a cluster) is its average speed. The pattern speed of a ring of particles is not the average circular speed of its particles, but the speed of a density gradient moving as a single wave. Depending on the subset selection, the speed of a compact subset of stars may therefore reflect the local streaming motion of this subset more than the motion of a pattern of larger size.

In numerical work, most often, in N-body simulations, full information about the particle distribution and properties is available, so that the assumptions made for galaxies can in principle be checked in detail. A popular approach to determine an average pattern speed is to Fourier transform, over angle and time, the particle distribution in radial bins (Sellwood & Athanassoula 1986; Masset & Tagger 1997; Rautiainen & Salo 1999). Because the simulations are not time periodic, a finite time window is used that is longer than the pattern period, but shorter than the simulation length. Power peaks of a given Fourier mode in the time domain at a given radius reveal the average mode pattern speed(s) over the chosen time window. The power spectrum indicates whether a single or multiple distinct patterns coexist on average at any given radius. The sampling of shortest and longest time-intervals bounds the examined average frequencies, but also introduce aliases. Strictly, any power spectrum different from a single sharply peaked function (the width of which should be consistent with the particle noise) then indicates time variability.

A direct straight method for finding the almost instantaneous bar pattern speed has been used, for instance, by Pfenniger & Friedli (1991). The method consists of measuring the angle of its inertia tensor major axis separated by a short time interval, where the inertia tensor is calculated with the particles within a cylinder as large as the bar. This method provides an almost instantaneous pattern speed averaged over a time step, but the result may depend on the adopted cylinder size. Its drawback is that it requires two particle snapshots, which is unnecessary, as we show in Sect. 6.

The current most popular method to find the instantaneous Ω_p has been proposed by Tremaine & Weinberg (1984, hereafter TW), with some modifications in Merrifield et al. (2006) and Meidt et al. (2008), where radial variations are considered. The conservation of the observed population projected in the galactic plane is assumed. Integrating spatial gradients over lines allows recovering the bar pattern speed when only the line-of-sight velocity component and the disk inclination angle are known. The problem with the TW method is that the pattern speed is assumed to be uniform over the whole observed disk, which is questionable in view of the known pattern speed differences between the bar and the spiral arms in N-body simulations (Sparke & Sellwood 1987; Rautiainen & Salo 1999). The possible interaction between the bars and spiral arms might also affect the pattern speed measurement (Hilmi et al. 2020). Furthermore, nested double bars are common (Erwin 2004), and are expected to rotate at very distinct speeds (Pfenniger & Norman 1990), and was indeed found in simulations (e.g., Friedli & Martinet 1993; Debattista & Shen 2007; Saha & Maciejewski 2013; Wu et al. 2016, 2018). A more flexible method giving the instantaneous pattern speed over a spatial domain relevant for the addressed problem appears to be required.

Methods based on the kinematic signature of specific resonances, such as the Milky Way bar outer Lindblad resonance observed in a star sample around the Sun (e.g., Monari et al. 2017), are global and implicitly assume that the disk dynamics at the Sun is not strongly perturbed by the local spiral arms, which rotate possibly with a distinct pattern speed, as is frequently observed in N-body simulations. The mere presence of self-gravitating spiral arms around the Sun suggests the opposite, however: Local forces are substantially dominated by the influence of local disk matter.

Offset bars, now seen in more than a dozen galaxies (de Swardt et al. 2015; Kruk et al. 2017), present another kind of problem. Because the offset bar may not necessarily be at the kinematical center of the host galaxy or at the geometrical center of the outer disk, the bar pattern, including its central density peak, may rotate both about itself and as a whole about the galaxy center, performing a kind of epicyclic motion around the center of mass, a superposition of two or more rotational motions.

Here we concentrate on methods providing a single instantaneous pattern speed and possibly local pattern speeds in N-body models because the information contained in a particle snapshot should be sufficient to find them in principle, and, in some cases, their pattern acceleration as well. In case of multiple rotational components, to first order, we search only for the fastest rotational component in the system and its instantaneous rotation center. The slower additional rotational components are then approximated to first order with a constant translation of the rotation center. Proceeding beyond first order may be vain as the pattern speed invariance hypothesis is likely to be invalid over longer time intervals.

Some of these methods may hopefully be useful for interpreting astrometric data of the Milky Way. Precise astrometric surveys, such as the Gaia survey (Gaia Collaboration 2018), start to provide full phase-space coordinate information for large Milky Way star samples. Alternative methods to determine a pattern speed, in particular of the Milky Way bar, for which a number of methods and tracers already exist (e.g., Debattista et al. 2002; Gerhard 2011; Minchev et al. 2007; Sanders et al. 2019; Clarke & Gerhard 2022), may be found to be useful and complementary.

It might be thought that any method should give the same pattern speed value. Using three of methods exposed here, Wu et al. (2018) showed that according to the adopted shape definition, somewhat different pattern speeds are found in the same model because different methods measure different characteristics of the not strictly rigidly and constantly rotating pattern. This is especially relevant for the double-bar galaxy models used in Wu et al. (2018), or for non-rigidly rotating spiral arms. In regions in which the system is strongly time dependent with respect to the local dynamical time, for example, at the interface between two nested bars or between a bar and the surrounding spiral arms around the corotation radius, methods based on the assumption of a constant pattern speed obviously break down, as shown in Wu et al. (2016, 2018), because the single steady pattern assumption is incorrect.

In Sect. 2 we discuss a naive method along with its shortcomings. In Sect. 3 we derive the scalar or vector pattern speed in space and velocity space from any observable function of space and/or velocity space that is assumed invariant by rotation and translation, with a possibly unknown rotation center. In Sect. 4, as a generalization of the TW method, we examine the local and regional TW methods (herafter LTW and RTW), and describe the criterion where in observations, the derived pattern speed is a priori the least affected by measurement errors. Then, in Sect. 5, we use the Jacobi integral, or equivalently, the rotational invariance of the mean potential, to determine the potential pattern speed. In Sect. 6 we consider how second-moment tensors in ordinary space, or velocity space, allow us to determine Ω_p, when the rotation axis direction is known or not. In Sect. 7 we use the Fourier space phase information to analyze how circular wave modes rotate, which is for mode 2 a similar but distinct approach from the moment methods. Finally, we conclude in Sect. 8.

2. Morphological feature method

A direct method (not counted in the five methods) is to use a clear morphological feature, assumed to be corotating with the figure. For example, a barred spiral can display a straight bar and a symmetric two-arm spiral apparently starting at the ends of the bar. The intersection points of a logarithmic spiral fitted along the spiral arms and a straight line along the bar major axis define two symmetric points that can be assumed to corotate with the bar. Dynamically, the most obvious explanation for this feature is that these points represent the stationary Lagrange points L_1, 2 of the bar pattern (Pfenniger 1990). Other examples in the literature attempted to use the Lagrange points L_4, 5 (e.g., Shimizu & Yoshii 2012), the location of which can be inferred from various tracers, such as a star-forming region or a local minimum of surface brightness in case of unstable L_4, 5 for strong bars.

When we know both the position (x, y) and the velocity in an inertial frame (v_x, v_y) of one of these points in the galaxy plane, the angle is found by

$$\begin{array}{c}\hfill \varphi =arctan(y(t),x(t)),\end{array}$$(1)

and the pattern speed is its time derivative,

$$\begin{array}{c}\hfill \mathrm{\Omega }\equiv \frac{\mathrm{d}}{\mathrm{d}t}arctan(y(t),x(t))=\frac{x{v}_y-y{v}_x}{x^2+y^2}=\frac{v_t}{R},\end{array}$$(2)

where v_t is the tangential velocity component of the point (x, y), and R = (x² + y²)^1/2 is its distance to origin.

This method clearly relies on the assumption that the chosen point is indeed stationary in some rotating frame. For N-body bars, Wu et al. (2016) have shown that the instantaneous equilibrium points between a bar and spiral arms may oscillate with a substantial amplitude in any rotating frame, including the bar frame. No strict Lagrange points therefore exist because of the perturbation by the adjacent spiral pattern that rotates at a different speed. The assumptions made about the pattern that we wish to measure are therefore clearly crucial. In some cases, no strictly constant pattern speed even exists.

3. Invariant function method

A pattern can be defined by the mean of a function f(x, v, t) (e.g., the phase-space distribution function of stars, the mass density of a particular stellar population, or the mean gravitational potential). This function is assumed to be invariant either (1) after an infinitesimal time dt, or (2) at fixed time t after an infinitesimal rigid rotation dt Ω about some rotation center in configuration space {x}, and likewise an infinitesimal rigid rotation dt ω about some rotation center in velocity space {v}.

3.1. Rotation about the origin

First, we start with the simplest case. The rotation center is assumed to be fixed and located at the origin. We generalize the derivation used by TW in their pioneer paper to more cases. The function f(x, v, t) invariance, either after an infinitesimal time displacement or after an infinitesimal solid rotation about the origin, reads

$$\begin{array}{c}\hfill f(\mathit{x},\mathit{v},t+\mathrm{d}t)=f(\mathit{x}+\mathrm{d}t\mathbf{\Omega }\wedge \mathit{x},\mathit{v}+\mathrm{d}t\mathit{\omega }\wedge \mathit{v},t).\end{array}$$(3)

A priori, the pattern speeds Ω and ω do not need to be equal. For the sake of rigor, we verify in Appendix B the mere existence of such functions by displaying an explicit class of these.

Expanding both sides of the equation to first order in time leads to a linear equation for Ω and ω that can be reformulated more conveniently by using the dot and cross products property a ⋅ (b ∧ c) = (c ∧ a)⋅b,

$$\begin{array}{c}\hfill {\partial }_{\mathrm{t}}f={\partial }_{\mathit{x}}f·(\mathbf{\Omega }\wedge \mathit{x})+{\partial }_{\mathit{v}}f·(\mathit{\omega }\wedge \mathit{v})=\mathit{x}\wedge {\partial }_{\mathit{x}}f·\mathbf{\Omega }+\mathit{v}\wedge {\partial }_{\mathit{v}}f·\mathit{\omega }.\end{array}$$(4)

If the direction unit vector n of the rotation axis is known, Ω = Ωn, and the two frequencies are the same, ω = Ω, we can then attempt to solve this linear equation for Ω provided the partial derivatives of f are known at a given phase space point (x, v) at time t. Rotating the axes such that n = (0, 0, 1), we find

$$\begin{array}{c}\hfill \mathrm{\Omega }=\frac{{\partial }_{t}f}{x{\partial }_{y}f-y{\partial }_{x}f+v_x{\partial }_{v_y}-v_y{\partial }_{v_x}f}.\end{array}$$(5)

If n is not known and ω = Ω, we need n ≥ 3 independent positions (x_i, v_i), i = 1, 2, …n to solve by linear optimization, such as the least-squares method, the linear system of equations

$$\begin{array}{c}\hfill \underset{n\times 3\text{matrix}}{\underbrace{\left(\begin{array}{c}{\mathit{x}}_1\wedge {\partial }_{\mathit{x}}{f}_1+{\mathit{v}}_1\wedge {\partial }_{\mathit{v}}{f}_1\\ ⋮\\ {\mathit{x}}_n\wedge {\partial }_{\mathit{x}}{f}_n+{\mathit{v}}_n\wedge {\partial }_{\mathit{v}}{f}_n\end{array}\right)}}\underset{3\times 1\text{vector}}{\underbrace{\mathbf{\Omega }}}\approx \underset{n\times 1\text{vector}}{\underbrace{\left(\begin{array}{c}{\partial }_t{f}_1\\ ⋮\\ {\partial }_t{f}_n\end{array}\right)}},\end{array}$$(6)

where f_i means f(x_i, v_i, t). In the presence of noise, no reliable solution exists if the n × 3 left-hand side matrix is degenerate, however.

The least-squares solution of a linear system Ax ≈ b, where A is the n × m matrix (n > m), and b the n × 1 vector, is the unique solution for the m vector x minimizing the Euclidian L₂ norm ||Ax − b||₂, and, when A is degenerate, the minimum ||x||₂. In addition to the traditional least-squares method, methods that minimize other norms are available, such as the L₁ norm, which may be less sensitive to outlier data points (Boyd & Vandenbergh 2004). In contrast to the L₂ norm, however, a L₁ norm solution may not always exist.

If ω ≠ Ω, we need n ≥ 6 independent positions (x_i, v_i), i = 1, 2, …n to solve the linear system as above,

$$\begin{array}{c}\hfill \underset{n\times 6\text{matrix}}{\underbrace{\left(\begin{array}{cc}{\mathit{x}}_1\wedge {\partial }_{\mathit{x}}{f}_1& {\mathit{v}}_1\wedge {\partial }_{\mathit{v}}{f}_1\\ ⋮& ⋮\\ {\mathit{x}}_n\wedge {\partial }_{\mathit{x}}{f}_n& {\mathit{v}}_n\wedge {\partial }_{\mathit{v}}{f}_n\end{array}\right)}}\underset{6\times 1\text{vector}}{\underbrace{\left(\begin{array}{c}\mathbf{\Omega }\\ \mathit{\omega }\end{array}\right)}}\approx \underset{n\times 1\text{vector}}{\underbrace{\left(\begin{array}{c}{\partial }_t{f}_1\\ ⋮\\ {\partial }_t{f}_n\end{array}\right)}}.\end{array}$$(7)

In these linear optimization problems, each line of Eqs. (6) and (7) can always be multiplied left and right by a weight w_i to weight the data points differently. For example, the weight can be inversely proportional to the error σ_i expected on the ith data point, $w_i={\sigma }_i^{-1}$. For the L₂ norm, this weight is equivalent to performing a χ² minimization.

3.2. Incompressible function f

If f conserves phase-space density, such as for Boltzmann’s collisionless equation, where f(x, v, t) is the phase-space distribution function,

$$\begin{array}{c}\hfill {\partial }_{t}f+\mathit{v}·{\partial }_{\mathit{x}}f+\mathit{a}·{\partial }_{\mathit{v}}f=0,\end{array}$$(8)

where $\mathit{a}\equiv \dot{\mathit{v}}$ is the acceleration, then ∂_tf can be eliminated.

When ω = Ω, we obtain

$$\begin{array}{c}\hfill \left(\begin{array}{c}{\mathit{x}}_1\wedge {\partial }_{\mathit{x}}{f}_1+{\mathit{v}}_1\wedge {\partial }_{\mathit{v}}{f}_1\\ ⋮\\ {\mathit{x}}_n\wedge {\partial }_{\mathit{x}}{f}_n+{\mathit{v}}_n\wedge {\partial }_{\mathit{v}}{f}_n\end{array}\right)\mathbf{\Omega }\approx -\left(\begin{array}{c}{\mathit{v}}_1·{\partial }_{\mathit{x}}{f}_1+{\mathit{a}}_1·{\partial }_{\mathit{v}}{f}_1\\ ⋮\\ {\mathit{v}}_n·{\partial }_{\mathit{x}}{f}_n+{\mathit{a}}_n·{\partial }_{\mathit{v}}{f}_n\end{array}\right),\end{array}$$(9)

and when ω ≠ Ω,

$$\begin{array}{c}\hfill \left(\begin{array}{cc}{\mathit{x}}_1\wedge {\partial }_{\mathit{x}}{f}_1& {\mathit{v}}_1\wedge {\partial }_{\mathit{v}}{f}_1\\ ⋮& ⋮\\ {\mathit{x}}_n\wedge {\partial }_{\mathit{x}}{f}_n& {\mathit{v}}_n\wedge {\partial }_{\mathit{v}}{f}_n\end{array}\right)\left(\begin{array}{c}\mathbf{\Omega }\\ \mathit{\omega }\end{array}\right)\approx -\left(\begin{array}{c}{\mathit{v}}_1·{\partial }_{\mathit{x}}{f}_1+{\mathit{a}}_1·{\partial }_{\mathit{v}}{f}_1\\ ⋮\\ {\mathit{v}}_n·{\partial }_{\mathit{x}}{f}_n+{\mathit{a}}_n·{\partial }_{\mathit{v}}{f}_n\end{array}\right).\end{array}$$(10)

Thus, when f depends on v, we need to also know the acceleration (or force) field a, which is typically minus the gradient of the gravitational potential Φ(x). When f does not depend on v, all the ∂_vf_i terms cancel, a is not needed, and as it should, ω then remains undetermined.

3.3. Mass-conserving function f(x, t)

Likewise, when f conserves mass in space, for instance, the mass density in a fluid, often called ρ, and u = u(x, t) is the velocity field, ∂_tf can be eliminated by the continuity equation

$$\begin{array}{c}\hfill {\partial }_{t}f+{\partial }_{\mathit{x}}·(f\mathit{u})=0.\end{array}$$(11)

Because u is a function of x, it rotates in space at the same speed as in velocity space, so that Ω = ω. Determining Ω then requires knowing only the spatial gradients of f and u at time t.

3.4. Translation and rotation about a known center

When f is invariant only when it rotates about the center at frequencies Ω and ω, and when in addition, the center moves at velocity V and acceleration A, Eq. (3) is modified as

$$\begin{array}{c}\hfill f(\mathit{x},\mathit{v},t+dt)=f(\mathit{x}+dt(\mathit{V}+\mathbf{\Omega }\wedge \mathit{x}),\mathit{v}+dt(\mathit{A}+\mathit{\omega }\wedge \mathit{v}),t).\end{array}$$(12)

Expanding both sides of the equation as before to first order in time leads to a linear equation system for V, A, Ω, and ω,

$$\begin{array}{c}\hfill \underset{n\times 12\text{matrix}}{\underbrace{\left(\begin{array}{cccc}{\mathit{x}}_1\wedge {\partial }_{\mathit{x}}{f}_1& {\mathit{v}}_1\wedge {\partial }_{\mathit{v}}{f}_1& {\partial }_{\mathit{x}}{f}_1& {\partial }_{\mathit{v}}{f}_1\\ ⋮& ⋮& ⋮& ⋮\\ {\mathit{x}}_n\wedge {\partial }_{\mathit{x}}{f}_n& {\mathit{v}}_n\wedge {\partial }_{\mathit{v}}{f}_n& {\partial }_{\mathit{x}}{f}_n& {\partial }_{\mathit{v}}{f}_n\end{array}\right)}}\underset{12\times 1\text{vector}}{\underbrace{\left(\begin{array}{c}\mathbf{\Omega }\\ \mathit{\omega }\\ \mathit{V}\\ \mathit{A}\end{array}\right)}}\approx \underset{n\times 1\text{vector}}{\underbrace{\left(\begin{array}{c}{\partial }_t{f}_1\\ ⋮\\ {\partial }_t{f}_n\end{array}\right)}}.\end{array}$$(13)

As above, we can replace the ∂_tf_i terms if f follows a conservation law. In principle, we can therefore quantify by how much the system center of rotation moves or accelerates with the same data required to solve Eq. (10) by solving Eq. (13) for Ω, V, and A with at least n ≥ 12 data points by linear optimization.

Because Eq. (13) is dimensionally inhomogeneous, it may be useful to scale the left-hand side matrix columns and vector lines by appropriate constants. For example, when we use length and velocity scale factors R₀ and V₀, the system becomes

$$\begin{array}{c}\hfill \left(\begin{array}{cccc}{\mathit{x}}_1\wedge {\partial }_{\mathit{x}}{f}_1& {\mathit{v}}_1\wedge {\partial }_{\mathit{v}}{f}_1& R_0{\partial }_{\mathit{x}}{f}_1& V_0{\partial }_{\mathit{v}}{f}_1\\ ⋮& ⋮& ⋮& ⋮\\ {\mathit{x}}_n\wedge {\partial }_{\mathit{x}}{f}_n& {\mathit{v}}_n\wedge {\partial }_{\mathit{v}}{f}_n& R_0{\partial }_{\mathit{x}}{f}_n& V_0{\partial }_{\mathit{v}}{f}_n\end{array}\right)\left(\begin{array}{c}\mathbf{\Omega }\\ \mathit{\omega }\\ \mathbf{{\rm Y}}\\ \mathbf{\Lambda }\end{array}\right)\approx \left(\begin{array}{c}{\partial }_t{f}_1\\ ⋮\\ {\partial }_t{f}_n\end{array}\right),\end{array}$$(14)

and V = R₀Υ, A = V₀Λ.

3.5. Translation and rotation about an unknown center

When the center of rotation is unknown both in space and velocity space, the true centers of rotation would be x₀ and v₀, respectively. The invariant f would follow the following constraint to first order in dt,

$$\begin{array}{c}\hfill f(\mathit{x},\mathit{v},t+dt)=f(\mathit{x}+dt(\mathit{V}+\mathbf{\Omega }\wedge (\mathit{x}-{\mathit{x}}_0)),\mathit{v}+dt(\mathit{A}+\mathit{\omega }\wedge (\mathit{v}-{\mathit{v}}_0)),t).\end{array}$$(15)

Expanding both sides of the equation as before to first order in time leads to a quadratic (hyperbolic) equation system for x₀, v₀, Ω, and ω, and linear for V, A,

$$\begin{array}{cc}\hfill \underset{n\times 12\text{matrix}}{\underbrace{\left(\begin{array}{cccc}({\mathit{x}}_1-{\mathit{x}}_0)\wedge {\partial }_{\mathit{x}}{f}_1& ({\mathit{v}}_1-{\mathit{v}}_0)\wedge {\partial }_{\mathit{v}}{f}_1& {\partial }_{\mathit{x}}{f}_1& {\partial }_{\mathit{v}}{f}_1\\ ⋮& ⋮& ⋮& ⋮\\ ({\mathit{x}}_n-{\mathit{x}}_0)\wedge {\partial }_{\mathit{x}}{f}_n& ({\mathit{v}}_n-{\mathit{v}}_0)\wedge {\partial }_{\mathit{v}}{f}_n& {\partial }_{\mathit{x}}{f}_n& {\partial }_{\mathit{v}}{f}_n\end{array}\right)}}\underset{12\times 1\text{vector}}{\underbrace{\left(\begin{array}{c}\mathbf{\Omega }\\ \mathit{\omega }\\ \mathit{V}\\ \mathit{A}\end{array}\right)}}& \approx \underset{n\times 1\text{vector}}{\underbrace{\left(\begin{array}{c}{\partial }_t{f}_1\\ ⋮\\ {\partial }_t{f}_n\end{array}\right)}}.\hfill \end{array}$$(16)

Rearranging the terms, we can separate the linear and quadratic parts,

$$\begin{array}{cc}\hfill \underset{n\times 6\text{matrix}}{\underbrace{\left(\begin{array}{cc}\mathbf{\Omega }\wedge {\partial }_{\mathit{x}}{f}_1& \mathit{\omega }\wedge {\partial }_{\mathit{v}}{f}_1\\ ⋮& ⋮\\ \mathbf{\Omega }\wedge {\partial }_{\mathit{x}}{f}_n& \mathit{\omega }\wedge {\partial }_{\mathit{v}}{f}_n\end{array}\right)}}\underset{6\times 1\text{vector}}{\underbrace{\left(\begin{array}{c}{\mathit{x}}_0\\ {\mathit{v}}_0\end{array}\right)}}& \approx \underset{n\times 1\text{vector}}{\underbrace{\left(\begin{array}{c}{\partial }_t{f}_1\\ ⋮\\ {\partial }_t{f}_n\end{array}\right)}}-\underset{n\times 12\text{matrix}}{\underbrace{\left(\begin{array}{cccc}{\mathit{x}}_1\wedge {\partial }_{\mathit{x}}{f}_1& {\mathit{v}}_1\wedge {\partial }_{\mathit{v}}{f}_1& {\partial }_{\mathit{x}}{f}_1& {\partial }_{\mathit{v}}{f}_1\\ ⋮& ⋮& ⋮& ⋮\\ {\mathit{x}}_n\wedge {\partial }_{\mathit{x}}{f}_n& {\mathit{v}}_n\wedge {\partial }_{\mathit{v}}{f}_n& {\partial }_{\mathit{x}}{f}_n& {\partial }_{\mathit{v}}{f}_n\end{array}\right)}}\underset{12\times 1\text{vector}}{\underbrace{\left(\begin{array}{c}\mathbf{\Omega }\\ \mathit{\omega }\\ \mathit{V}\\ \mathit{A}\end{array}\right)}}.\hfill \end{array}$$(17)

An iterative approach to determine a least-squares solution for this system might be to start by guessing initial values for x₀, v₀, and then to solve Eq. (16) for Ω, ω, V, and A by linear least-squares. Using these values, we could then solve for x₀, v₀ Eq. (17) by linear least-squares. Iterating this process a number of times might converge to a solution. Otherwise, quadratic or nonlinear least-squares algorithms could be required. In general, we expect either no, one, or several real solutions with quadratic problems.

4. Local and regional 3D Tremaine-Weinberg methods (LTW and RTW methods)

The TW method is based on a particular case of invariant function, as discussed in Sect. 3, where the surface density

$$\begin{array}{c}\hfill \mathrm{\Sigma }(x,y,t)={\displaystyle \int }\mathrm{d}z{\displaystyle \int }{\mathrm{d}}^{3}vf(x,y,z,\mathit{v},t)\end{array}$$(18)

is a projection of a thin distribution function f in the galactic plane of a tracer assumed to verify the continuity equation ∂_tΣ + ∂_x ⋅ (Σu) = 0, where u(x, y) is the local mass-weighted average velocity in the plane x − y. Because Σ does not depend on v, this method does not require knowledge of the force field. In the original method, the direction of the rotation axis n is assumed to be known. Using Eq. (9) at a single position, and integrating over all velocities,

$$\begin{array}{c}\hfill (\mathit{x}\wedge {\partial }_{\mathit{x}}\mathrm{\Sigma })·\mathit{n}{\mathrm{\Omega }}_{\mathrm{LTW}}=-{\partial }_{\mathit{x}}·\left(\mathrm{\Sigma }\mathit{u}\right).\end{array}$$(19)

In principle, a single local calculation at x is therefore sufficient to find Ω_LTW. Because this calculation involves gradients, however, observational errors may be an issue, as taking the gradient of data, for example, by finite differences, typically amplifies errors. Furthermore, even when precise gradients are available, the coefficient x ∧ ∂_xΣ in Eq. (19) may vanish (e.g., when the azimuthal density profile is an extremum, or more generally, when it is tangent to a circle centered at the origin), so that large errors are contributed by all the regions that are locally rotationally symmetric. This occurs, for example, close to the main axes of a bar, but also in the disk. By weighting each data point line in Eq. (19) by a coefficient w_i proportional to the error inverse, we can thus achieve a χ² minimization.

The original TW method attempts to minimize noise by using weighted averages of Eq. (19) along lines. In these line integrals, the high-noise segments mostly contribute to increasing errors, the above-noise signal segments alone do contribute to finding Ω. Extending the integration interval to infinity and performing integration by parts avoids formally dealing with gradients in the integrals, but this does not remove the measurement noise, and it necessarily includes contributions outside the relevant region and in this way, it may add systematic errors in case of multiple pattern speeds. Since N-body bars and real galaxies have different pattern speeds in their bar(s) and surrounding spiral arms (Sellwood & Sparke 1988; Wu et al. 2016, 2018), and double-bar systems are frequent (Erwin 2004), it is therefore necessary to consider local or regional evaluations of pattern speeds.

Local methods use data gradients, however, which, if done without attention, are prone to amplifying noise. In the literature, many approaches exist for evaluating gradients of noisy data. For the systems considered here, methods that smooth data locally have to be preferred over methods that use global characteristics of noise. One of the oldest local method is the Savitzky-Golay filter (Savitzky & Golay 1964). Spectral methods that smooth the gradient operator instead of smoothing data are intrinsically nonlocal. In astrophysical particle systems, the smoothed-particle hydrodynamics technique (Gingold & Monaghan 1977; Lucy 1977) is frequently used to evaluate the density gradients and other physical quantities attached to particles. A smoothing kernel convolved spatially over the particle data provides smooth density and velocity functions. In numerical mathematics, more advanced methods have been proposed that not only provide accurate gradients, but also preserve discontinuities (e.g., Chartrand 2011, 2017). We refer to Van Breugel et al. (2020) for a recent short review of the problem. In view of the general need in science and engineering to extract derivatives of data, the question will continue to progress in the following years. Below we assume that the problem of taking gradients of data has been sufficiently well mastered.

4.1. Local method (LTW)

An alternative approach to TW is to choose a number of selected locations with a high signal-to-noise ratio instead of lines, which consequently excludes the regions in which azimuthal gradients are too small, below the observational noise.

When we do not assume n as known, we can determine the vector Ω_LTW by least-squares using three or more distinct positions x_i, i = 1, …n, for the 3D space density, and solve the equation system by linear optimization,

$$\begin{array}{c}\hfill \left(\begin{array}{c}{\mathit{x}}_1\wedge {\partial }_{\mathit{x}}\rho ({\mathit{x}}_1)\\ ⋮\\ {\mathit{x}}_n\wedge {\partial }_{\mathit{x}}\rho ({\mathit{x}}_n)\end{array}\right){\mathbf{\Omega }}_{\mathrm{LTW}}\approx -\left(\begin{array}{c}{\partial }_{\mathit{x}}·\left(\rho ({\mathit{x}}_1)\mathit{u}({\mathit{x}}_1)\right)\\ ⋮\\ {\partial }_{\mathit{x}}·\left(\rho ({\mathit{x}}_n)\mathit{u}({\mathit{x}}_n)\right)\end{array}\right),\end{array}$$(20)

where

$$\begin{array}{c}\hfill \rho (\mathit{x},t)={\displaystyle \int }{\mathrm{d}}^{3}vf(\mathit{x},\mathit{v},t),\end{array}$$(21)

and

$$\begin{array}{c}\hfill \rho (\mathit{x},t)\mathit{u}(\mathit{x},t)={\displaystyle \int }{d}^{3}v\mathit{v}f(\mathit{x},\mathit{v},t).\end{array}$$(22)

In principle, if these positions provide a nondegenerate set of linear equations, the full rotation vector Ω_LTW can be retrieved.

Thus the LTW method consists in choosing a finite number of locations well in which the azimuthal density gradient ∝|x∧∂_xρ(x)| has a high signal-to-noise ratio. In bars, the maximum signal is to be sought in the off-axis regions, a butterfly shape, determined by preferentially accepting the data points where |x∧∂_xρ(x)| is above the noise. Typically, locations near the principal axes of a bar should be rejected: because the density gradient is radial, these regions add noise to the linear system. Because this involves measuring the gradient of ρ and ρu, it can be used in N-body bars where errors are under control.

4.2. Regional method (RTW)

However, if gradients are too noisy, we can proceed as in TW by integrating data over a volume. We do this here in 3D. Instead of the TW infinite line in the galaxy plane, we generalize by integrating over an arbitrary finite contiguous volume V bounded by an oriented surface S in one or several favorable high signal-to-noise regions.

The 3D local equation for the mass density equation ρ(x) in a velocity file $\overline{\mathit{v}}(\mathit{x})$ is

$$\begin{array}{c}\hfill (\mathit{x}\wedge {\partial }_{\mathit{x}}\rho )·{\mathbf{\Omega }}_{\mathrm{RTW}}=-{\partial }_{\mathit{x}}·\left(\rho \overline{\mathit{v}}\right).\end{array}$$(23)

Integrating over V, noting that ∇ρ ∧ x = ∇ ∧ (ρx), and using the divergence and curl theorems of vector calculus¹, we obtain an integral over the surface S on both sides of the equation,

$$\begin{array}{c}\hfill {\displaystyle ({\int }_S\rho \mathit{x}\wedge \mathrm{d}\mathit{S})·{\mathbf{\Omega }}_{\mathit{RTW}}\approx {\int }_S\rho \overline{\mathit{v}}·\mathrm{d}\mathit{S}.}\end{array}$$(24)

These integrals amount to the total mass flux through surface S, one expressed with the supposed unknown uniform rotation of the system², the other by the measured mass flux parallel to S. When S is chosen poorly, this flux may vanish in all cases, which is thus a degenerate configuration to be avoided. For example, a sphere centered on the origin sees all x being parallel to the surface elements dS, so (ρx)∧dS = 0 for any ρ, Eq. (24) reduces to 0 ⋅ Ω_RTW = 0.

The RTW method can be considered to be a method that replaces the gradient of data evaluation at a point by directly integrating data over a finite surface. The optimal volume size is determined by the noise in the available data. Since only the information present on the surface of V is used, there is no interest in taking overly large regions.

4.3. Multiple regional method (MRTW)

In principle, we can therefore combine the two previous approaches by a multiple RTW method (MRTW), for which several distinct signal-rich regions inside volumes bounded by surfaces S₁, S₂, …S_n are used, allowing us to solve a series of Eq. (24) by linear optimization for Ω_MRTW,

$$\begin{array}{c}\hfill \left(\begin{array}{c}{\int }_{S_1}\rho \mathit{x}\wedge \mathrm{d}{\mathit{S}}_1\\ ⋮\\ {\int }_{S_n}\rho \mathit{x}\wedge \mathrm{d}{\mathit{S}}_n\end{array}\right){\mathbf{\Omega }}_{\mathrm{MRTW}}\approx \left(\begin{array}{c}{\int }_{S_1}\rho \overline{\mathit{v}}·\mathrm{d}{\mathit{S}}_1\\ ⋮\\ {\int }_{S_n}\rho \overline{\mathit{v}}·\mathrm{d}{\mathit{S}}_n\end{array}\right).\end{array}$$(25)

For the Milky Way, using full 6D astrometric data from Gaia in volume-limited regions around the Sun where extinction is under control, we might in principle attempt to evaluate the pattern speed of local spiral arms, and possibly more.

5. Jacobi integral, or potential method

Here we develop a method that does not require density gradients, but requires the individual kinematic data of test particles, or stars. This method was already briefly described in Wu et al. (2018).

In a gravitational mean potential Φ that rotates uniformly with angular frequency Ω, a test particle moving in this potential always has the Jacobi integral H_J as the global integral of motion (Binney & Tremaine 2008),

$$\begin{array}{c}\hfill H_{\mathrm{J}}=\frac{1}{2}{\mathit{v}}^2(t)+\mathrm{\Phi }(\mathit{x}(t))-\mathbf{\Omega }·(\mathit{x}(t)\wedge \mathit{v}(t))=e(t)-\mathbf{\Omega }·\mathit{l}(t),\end{array}$$(26)

where x(t) is the position of the particle, and v(t) the velocity in the inertial frame instantaneously aligned with the rotating frame at time t. The quantity $e(t)=\frac{1}{2}{\mathit{v}}^2(t)+\mathrm{\Phi }(\mathit{x}(t))$ is the particle specific energy in the inertial frame, and l(t) = x(t)∧v(t) is its specific angular momentum, also in the inertial frame.

5.1. Pattern speed knowing its direction

When we know the direction of the rotation axis n, Ω = Ω n, we can in principle determine Ω from the motion of a single particle because along an orbit

$$\begin{array}{c}\hfill 0={\dot{H}}_{\mathrm{J}}(\mathit{x},\mathit{v})=\dot{e}-\mathrm{\Omega }\mathit{n}·\dot{\mathit{l}}=\mathit{v}·\mathit{a}+\dot{\mathrm{\Phi }}(\mathit{x})-\mathrm{\Omega }\mathit{n}·(\mathit{x}\wedge \mathit{a}),\end{array}$$(27)

that is,

$$\begin{array}{c}\hfill (\mathit{x}\wedge \mathit{a})·\mathit{n}\mathrm{\Omega }=\mathit{v}·\mathit{a}+\dot{\mathrm{\Phi }}(\mathit{x}).\end{array}$$(28)

In principle, provided the torque (x∧a) is nonzero and not orthogonal to n, a single particle has the information about Ω in the values of its instantaneous position, velocity, acceleration and the time-derivative of its potential along the trajectory, $\dot{\mathrm{\Phi }},$ at any time along its trajectory. At some locations along a trajectory, the tangential acceleration may become low (the torque term x ∧ a in Eq. (28) almost vanishes). Then Ω is ill defined. In nonaxisymmetric systems, these locations are exceptional. Thus the best moments for determining Ω along a trajectory occur when angular momentum varies the most, that is, when the torque is strong.

Alternatively, the above derivation leading to Eq. (28) can be obtained without invoking Jacobi’s integral by assuming that the potential Φ is an invariant function in some rotating frame, and using the equations developed in Sect. 3. Replacing f by Φ in Eq. (9), and noting that $\dot{\mathrm{\Phi }}={\partial }_t\mathrm{\Phi }+\mathit{v}·{\partial }_{\mathit{x}}\mathrm{\Phi }$, so ${\partial }_t\mathrm{\Phi }=\mathit{v}·\mathit{a}+\dot{\mathrm{\Phi }}$, we retrieve Eq. (28), with nΩ replaced by Ω.

5.2. Pattern speed vector

Using three or more independent particles, we can form the linear system

$$\begin{array}{c}\hfill \underset{n\times 3\text{matrix}}{\underbrace{\left(\begin{array}{c}{\mathit{x}}_1\wedge {\mathit{a}}_1\\ ⋮\\ {\mathit{x}}_n\wedge {\mathit{a}}_n\end{array}\right)}}\mathbf{\Omega }\approx \underset{n\times 1\text{vector}}{\underbrace{\left(\begin{array}{c}{\mathit{v}}_1·{\mathit{a}}_1+\dot{\mathrm{\Phi }}({\mathit{x}}_1)\\ ⋮\\ {\mathit{v}}_n·{\mathit{a}}_n+\dot{\mathrm{\Phi }}({\mathit{x}}_n)\end{array}\right)}},\end{array}$$(29)

which can be solved for Ω by linear optimization. Furthermore, as in Eq. (13), we can extend the search for the motion of the system rotation center at speed V,

$$\begin{array}{c}\hfill \underset{n\times 6\text{matrix}}{\underbrace{\left(\begin{array}{cc}{\mathit{x}}_1\wedge {\mathit{a}}_1& {\mathit{a}}_1\\ ⋮& ⋮\\ {\mathit{x}}_n\wedge {\mathit{a}}_n& {\mathit{a}}_n\end{array}\right)}}\underset{6\times 1\text{matrix}}{\underbrace{\left(\begin{array}{c}\mathbf{\Omega }\\ \mathit{V}\end{array}\right)}}\approx \underset{n\times 1\text{vector}}{\underbrace{\left(\begin{array}{c}{\mathit{v}}_1·{\mathit{a}}_1+\dot{\mathrm{\Phi }}({\mathit{x}}_1)\\ ⋮\\ {\mathit{v}}_n·{\mathit{a}}_n+\dot{\mathrm{\Phi }}({\mathit{x}}_n)\end{array}\right)}}.\end{array}$$(30)

The particularity of this method is that $\dot{\mathrm{\Phi }}$ must be determined. In usual N-body simulation codes, all the other values except for $\dot{\mathrm{\Phi }}$ are calculated at each time step. With a small increase in computation cost, the force algorithm can also calculate Φ, so that $\dot{\mathrm{\Phi }}(t)$ can be approximated to second order in Δt by finite difference in time and two additional potential calculations,

$$\begin{array}{c}\hfill \dot{\mathrm{\Phi }}(\mathit{x}(t))=\frac{\mathrm{\Phi }(\mathit{x}(t+\mathrm{\Delta }t))-\mathrm{\Phi }(\mathit{x}(t-\mathrm{\Delta }t))}{2\mathrm{\Delta }t}+\mathcal{O}(\mathrm{\Delta }{t}^2)\approx \frac{\mathrm{\Phi }(\mathit{x}(t)+\mathit{v}(t)\mathrm{\Delta }t)-\mathrm{\Phi }(\mathit{x}(t)-\mathit{v}(t)\mathrm{\Delta }t)}{2\mathrm{\Delta }t},\end{array}$$(31)

where Δt is a positive time step. In a pure N-body system, and similarly, for other potential solvers, the time derivative of the potential can be calculated directly with little additional expense simultaneously with the potential and forces by

$$\begin{array}{c}\hfill \dot{\mathrm{\Phi }}({\mathit{x}}_j(t))=\frac{\mathrm{d}}{\mathrm{d}t}(-G{\displaystyle \sum _{i\ne j}}\frac{m_i}{|{\mathit{x}}_i(t)-{\mathit{x}}_j(t)|})=G{\displaystyle \sum _{i\ne j}{m}_i\frac{({\mathit{x}}_i(t)-{\mathit{x}}_j(t))·({\mathit{v}}_i(t)-{\mathit{v}}_j(t))}{|{\mathit{x}}_i(t)-{\mathit{x}}_j{(t)|}^3}.}\end{array}$$(32)

In an N-body system, however, the potential is clearly never strictly time independent. It differs from a mean potential because nearby particles eventually contribute to strong fluctuations of the potential over the background of slowly varying contributions from particles that are far away. The distribution of the found Ω obtained from a set of selected particles can help to quantify the force fluctuations from the nearest particles with respect to the smoother mean field force: if the distribution of found Ω is broad, it means that there is no single frame with a uniformly rotating potential, or that the fast local fluctuating forces dominate the more constant force from the mean rotating potential. In a smooth mass distribution, the contribution to force fluctuations from nearby particles becomes particularly important near the potential minimum (near the center), or everywhere in extended nearly homogeneous distributions like in cosmology.

6. Second-order moment methods

6.1. Moment of inertia tensor

In mechanics, the polar moment of inertia tensor I is a simple 3 × 3 tensor characterizing the second moments of the mass distribution with respect to the origin. It is set at the location at which the center of mass and center of rotation coincide. It can be seen as the covariance matrix weighted by the mass of the particle positions. For an N particle distribution, it is defined as

$$\begin{array}{c}\hfill \mathit{I}\equiv {\mathit{x}}^T\mathit{M}\mathit{x}=\left(\begin{array}{ccc}{\sum }_i{m}_i{x}_i^2& {\sum }_i{m}_i{x}_i{y}_i& {\sum }_i{m}_i{x}_i{z}_i\\ {\sum }_i{m}_i{x}_i{y}_i& {\sum }_i{m}_i{y}_i^2& {\sum }_i{m}_i{y}_i{z}_i\\ {\sum }_i{m}_i{x}_i{z}_i& {\sum }_i{m}_i{y}_i{z}_i& {\sum }_i{m}_i{z}_i^2\end{array}\right),\end{array}$$(33)

where x is the N × 3 matrix composed of the stacked particles in Cartesian coordinates x_i,

$$\begin{array}{c}\hfill \mathit{x}\equiv \left(\begin{array}{ccc}{x}_1& y_1& z_1\\ ⋮& ⋮& ⋮\\ x_N& y_N& z_N\end{array}\right),\end{array}$$(34)

and M is the N × N mass diagonal matrix, with its diagonal filled with the particle masses m_i,

$$\begin{array}{c}\hfill \mathit{M}\equiv \left(\begin{array}{ccc}{m}_1& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & m_N\end{array}\right).\end{array}$$(35)

When only the average particle distribution is rigid except for a uniform rotation about an axis, we can still use the axial moment of inertia tensor in use in solid body mechanics (e.g., Goldstein 1971) even though each particle moves nonrigidly with respect to each other to relate the rotation of the particle distribution to its constant total angular momentum vector L = ∑_im_ix_i ∧ v_i and its axial moment of inertia tensor, which is given by

$$\begin{array}{c}\hfill \mathbb{I}\equiv \mathrm{Tr}(\mathit{I})\mathbf{1}-\mathit{I}=\left(\begin{array}{ccc}{\sum }_i{m}_i(y_i^2+z_i^2)& -{\sum }_i{m}_i{x}_i{y}_i& -{\sum }_i{m}_i{x}_i{z}_i\\ -{\sum }_i{m}_i{x}_i{y}_i& {\sum }_i{m}_i(z_i^2+x_i^2)& -{\sum }_i{m}_i{y}_i{z}_i\\ -{\sum }_i{m}_i{x}_i{z}_i& -{\sum }_i{m}_i{y}_i{z}_i& {\sum }_i{m}_i(x_i^2+y_i^2)\end{array}\right).\end{array}$$(36)

Then the angular speed Ω_s of the solid body is related to L by

$$\begin{array}{c}\hfill \mathit{L}=\mathbb{I}{\mathbf{\Omega }}_{\text{s}}.\end{array}$$(37)

As the system rotates, I and therefore 𝕀 vary: the non-diagonal terms, for example, vanish when the Cartesian coordinates are aligned with the principal axes of the mass distributions. Ω_s must still be constant, however, because L is conserved and its value should not change under a rotation of the coordinates about L. This means that knowing the particle distribution and its total angular momentum allows us to find a unique Ω_s that describes the instantaneous rotation rate as if the particles were moving as a solid. When the particle distribution does not rotate rigidly, then Ω_s can clearly change, still respecting Eq. (37).

The only technical problem is to find a well-behaved invert of 𝕀 for all possible configurations. Pfenniger (2019) gave an explicit procedure for finding the inverse 𝕀⁻¹, or more generally, the pseudo-inverse 𝕀^† in degenerate cases,

$$\begin{array}{c}\hfill {\mathbf{\Omega }}_{\text{s}}={\mathbb{I}}^{†}\mathit{L}.\end{array}$$(38)

For a body with internal streaming, however, such as a bar in a galaxy, Ω_s is not necessarily the pattern speed because its value does not depend on any feature of the mass distribution and also works with axisymmetric systems rotating about the symmetry axis.

The polar moment of inertia tensor I described above is a simple 3 × 3 tensor characterizing the second moments of the mass distribution with respect to the origin, set at the center of mass. If the system globally keeps its shape over time while rotating uniformly, the eigenvalues of I remain constant, and its rotating eigenvectors may be used to determine the pattern speed.

6.1.1. Pattern speed of the 3D system vector

In principle, the real eigenvalues and eigenvectors can be calculated analytically for a 3 × 3 symmetric matrix, but the expressions are huge because the eigenvalues are solutions of a cubic polynomial. It is more practical to use well-established numerical methods to determine them (Anderson et al. 1999; Trefethen & Bau 1997).

For a symmetric positive definite matrix A, the eigenvalue problem, which consists of finding the orthogonal eigenvector matrix V and the diagonal eigenvalue matrix Λ for which A = VΛV^T, is identical to the singular value decomposition (SVD) problem for a general (nonsquare) matrix: A = USV^T-pagination, where S is the diagonal matrix of (real) singular values and U and V are two orthogonal matrices. For positive definite symmetric matrices U = V, the eigenvalue problem is thus equivalent to the SVD problem. In contrast to standard SVD algorithms, however, common eigenvalue–eigenvector algorithms do not sort eigenvalues and eigenvectors in a standard way because they may be complex or real negative. It is therefore more convenient to use the SVD algorithm. For small matrices, the difference in computational cost is negligible.

After decomposing I with SVD: I = USU^T-pagination, we obtain a real orthogonal rotation matrix U made of the three orthogonal unit vectors n_i, where i = 1, 2, 3, U = [n₁, n₂, n₃], which gives the direction of the eigenvectors, and a diagonal real matrix S with singular values s₁ ≥ s₂ ≥ s₃ ≥ 0 in decreasing order, giving the second-order moment amplitudes in the respective directions. If the second-order moments are invariant, S does not depend on time, but $\dot{U}$, the time derivative of U, is a matrix giving the speed of the eigenvectors.

To evaluate $\dot{U}$, we first need İ(t), which is obtained using the particle positions x and velocities V (assuming mass conservation, Ṁ = 0),

$$\begin{array}{c}\hfill \dot{\mathit{I}}={\dot{\mathit{x}}}^T\mathit{M}\mathit{x}+{\mathit{x}}^T\mathit{M}\dot{\mathit{x}}={\mathit{V}}^T\mathit{M}\mathit{x}+{\mathit{x}}^T\mathit{M}\mathit{V}=\left(\begin{array}{ccc}2{\sum }_i{m}_i{x}_i{v}_{x_i}& {\sum }_i{m}_i(x_i{v}_{y_i}+y_i{v}_{x_i})& {\sum }_i{m}_i(x_i{v}_{z_i}+z_i{v}_{x_i})\\ {\sum }_i{m}_i(x_i{v}_{y_i}+y_i{v}_{x_i})& 2{\sum }_i{m}_i{y}_i{v}_{y_i}& {\sum }_i{m}_i(y_i{v}_{z_i}+z_i{v}_{y_i})\\ {\sum }_i{m}_i(x_i{v}_{z_i}+z_i{v}_{x_i})& {\sum }_i{m}_i(y_i{v}_{z_i}+z_i{v}_{y_i})& 2{\sum }_i{m}_i{z}_i{v}_{z_i}\end{array}\right).\end{array}$$(39)

A safe, accurate method for calculating he derivative of the SVD, knowing İ is described in Appendix A. To do this, knowing the SVD decomposition of İ, the symmetric matrix H = U^TİU is first calculated. Its diagonal is the diagonal of $\dot{S}$, while its nondiagonal elements allow us to build an auxiliary antisymmetric matrix $\stackrel{\sim }{U}$-pagination. The time derivative of U is then found by $\dot{U}=U\stackrel{\sim }{U}$-pagination. The singular value derivative matrix $\dot{S}$ allows us to check the radial rigidity of the rotating pattern, which is small by hypothesis. Along the respective eigenvector i, each ${\dot{s}}_i$ provides a check of rigidity, and $s_i/{\dot{s}}_i$ a timescale for substantial radial change.

Then, using the unit vectors n_i and ṅ_i, i = 1, …3, of U = (n₁,n₂,n₃) and $\dot{U}=({\dot{\mathit{n}}}_1,{\dot{\mathit{n}}}_2,{\dot{\mathit{n}}}_3),$ respectively, corresponding to the singular values s_i, we derive the instantaneous rotation frequency vector of each of the principal axes along n_i,

$$\begin{array}{c}\hfill {\mathbf{\Omega }}_{I_i}={\mathit{n}}_i(t)\wedge \dot{{\mathit{n}}_i}(t),i=1,\dots 3.\end{array}$$(40)

If the particle ensemble rotates as a rigid figure, and if the axis of rotation is not aligned with any ṅ_i, the three Ω_Ii should be almost identical within the numerical errors. However, figures frequently rotate close to around one of the principal axes, in which case the corresponding |ṅ_j|, thus |Ω_Ij|, is small or zero. In the frequent case that the figure rotates around the shortest axis, Ω_I1 and Ω_I2 are almost identical and provide the sought pattern speed. Determining the instantaneous vector Ω_I of distinct subsets of particles allows us to quantify as a function of time by how much the rotation axes of these subsets may tilt with respect to each other.

6.1.2. Pattern speed of the 2D system

When the direction of the rotation axis is known and constant, for instance, along the z-axis, we can evaluate Ω explicitly avoiding SVDs. Noting

$$\begin{array}{cc}\hfill I_{\mathit{xx}}& ={\displaystyle \sum _i{m}_i{x}_i^2,I_{\mathit{yy}}=\sum _i{m}_i{y}_i^2,I_{\mathit{xy}}=\sum _i{m}_i{x}_i{y}_i,}\hfill \\ \hfill I_{\pm }& =\frac{1}{2}(I_{\mathit{xx}}\pm I_{\mathit{yy}})=\frac{1}{2}{\displaystyle \sum _i{m}_i(x_i^2\pm y_i^2),}\hfill \\ \hfill {\dot{I}}_{\mathit{xy}}& ={\displaystyle \sum _i{m}_i(x_i{v_y}_i+y_i{v_x}_i),{\dot{I}}_{-}=\sum _i{m}_i(x_i{v_x}_i-y_i{v_y}_i),}\hfill \\ \hfill P& =\sqrt{I_{-}^2+I_{\mathit{xy}}^2},\hfill \end{array}$$(41)

we calculate the eigenvalues λ_± by solving the characteristic polynomial of I,

$$\begin{array}{c}\hfill 0=|\mathit{I}-\lambda \mathbf{1}|={\lambda }^2-2I_{+}\lambda +I_{+}^2-P^2,\end{array}$$(42)

which gives two solutions for λ: λ_± = I₊ ± P. The eigenvalue ratio is a measure of the moment-of- inertia ellipsoid axis-ratio squared, thus the ellipticity reads

$$\begin{array}{c}\hfill ϵ=\sqrt{\frac{{\lambda }_{-}}{{\lambda }_{+}}}=\sqrt{\frac{I_{+}-P}{I_{+}+P}}.\end{array}$$(43)

The eigenvectors ξ_± only depend on the ratio ζ ≡ I_xy/I₋ ,

$$\begin{array}{c}\hfill {\mathit{\xi }}_{\pm }=\{\pm \zeta ,\sqrt{1+{\zeta }^2}\mp 1\}.\end{array}$$(44)

The two vectors can be verified to be orthogonal, that is, ξ₊ ⋅ ξ₋ = 0. Using identity 1.626.3 in Gradshteyn et al. (2007)³, the angle ϕ at time t between one of the eigenvectors and the x-axis in the x − y plane is

$$\begin{array}{c}\hfill \varphi (t)=\frac{1}{2}arctan\left(\zeta \right)+\mathrm{constant}.\end{array}$$(45)

The constant is 0, or ±π/2, which is irrelevant for our aim because it disappears when we differentiate ϕ with respect to time. We find

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathit{Iz}}\equiv \dot{\varphi }(t)=\frac{1}{2}\frac{\dot{\zeta }}{1+{\zeta }^2}=\frac{1}{2}\frac{I_{-}{\dot{I}}_{\mathit{xy}}-{\dot{I}}_{-}{I}_{\mathit{xy}}}{I_{-}^2+I_{\mathit{xy}}^2}.\end{array}$$(46)

The instantaneous angular speed of the moment of inertia at any time can therefore be calculated directly with the position and position-velocity moments given in Eq. (41).

When the particle distribution is close to axisymmetric, for instance, I_xy ≈ 0 and I₋ ≈ 0, Ω_z is undetermined, as it should. This also occurs when the distribution, while not axisymmetric, has strictly no bisymmetric component, however, such as a pure four−arm spiral pattern. This method was briefly described in Wu et al. (2018).

6.1.3. Pattern acceleration of the 2D system

The time derivative of Ω_Iz can be calculated similarly, knowing the particle accelerations, which are calculated in N-body simulations in any case. This can be useful for estimating the timescale over which the pattern speed can be approximated as constant. Differentiating Eq. (46), we obtain the relative pattern acceleration after simplifications,

$$\begin{array}{c}\hfill \frac{{\dot{\mathrm{\Omega }}}_{\mathit{Iz}}}{{\mathrm{\Omega }}_{\mathit{Iz}}}=\frac{I_{-}{\ddot{I}}_{\mathit{xy}}-I_{\mathit{xy}}{\ddot{I}}_{-}}{I_{-}{\dot{I}}_{\mathit{xy}}-I_{\mathit{xy}}{\dot{I}}_{-}}-2\frac{I_{-}{\dot{I}}_{-}+I_{\mathit{xy}}{\dot{I}}_{\mathit{xy}}}{I_{-}^2+I_{\mathit{xy}}^2},\end{array}$$(47)

where

$$\begin{array}{cc}& {\ddot{I}}_{\mathit{xy}}={\displaystyle \sum _i{m}_i(2{v_x}_i{v_y}_i+x_i{a_y}_i+y_i{a_x}_i),}\hfill \\ \hfill & {\ddot{I}}_{-}={\displaystyle \sum _i{m}_i({v_x}_i^2-{v_y}_i^2+x_i{a_x}_i-y_i{a_y}_i).}\hfill \end{array}$$(48)

Dimensionless estimators of the pattern speed invariance can be obtained by dividing Eq. (47) by a reference frequency, such as the solid-body speed Ω_s, or the pattern speed Ω_Iz itself. When the dimensionless estimator is of order unity or more, the constant pattern speed assumption is no longer valid.

6.2. Kinetic tensor rotation

The velocity distribution can rotate in velocity space with an angular speed that is distinct from the angular speed in position space. Exactly the same discussion as in Sect. 6.1.1 for the inertia tensor method can be carried out for the second-order moment tensor of velocities, the kinetic tensor (where the velocity origin is set at the velocity center of mass),

$$\begin{array}{ccc}\hfill \mathit{K}& \equiv \hfill & \hfill {\mathit{V}}^T\mathit{M}\mathit{V}=\left(\begin{array}{ccc}{\sum }_i{m}_i{v_x}_i^2& {\sum }_i{m}_i{v_x}_i{v_y}_i& {\sum }_i{m}_i{v_x}_i{v_z}_i\\ {\sum }_i{m}_i{v_x}_i{v_y}_i& {\sum }_i{m}_i{v_y}_i^2& {\sum }_i{m}_i{v_y}_i{v_z}_i\\ {\sum }_i{m}_i{v_x}_i{v_z}_i& {\sum }_i{m}_i{v_y}_i{v_z}_i& {\sum }_i{m}_i{v_z}_i^2\end{array}\right).\end{array}$$(49)

All the calculations for finding the 3D pattern speed with SVD decomposition can be made to obtain the pattern speed of the velocity ellipsoid by replacing positions by velocities, and velocities by accelerations, as in Sect. 6.1.1. We therefore do not elaborate further.

Likewise, for the 2D case, using these notations,

$$\begin{array}{cc}\hfill K_{\mathit{xx}}& ={\displaystyle \sum _i{m}_i{v_x}_i^2,K_{\mathit{yy}}=\sum _i{m}_i{v_y}_i^2,K_{\mathit{xy}}=\sum _i{m}_i{v}_x{v_y}_i,}\hfill \\ \hfill {\dot{K}}_{\mathit{xy}}& ={\displaystyle \sum _i{m}_i({v_x}_i{a_y}_i+{v_y}_i{a_x}_i),}\hfill \\ \hfill E_{\mathit{xy}}& =\frac{1}{2}(K_{\mathit{xx}}-K_{\mathit{yy}})=\frac{1}{2}{\displaystyle \sum _i{m}_i({v_x}_i^2-{v_y}_i^2),}\hfill \\ \hfill {\dot{E}}_{\mathit{xy}}& ={\displaystyle \sum _i{m}_i({v_x}_i{a_x}_i-{v_y}_i{a_y}_i),}\hfill \end{array}$$(50)

we obtain the pattern speed of the velocity ellipsoid,

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathit{Kz}}=\frac{1}{2}\frac{E_{\mathit{xy}}{\dot{K}}_{\mathit{xy}}-{\dot{E}}_{\mathit{xy}}{K}_{\mathit{xy}}}{E_{\mathit{xy}}^2+K_{\mathit{xy}}^2}.\end{array}$$(51)

To evaluate the pattern acceleration, we would need to know the particle jerk (ȧ).

7. Fourier method

A 2D mass distribution can have a simple angular Fourier decomposition where a low-order mode may capture the nonaxisymmetric shape of the distribution. This is commonly used to measure the strength, phase, and length of a bar, for instance, where the dominant m = 2 mode is often exclusively considered, while the m = 4, 6 and 8 modes are still non-negligible. The m = 2 and higher-order modes are thought to rotate at the same angular velocity as the mass distribution. In this case, we can attempt to determine the pattern speed by measuring just the mode phase velocity.

Below we describe the Fourier method introduced in Wu et al. (2018) more extensively. An inverse Fourier decomposition in the plane x − y of the mass-weighted azimuthal positions of a distribution of particles sampled in a ring reads

$$\begin{array}{c}\hfill F_m={\displaystyle \sum _j{m}_{j}exp(\mathrm{i}m{\theta }_j)=\sum _j{m}_j(cos(m{\theta }_j)+\mathrm{i}sin(m{\theta }_j)),}\end{array}$$(52)

where the azimuth angle θ_j of particle j depends on the Cartesian coordinates x_j, and y_j, θ_j = arctan(y_j, x_j), and m is the Fourier mode. The summation runs over the number of particles in the ring. The phase of the Fourier transform is

$$\begin{array}{c}\hfill {\varphi }_m=arctan(\mathfrak{I}(F_m),\mathfrak{R}(F_m)),\end{array}$$(53)

where

$$\begin{array}{cc}\hfill \mathfrak{R}(F_m)={\displaystyle \sum _j}{m}_{j}cos(m{\theta }_j),& \mathfrak{I}(F_m)={\displaystyle \sum _j{m}_{j}sin(m{\theta }_j).}\hfill \end{array}$$(54)

7.1. 2D Fourier mode velocity

Because the mode phase goes from 0 to 2π when the particle positions rotate in space by 2π/m, the mode velocity in space is the phase velocity divided by m,

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_m\equiv \frac{{\dot{\varphi }}_m}{m}=\frac{1}{m}\frac{\dot{\mathfrak{I}}(F_m)\mathfrak{R}(F_m)-\dot{\mathfrak{R}}(F_m)\mathfrak{I}(F_m)}{\mathfrak{R}{(F_m)}^2+\mathfrak{I}{(F_m)}^2},\end{array}$$(55)

where

$$\begin{array}{cc}\hfill \frac{\dot{\mathfrak{R}}(F_m)}{m}=-{\displaystyle \sum _j}{m}_{j}sin(m{\theta }_j){\dot{\theta }}_j,& \frac{\dot{\mathfrak{I}}(F_m)}{m}={\displaystyle \sum _j{m}_{j}cos(m{\theta }_j){\dot{\theta }}_j,}\hfill \end{array}$$(56)

and

$$\begin{array}{c}\hfill {\dot{\theta }}_j=\frac{x_j{\dot{y}}_j-y_j{\dot{x}}_j}{x_j^2+y_j^2}=\frac{x_j{v_y}_j-y_j{v_x}_j}{x_j^2+y_j^2}.\end{array}$$(57)

Explicitly, defining the following trigonometric sums,

$$\begin{array}{c}\hfill C\equiv {\displaystyle \sum _j{m}_{j}cos(m{\theta }_j),S\equiv \sum _j{m}_{j}sin(m{\theta }_j),C_1\equiv \sum _j{m}_{j}cos(m{\theta }_j){\dot{\theta }}_j,S_1\equiv \sum _j{m}_{j}sin(m{\theta }_j){\dot{\theta }}_j,}\end{array}$$(58)

we obtain the compact expression,

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_m=\frac{C{C}_1+S{S}_1}{C^2+S^2}.\end{array}$$(59)

When several modes superpose with different mode velocities, the mass distribution is time-dependent and no single pattern speed characterizes the distribution well. Determining how well the large amplitude (C² + S²)^1/2 mode velocities coincide is therefore a way to verify the time-independence of the pattern.

7.2. Fourier mode acceleration

Another way to confirm time dependence is to calculate the instantaneous mode acceleration. Differentiating Eq. (55), and defining

$$\begin{array}{c}\hfill C_2\equiv {\displaystyle \sum _j{m}_j[cos(m{\theta }_j){\ddot{\theta }}_j-msin(m{\theta }_j){\dot{\theta }}_j^2],S_2\equiv \sum _j{m}_j[sin(m{\theta }_j){\ddot{\theta }}_j+mcos(m{\theta }_j){\dot{\theta }}_j^2],}\end{array}$$(60)

with

$$\begin{array}{c}\hfill {\ddot{\theta }}_j=\frac{(x_j{a_y}_j-y_j{a_x}_j)-2{\dot{\theta }}_j(x_j{v_x}_j+y_j{v_y}_j)}{x_j^2+y_j^2},\end{array}$$(61)

we obtain the relative mode acceleration,

$$\begin{array}{c}\hfill \frac{{\dot{\mathrm{\Omega }}}_m}{{\mathrm{\Omega }}_m}=\frac{C{C}_2+S{S}_2}{C{C}_1+S{S}_1}+2m\frac{C{S}_1-S{C}_1}{C^2+S^2}.\end{array}$$(62)

Thus we can also calculate the instantaneous variation of the mode phase velocity from the particle coordinates, velocities, and accelerations at any time.

7.3. 3D extension for the m = 2 mode

Extending a Fourier-like decomposition for 3D distributions in order to determine a given mode speed vector is not straightforward. For example, choosing a decomposition in spherical harmonics or a double Fourier expansion over the sphere predetermines the positions of the poles and therefore breaks isotropy.

For the m = 2 mode in 2D, however, the Fourier method is equivalent to the moment of inertia method, as shown by substitution, where the positions x_i are replaced by unit vectors n_i = x_i/∥x_i∥. We showed above that the moment of inertia method works for 3D naturally using a SVD decomposition and differentiation. We can therefore use this method to replace the position vectors by unit vectors weighted by the mass.

The angular moment of inertia tensor I_n reads

$$\begin{array}{c}\hfill {\mathit{I}}_n\equiv {\mathit{n}}^T\mathit{M}\mathit{n}=\left(\begin{array}{ccc}{\sum }_i{m}_i{n}_{x_i}^2& {\sum }_i{m}_i{n}_{x_i}{n}_{y_i}& {\sum }_i{m}_i{n}_{x_i}{n}_{z_i}\\ {\sum }_i{m}_i{n}_{x_i}{n}_{y_i}& {\sum }_i{m}_i{n}_{y_i}^2& {\sum }_i{m}_i{n}_{y_i}{n}_{z_i}\\ {\sum }_i{m}_i{n}_{x_i}{n}_{z_i}& {\sum }_i{m}_i{n}_{y_i}{n}_{z_i}& {\sum }_i{m}_i{n}_{z_i}^2\end{array}\right),\end{array}$$(63)