<blockquote>
  <p>Convert the 3D luminosity distribution to mass using a bunch of mass particles.</p>
</blockquote>

<p>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.</p>

<p>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.</p>

<p>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.</p>

<blockquote>
  <p>To my understanding, the orbits of the particles in a Schwarzschild model can be counter-rotating.</p>
</blockquote>

<p>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.</p>

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

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

<h2>1. Converting spectroscopic data to stellar-kinematic data</h2>

<p>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. </p>

<p>Then you fit each spectrum with a weighted combination of stellar spectra (observed or from models), broadened by a model<a href="https://i.sstatic.net/lGeit.png" rel="nofollow noreferrer">1</a> line-of-sight velocity distribution (LOSVD). In very simple cases, the LOSVD is assumed to be just a Gaussian with a specific central velocity <span class="math-container">$V$</span> and velocity dispersion <span class="math-container">$\sigma$</span>. 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 <span class="math-container">$h_{3}$</span> and measures skew; the 4th-order polynomial is parameterized by <span class="math-container">$h_{4}$</span> and measures kurtosis. For galaxies with counter-rotating components, <span class="math-container">$h_{3}$</span> and <span class="math-container">$h_{4}$</span> provide <em>some</em> information about the non-single-Gaussian nature of the LOSVD, although it's rather fuzzy and poorly constrained.</p>

<p>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:</p>

<p>The first plot, from Fig. B1 of <a href="http://cdsads.u-strasbg.fr/abs/2011MNRAS.414.2923K" rel="nofollow noreferrer">Krajnovic et al. 2011 (<em>MNRAS</em>, 414, 2923)</a>, shows <span class="math-container">$V$</span> and <span class="math-container">$\sigma$</span> maps (from the SAURON IFU). In the velocity field, you can see that the velocity reverses between <span class="math-container">$r \sim 5$</span> and 10 arc seconds; at roughly the same radius, <span class="math-container">$\sigma$</span> peaks. The paired off-center peaks in <span class="math-container">$\sigma$</span> are suggestive of a region where the two counter-rotating components are about equal (as argued by Krajnovic et al. in their paper).</p>

<p><a href="https://i.sstatic.net/lGeit.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/lGeit.png" alt="Based on Fig.B1 of Krajnovic et al. (2011)"></a></p>

<p>The second plot, based on Fig. 2 of <a href="http://cdsads.u-strasbg.fr/abs/2015A%26A...581A..65C" rel="nofollow noreferrer">Coccato et al. 2015 (<em>A&amp;A</em>, 581, A65)</a>, shows the <span class="math-container">$\sigma$</span> 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.</p>

<p><a href="https://i.sstatic.net/Lmbew.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Lmbew.png" alt="Based on Fig.2 of Coccato et al. (2015)"></a></p>

<h2>2. Comparing the Schwarzschild model with the kinematic data</h2>

<p>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 <span class="math-container">$V$</span>, <span class="math-container">$\sigma$</span>, <span class="math-container">$h_3$</span>, and <span class="math-container">$h_4$</span>, then you fit the resulting model LOSVD for that bin with the Gauss-Hermite function to get the model bin's <span class="math-container">$V$</span>, <span class="math-container">$\sigma$</span>, <span class="math-container">$h_3$</span>, and <span class="math-container">$h_4$</span>. 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.</p>

