# Dynamical modelling of galaxies with counter-rotating components using the Schwarzschild method

I am currently trying to understand dynamical modelling using Schwarzschild's method. I have been doing some reading the past days, but I suspect that I am getting something really wrong. According to my current understanding, the Schwarzschild method for dynamical modelling using spatially-resolved kinematic data works like this:

• Fit one or more 2D profiles (e.g., Gaussian) to the image or zeroth moment of the observation.
• De-project the above fitted model to 3D space and create a 3D luminosity distribution.
• Convert the 3D luminosity distribution to mass using a bunch of mass particles.
• Create a series of possible orbits for the above mass particles and calculate the resulting velocity vectors.
• From the above velocity vectors create velocity and velocity dispersion maps and calculate the chi-squared between them and the observed velocity and velocity dispersion maps.

The last step is what troubles me. To my understanding, the orbits of the particles in a Schwarzschild model can be counter-rotating. Right? In a galaxy observation, when we have two collisionless components (e.g., stellar components) that counter-rotate, it is possible that we will observe double-peaked spectral lines. Right? This information is impossible to be stored using a velocity and a velocity dispersion map, but it can be stored in a spectral cube. Yet, all papers on dynamical modelling using the Schwarzschild method fit models to velocity and velocity dispersion maps instead of a spectral cube.

So, from the above I have concluded that my understanding of the Schwarzschild method is not correct. I would really appreciate if someone could help me to understand what I am getting wrong.

Convert the 3D luminosity distribution to mass using a bunch of mass particles.

Well, you probably could do that, but I suspect it would be inaccurate and extremely time-consuming. Most Schwarzschild modeling approaches parameterize the 3D luminosity-density distribution using functions from which the gravitational potential can be relatively easily calculated.

