Import standard modules:
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import HTML
HTML('../style/course.css') #apply general CSS
Import section specific modules:
In [ ]:
pass
In [ ]:
HTML('../style/code_toggle.html')
The HI disks typically observed in spiral galaxies often show complex features that are difficult to parameterise and identify by simply studying different representations of the data by eye. However, obtaining such parameters is crucial in understanding how galaxies have formed and evolved. In order to extract the more complex information about these HI disks methods have been developed to quantify and extract the many parameters that describe the velocity structure and distribution of the neutral hydrogen in these disks. Initially this was predominantly done in order to extract reliable rotation curves of spiral galaxies but these days such methods are used to investigate many other open questions.
The foremost method of this higher order analysis is the tilted ring model. Originally introduced by Rogstadt et al. (1974)⤴ this method is now used in many different forms and shapes to facilitate analysis of HI data but, with the advancement of IFUs, also optical data that contains kinematical information. Its basic premise is that the disk of a galaxy can be represented by a series of concentric rings that are each described by their own kinematical and orientational parameters (See Fig 10.1).
In its simplest form each ring is described by a Position Angle (PA), Inclination (INCL) and a rotational velocity (VROT). It can be fitted to a velocity field, 2D fitting, or directly to the the datacube, 3D fitting, and both methods have its advantages and disadvantages. This chapter will mostly focus on the basics of tilted ring modelling of HI observations in both 2D and 3D. Along the way we will also brush on applications of the tilted ring model in IFUs as well as other advanced methods of analysis.
The most common way to fit a tilted ring model to the HI observations is to fit it to the velocity field (see Chapter 9). Once the assumption is made that the gas in the disk can be represented by the tilted ring model it becomes straightforward to fit the model to the velocity field. This is because the velocities of a rotating galaxy disk or ring in an image can easily be described with the simple formula: \begin{equation} v(x,y)=v_{sys} + v_{rot} \times sin(i) \times cos(\theta) \end{equation} With x and y the pixel coordinate, ${\rm v_{sys}}$ the systemic velocity of the galaxy, ${\rm v_{rot}}$ the rotational velocity, i the inclination and $\theta$ describing the angle between the major axis and the pixel with coordinate x,y.
This is related to the PA and image coordinates through: \begin{equation} cos(\theta) =\frac{-(x-XPOS) \times sin(PA) + (y-YPOS)\times cos(PA)}{r} \end{equation} With XPOS, YPOS the central coordinates of the rings and r their radius. This analytical representation can then be fitted to the velocity field in order to retrieve the parameters XPOS, YPOS, ${\rm v_{rot}}$, ${\rm v_{sys}}$, inclination and PA as function of radius.
The advantages of fitting the tilted ring model through this analytical representation are that simple fitting techniques can be used in order to retrieve the parameters. This also means that the fitting is fast and stable. Additionally the parameters that are being fitted can be expanded by adding other velocity components, such as radial velocities or streaming motions, to the analytical representation.
However, as the fit is performed on the velocity field it does mean that all information on the surface density distribution is lost. If one would require this distribution it would have to be retrieved separately from the integrated moment map. This could be done by performing an ellipse fit on the integrated moment map, possible with the orientation of the ellipses defined by the fit to the velocity field.
Additionally, the quality of the fit will directly depend on the quality of the velocity field. This means that any errors present in the velocity field will be propagated into the tilted ring fit. Some of these errors are well known. For example, as the line profile in each pixel is represented by a single velocity the more complex kinematical structures are not accurately represented. An example of such a complex structure is one where the disk crosses the sightline multiple times. In a case where the disk crosses the sightline twice, for example due to a warp, the line profile will show peaks intensities at two different velocities. However, in the velocity field these will be represented by a single value, thus losing vital information that will complicate the fitting. The clearest case of this is an edge-on galaxy where the disk crosses each sight line an infinite number of times. As such, the tilted ring model should not be fitted to such galaxies using the velocity field. There is no hard upper inclination limit for the usage of the velocity field to fit tilted ring models, however it is generally accepted that for i>75 degrees the projection effects become too severe to still fit the models reliably.
There is also a lower inclination limit at which the tilted ring model can no longer be fitted accurately to velocity field. That such a limit must exist is clear when thinking about a disk that is exactly face-on. As the observations only measure velocity along the line of sight and all the rotational movement in the disk is exactly in the plane of the sky an observation of a face-on galaxy contains no information about the rotation curve.
In reality very few galaxies will be exactly face-on, however the example shows that there must be a point where the rotation curve can no longer be extracted from the velocity field reliably. Begeman (1989)⤴ has shown that for i < 40 degrees the chi square landscape becomes such that the minimum is ill defined in velocity. Hence it is generally accepted that tilted ring modelling through fitting the velocity field should only be applied to galaxies in the inclination range 40-75 degrees. Note that these limits are set due to the method of fitting the analytical formula to the velocity field and not, as sometimes thought, due to the intrinsic idea of the tilted ring model.
Velocity field are also hampered by instrumental effects. The most dominant of these are resolution effects. If the observed galaxy is not well-resolved the low resolution will lead to the line profiles having a long tail towards the systemic velocity. If not corrected for this tail will lead to lower apparent velocities. This effect is called “beam-smearing” and will bias the fitted velocities to lower values than present in the galaxy. Bosma (1978)⤴ already showed early on that the slope of the retrieved rotation curve will be biased low when less than 7 resolution elements are present within a Holmberg radius.
Due to the many parameters involved velocity field fitting, the fit is typically underconstrained and hence the fitting process has to be guided interactively. Generally a fit would start with a set of initial estimates for all parameters and in the first fit one would fix the most reliable initial estimates (usually the central coordinates) while fitting the least reliable. Through a series of fits in which different parameters are allowed to vary each time one then finally arrives at the fitted model. The order in which parameters are allowed to vary and fixed often differs between researchers and also depends on the program that is used. An exact example of this procedure for the most commonly used program, ROTCUR in GIPSY (Begeman 1987⤴; van der Hulst et al. 1992⤴), can be found in de Blok et al. (2008)⤴
From equation 10.1 it is easy to see that around the minor axis ($\theta \sim 0$) the velocity field carries very little information about the rotation curve. Therefore, usually one would exclude a given angle around the minor axis from being included in the fit. The angle is excluded is usually in the range 10$^{\circ}$-30$^{\circ}$.
Due to the iterative process the parameterisation of the HI observations is still a manual process. This means that parameterised samples are still limited to tens of galaxies. With the imminent arrival of large high resolution surveys such as WALLABY, WNSHS and MALS an automated procedure will be required. Efforts to do this will be described in section 10.4.
To overcome the limitations of 2D fitting, e.g. limited inclination range and beam smearing, in the last decade or so a push has been made to fit the tilted ring model directly to the data cube. This has the huge advantage that all information of the observation can be utilised in the fitting procedure leading to the possibility of fitting more complex HI distributions.
In 3D there is no simple analytical formula that can map the tilted ring model to the data cube. This means that in order to fit the model to the data a full mock observation, based on a set of parameters, is required. This mock observation can subsequently be subtracted from the data and the residuals can then be compared and minimised. The common approach is to populate the tilted ring model, as described by the input parameters, randomly with point sources of a given flux and subsequently convolve the model to the resolution of the data, making sure that instrumental effects are included. The fact that a full model needs to be created and convolved with a beam at every iteration creates a bottleneck in speed, and thus 3D methods are slower than 2D methods. On the other hand, because the actual minimiser compares a mock observation to the data there are no limits on the complexity of the model.
Additionally as the model is built in the frame of the galaxy with point sources, based on the circular ring premise of the tilted ring model, the fitting is much less affected by projection effects or beam smearing. This is easily seen once more taking the example of an edge-on galaxy. As effectively the full velocity structure of the line profile is modelled and compared in 3D fitting the long tails towards the systemic velocity caused by projection are automatically reproduced in the model. In order to populate the model with flux the surface brightness and the scale height of the disk need to be defined as well. This means that true 3D models require these parameters to be fitted as well. This makes the fitting process more complex but does allow for the extraction of the surface brightness profiles and scale heights from the data.
In it has been shown that edge-on galaxies can be reliably fitted in an iterative method (Kamphuis et al. 2011⤴; Zschaechner et al. 2011⤴) that is analogous to the 2D fitting method described above. Additionally, it allows for the fitting of disks that cross the line of sight multiple times (Józsa 2007⤴), multiple disks (Gentile et al. 2013⤴) and non-cylindrically symmetric models (Kamphuis et al. 2013⤴). These features have made 3D tilted ring modelling crucial in the search for extra-planar gas and determining its distribution and kinematics (Oosterloo et al. 2007⤴; Heald et al. 2011⤴; Gentile et al. 2013⤴; de Blok et al. 2014⤴).
Although it is now clear that 3D tilted ring modelling allows for the fitting of galaxies with i>75$^{\circ}$, the lower limit in inclination as well the minimum required size is not yet clear. Early work on automated fitting indicates that inclinations as low as 20$^{\circ}$ can be reliably fitted (Kamphuis et al. 2015⤴) as well as galaxies with only 2 resolution elements with the observed radius (di Teodoro & Fraternali 2015⤴).
Until recently the only publicly available code for 3D fitting has been the Tilted Ring Fitting Code (TiRiFiC,Józsa 2007⤴, http://gigjozsa.github.io/tirific/). Besides the previously discussed advantages of 3D fitting over 2D, TiRiFiC has greatly expanded the type and amount of parameters that can be fitted. As, due to its 3D nature, full models need to be produced as part of the fitting, the code produces a model data cube as standard output. This means that the fitted models can easily be compared to the data by comparing contour levels or displaying the cubes next to each other in kvis.
TiRiFiC initially performed a $\chi^2$ minimization through sequential fitting of all parameters using the golden-section method (Press et al. 1992⤴). Even though in the current version other fitting paradigms are available this is still the most commonly used minimisation technique. The sequential fitting of the rings leads to several effects. First of all, for the higher inclined galaxies it is important to fit the rings in an outside-in sequence as the outer rings are projected onto the inner rings (Kregel & van der Kruit 2004⤴). Additionally, small errors in one ring can be compensated in the following ring in the sequence. This can lead to an artificial sawtooth pattern in the parameters which needs to be corrected.
Due to the additional parameters and the sequential ring fitting it is still necessary to manually adjust and guide the TiRiFiC fitting. However, also for 3D efforts are underway to automate this process.
In the 2D fitting of the velocity field an analytical formula is fitted, and therefore the formal errors of the fit can be calculated. However, due to various sources of errors in the data, e.g. radial motions or errors in the construction of the velocity field, these formal errors underestimate the true uncertainties. Currently, the most commonly accepted way to obtain errors on the rotation curve is to divide the galaxy along the minor axis and to fit both sides independently. It is then assumed that the potential of the galaxy is spherical and hence that the rotation curve on both sides is the same. This means that the difference between the rotation curves retrieved when both sides are fitted independently and the rotation curve that is retrieved for the full galaxy, either by fitting the full velocity field or by averaging both sides, would then contain all possible errors except systematics. This still leaves the other parameters without errors though.
In 3D fitting the situation is even worse as no formal errors are available. Even though the residuals of data and model can be used for the $\chi^2$ minimisation it is difficult to translate these into errors for the individual parameters. This is because many of the parameters in the model are interdependent and therefore there is no straightforward way to translate the residuals to individual errors. In the literature the errors have been estimated by modifying the individual parameters in the best fit model one by one until they clearly deviate from the data (Gentile et al. 2013⤴; Kamphuis et al. 2013⤴; de Blok et al. 2014⤴). This method gives realistic errors even though it still ignores the interdependence of the parameters.
In order to obtain the errors in a formally correct way that takes the interdependence of the parameters as well as the various sources into account one would need to develop a bootstrapping or Monte Carlo method. The problem there is that for these methods to give reliable estimates one would need to create thousands of realisations of the model. Especially for 3D models this is currently far too computationally expensive.
Currently there are already several HI surveys that have observed tens to hundreds of galaxies. Surveys such as the Westerbork HI Survey of Irregular and Spiral galaxies (WHISP; van der Hulst 2002⤴">van der Hulst 2002</a>), the Local Volume HI Survey (LVHIS; Koribalski 2010⤴), DiskMass (Martinsson 2011⤴), ATLAS3D (Serra et al. 2012⤴) and the Blue disc Project (Wang et al. 2013⤴) test the limits of surveys for which models can still be obtained interactively.
In the near future surveys such as, The Widefield ASKAP L-band Legacy All-sky Blind surveY (WALLABY; Koribalski 2012⤴) will provide a widefield HI survey at ∼15–30 arcsec resolution, which will result in thousands of galaxies that are resolved to an extent that their kinematical parameters can be retrieved (Duffy et al. 2012⤴). This clearly shows the need for automated routines that can accurately parameterise the data.
Such routines will have the added effect that all models will be produced to the same accuracy and that any bias introduced in the interactive fitting is negated. 2015 has seen the release of several packages (di Teodoro & Fraternali 2015⤴; Kamphuis et al. 2015⤴; Bekiaris et al. 2016⤴) that aim to do automated fitting without the need for interaction. Notably all these codes are 3D efforts.
In order to be able to fit the tilted ring model without the need for interaction the codes need to determine the extent of the emission in the observations as well as to ensure that the problem is overconstrained. This means that these codes require built in source finders (see Chapter 8) and that the models are constrained to only a few parameters. In practice this means that these codes return the rotation curve and, in some cases, the PA, Inclination and surface brightness as function of radius. Additionally, they return the central coordinates, dispersion and possibly the scale height as single values. To add more parameters would currently seriously affect the stability of these codes.