<p>This set of predicted LOSVD values (either <span class="math-container">$V$</span>, <span class="math-container">$\sigma$</span>, <span class="math-container">$h_3$</span>, and <span class="math-container">$h_4$</span> 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.</p>

<p><a href="https://i.sstatic.net/lGeit.png" rel="nofollow noreferrer">1</a> "Model" in this case having nothing to do with the Schwarzschild modeling.</p>

CC BY-SA 4.0
2019-06-16T21:17:10.510Z
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 - Help Ukraine
CC BY-SA 4.0
2019-06-16T22:50:27.257Z
@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
CC BY-SA 4.0
2019-06-17T07:07:03.120Z
@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
CC BY-SA 4.0
2019-06-17T08:20:15.490Z
@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
CC BY-SA 4.0
2019-06-17T09:45:44.133Z
@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
CC BY-SA 4.0
2019-06-24T13:37:59.623Z
@PeterErwin thank you so much for the explanation, and please accept my apologies for the late response. Last week life struck me really hard and I didn't get the time to respond. :) I think I understood your explanation completely. I guess most of my confusion was less about Schwarzschild's method, and more about stellar kinematic modelling in general. So, essentially you deal with all the template-related stuff when you prepare the data for the fit. Then, your data and model (in the non-parametric approach) are just a bunch of (normalised) LOSVDs.
AstrOne
CC BY-SA 4.0
2019-06-24T13:40:07.653Z
(1/2) I think there is only only more thing left to wrap my head around before I completely understand what is going on. However, since my original question did not ask about it, I will accept your excellent answer regardless. Having that said, it would be highly appreciated if you could give me any pointers on this final question: I am not quite sure how such a model can account for beam smearing *correctly*.
AstrOne
CC BY-SA 4.0
2019-06-24T13:40:16.577Z
(2/2) A convolution of the LOSVD cube with a PSF won't be right, for the same reason the convolution of a velocity map with the PSF is not correct (unless the underlying surface brightness is completely flat). The convolution has to be between the surface brightness and the PSF. I wonder if this can easily be fixed by weighting/multiplying the modeled LOSVD cube with the modeled surface brightness image used to generate the 3D luminosity distribution at the beginning of the algorithm...
AstrOne
CC BY-SA 4.0
2019-06-16T14:56:06.727Z
2019-06-17T10:51:34.743Z
2019-06-17T10:51:34.743Z
Peter Erwin
<p>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:</p>

<ul>
<li>Fit one or more 2D profiles (e.g., Gaussian) to the image or zeroth moment of the observation.</li>
<li>De-project the above fitted model to 3D space and create a 3D luminosity distribution.</li>
<li>Convert the 3D luminosity distribution to mass using a bunch of mass particles.</li>
<li>Create a series of possible orbits for the above mass particles and calculate the resulting velocity vectors.</li>
<li>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. </li>
</ul>

<p>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.</p>

<p>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. </p>

<p>Thank you in advance.</p>

<galactic-dynamics><stellar-dynamics><modeling>
Dynamical modelling of galaxies with counter-rotating components using the Schwarzschild method
AstrOne
# Dynamical modelling of galaxies with counter-rotating components using the Schwarzschild method

(Asked by: username_0 on 2019-06-16)

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.

Thank you in advance.

**Tags:** galactic-dynamics, stellar-dynamics, modeling


---

## 1 answers

### Answer by username_1 on 2019-06-16. Score: 3 | Accepted answer

> 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 model[1](https://i.sstatic.net/lGeit.png) 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)](http://cdsads.u-strasbg.fr/abs/2011MNRAS.414.2923K), 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).

[![Based on Fig.B1 of Krajnovic et al. (2011)](https://i.sstatic.net/lGeit.png)](https://i.sstatic.net/lGeit.png)

The second plot, based on Fig. 2 of [Coccato et al. 2015 (*A&A*, 581, A65)](http://cdsads.u-strasbg.fr/abs/2015A%26A...581A..65C), 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.

[![Based on Fig.2 of Coccato et al. (2015)](https://i.sstatic.net/Lmbew.png)](https://i.sstatic.net/Lmbew.png)

## 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](https://i.sstatic.net/lGeit.png) "Model" in this case having nothing to do with the Schwarzschild modeling.

> **Comments:**
> **username_2**: 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.  
> **username_0**: @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?  
> **username_1**: @username_0 -- 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.  
> **username_0**: @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?  
> **username_1**: @username_0 -- 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.)  
> **username_0**: @PeterErwin thank you so much for the explanation, and please accept my apologies for the late response. Last week life struck me really hard and I didn't get the time to respond. :) I think I understood your explanation completely. I guess most of my confusion was less about Schwarzschild's method, and more about stellar kinematic modelling in general. So, essentially you deal with all the template-related stuff when you prepare the data for the fit. Then, your data and model (in the non-parametric approach) are just a bunch of (normalised) LOSVDs.  
> **username_0**: (1/2) I think there is only only more thing left to wrap my head around before I completely understand what is going on. However, since my original question did not ask about it, I will accept your excellent answer regardless. Having that said, it would be highly appreciated if you could give me any pointers on this final question: I am not quite sure how such a model can account for beam smearing *correctly*.  
> **username_0**: (2/2) A convolution of the LOSVD cube with a PSF won't be right, for the same reason the convolution of a velocity map with the PSF is not correct (unless the underlying surface brightness is completely flat). The convolution has to be between the surface brightness and the PSF. I wonder if this can easily be fixed by weighting/multiplying the modeled LOSVD cube with the modeled surface brightness image used to generate the 3D luminosity distribution at the beginning of the algorithm...  

---