For example, one popular approach is to approximate the observed 2D light distribution using the sum of multiple Gaussians (as you mention), which can be deprojected under certain assumptions into the sum of multiple 3D Gaussians to get the 3D luminosity-density distribution. Then you multiply this distribution by a mass-to-light ratio (usually one of the parameters of the model you're ultimately trying to fit to the data) and derive the gravitational potential from that. No need to approximate things with a bunch of particles.

This also allows you to easily add extra hypothetical "dark" components, like dark-matter halos or central supermassive black holes, which have simple analytical descriptions and don't need to modeled with particles.

To my understanding, the orbits of the particles in a Schwarzschild model can be counter-rotating.

Correct. In practice, if the system is (assumed to be) axisymmetric, you can calculate a set of orbits, and then easily generate their counter-rotating counterparts by just reversing the velocity vectors.

Yet, all papers on dynamical modelling using the Schwarzschild method fit models to velocity and velocity dispersion maps instead of a spectral cube.

You're clearly not reading the right papers ;-). Let me divide this into two parts:

## 1. Converting spectroscopic data to stellar-kinematic data

Before any modeling takes place, the observed spectroscopic data (long-slit and/or integral-field-unit "spectral cubes") must be analyzed to extract the stellar-kinematic data. This is done by dividing up the data into different spatial bins on the sky (bins along the slit, or 2D bins for a spectral cube) and summing the spectra within each bin, to improve the S/N.

Then you fit each spectrum with a weighted combination of stellar spectra (observed or from models), broadened by a model1 line-of-sight velocity distribution (LOSVD). In very simple cases, the LOSVD is assumed to be just a Gaussian with a specific central velocity $$V$$ and velocity dispersion $$\sigma$$. The more common present-day approach is to model the LOSVD with a Gaussian which is multiplied by third- and fourth-order Gauss-Hermite polynomials. The strength of the 3rd-order polynomial is parameterized by $$h_{3}$$ and measures skew; the 4th-order polynomial is parameterized by $$h_{4}$$ and measures kurtosis. For galaxies with counter-rotating components, $$h_{3}$$ and $$h_{4}$$ provide some information about the non-single-Gaussian nature of the LOSVD, although it's rather fuzzy and poorly constrained.

A more sophisticated approach (requiring high-S/N data) is to "non-parametrically" model the LOSVD with a histogram having twenty, thirty, or more velocity bins. To show how the non-parametric approach works, consider the following stellar-kinematic data for the S0 galaxy NGC 4191:

The first plot, from Fig. B1 of Krajnovic et al. 2011 (MNRAS, 414, 2923), shows $$V$$ and $$\sigma$$ maps (from the SAURON IFU). In the velocity field, you can see that the velocity reverses between $$r \sim 5$$ and 10 arc seconds; at roughly the same radius, $$\sigma$$ peaks. The paired off-center peaks in $$\sigma$$ are suggestive of a region where the two counter-rotating components are about equal (as argued by Krajnovic et al. in their paper). The second plot, based on Fig. 2 of Coccato et al. 2015 (A&A, 581, A65), shows the $$\sigma$$ map (colored panel, rotated 90 degrees with respect to the figure above) and non-parametric LOSVDs (31 velocity bins) for three of the bins (histograms) from VIRUS-W IFU observations of NGC 4191. You can see how, in the second histogram, the double peaks due to the overlapping, counter-rotating stellar components are clearly represented. ## 2. Comparing the Schwarzschild model with the kinematic data

After computing the library of stellar orbits for a given galaxy potential and assigning a set of weights to the individual orbits, you "observe" the model by orienting the galaxy to match the observed orientation, dividing it up into projected 2D spatial bins that match the 2D binning used for determining real galaxy's stellar kinematics. Then (for each bin) you compute the weighted LOSVD by integrating along the line of sight through that bin, adding up the radial-velocity contributions from each orbit that intersects the line of sight. If your data consists of $$V$$, $$\sigma$$, $$h_3$$, and $$h_4$$, then you fit the resulting model LOSVD for that bin with the Gauss-Hermite function to get the model bin's $$V$$, $$\sigma$$, $$h_3$$, and $$h_4$$. If your data are non-parametric LOSVDs, then you can simply bin the model LOSVD using the same velocity-binning scheme you used for the data.

This set of predicted LOSVD values (either $$V$$, $$\sigma$$, $$h_3$$, and $$h_4$$ or the non-parametric binned velocity histogram) is what you compare with data's LOSVD values for that spatial bin, repeating for all the spatial bins with kinematic data.

1 "Model" in this case having nothing to do with the Schwarzschild modeling.

• If you want to quote e.g. parts of the question I'd suggest using the > symbol (at the start of the line with a space after it) which the system will recognize as a quote for an entire paragraph. It also changes the displayed background color for that text block so it's easier to see. – StephenG Jun 16 '19 at 21:17
• @Peter, thank you for the response! However, I am confused with your final two paragraphs. You mean that the model outputs a cube with velocity values in spatial 3D space? How can you compare that with the data, since the usual data products are moment maps or spectral cubes? Maybe by converting that 3D spatial cube of velocities to a spectral cube of intensities by convolving an analytic profile (eg, negative Gaussian) with a stellar population model for each velocity bin? The mean of the analytic profile will be the line-of-sight velocity from the 3D spatial cube. Does this make sense? – AstrOne Jun 16 '19 at 22:50
• @AstrOne -- Sorry, I was a bit rushed and probably didn't explain it very clearly; I'll try to expand my answer later today. In the meantime: the spectroscopic data (long-slit or spectral cube) is converted into 2D maps of line-of-sight velocity measurements: e.g., V, $\sigma$, $h_{3}$, and $h_{4}$ in the "parameteric" case. Then the model is "observed" to produce similar maps, which are compared with the data-derived maps. – Peter Erwin Jun 17 '19 at 7:07
• @PeterErwin I really appreciate your help on this! To help you with your revised response, I will comment on your last comment. To my understanding, the spectroscopic data (eg, spectral cube) cannot be converted into 2D maps of line-of-sight velocity measurements if the observed spectral line profiles have two distinct peaks due to counter-rotating components. The only solution I can think of is to somehow perform decomposition and generate a set of V, σ, h3, and h4 maps for each peak/component. Then, run Schwarzschild's modelling multiple times, one for each set of maps. Does this make sense? – AstrOne Jun 17 '19 at 8:20
• @AstrOne -- You certainly can convert a datacube with multiply-peaked LOSVDs into 2D maps, it's just that the values in your maps may be very misleading. (In simple terms: you can fit a double-peaked profile with a single Gaussian, getting central value and dispersion -- it just won't be a good fit, and you'll get a misleading picture.) – Peter Erwin Jun 17 '19 at 9:45