<figure><img src="image_1.png" /></figure>

HAL
open science

# Constraining the dark matter content of NGC 1291 using hydrodynamic gas response simulations

F. Fragkoudi, E. Athanassoula, A. Bosma

► To cite this version:

F. Fragkoudi, E. Athanassoula, A. Bosma. Constraining the dark matter content of NGC 1291 using hydrodynamic gas response simulations. Monthly Notices of the Royal Astronomical Society, 2017, 466 (1), pp.474-488. 10.1093/mnras/stw3023. hal-01554697

HAL Id: hal-01554697

https://hal.science/hal-01554697v1

Submitted on 26 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Constraining the dark matter content of NGC 1291 using hydrodynamic gas response simulations

F. Fragkoudi,<fnref n="1" />,<fnref n="2" />* E. Athanassoula<fnref n="2" /> and A. Bosma<fnref n="2" />

Accepted 2016 November 17. Received 2016 November 16; in original form 2016 October 18

## ABSTRACT

We present a pilot study on the nearby massive barred galaxy NGC 1291, in which we use dynamical modelling to constrain the disc mass-to-light ratio (M/L), thus breaking the degeneracy between the baryonic and dark matter in its central regions. We use the gas, specifically the morphology of the dust lanes on the leading side of the bar, as a tracer of the underlying gravitational potential. We run a large number of hydrodynamic gas response simulations, in potentials obtained directly from near-infrared images of the galaxy, which have three free parameters: the M/L, the bar pattern speed and the height function. We explore the three-dimensional parameter space, by comparing the morphology of the shocks created in the gas response simulations with those of the observed dust lanes, and find the best-fitting models; these suggest that the M/L of NGC 1291 agrees with that predicted by stellar population synthesis models in the near-infrared (≈0.6 M⊙/L⊙), which leads to a borderline maximum disc for this galaxy. Furthermore, we find that the bar rotates fast, with a corotation radius that is ≤1.4 times the bar length.

**Key words:** galaxies: kinematics and dynamics – dark matter.

## 1 INTRODUCTION

Since the discovery in the 1970s that rotation curves derived from H1 observations – which extend well beyond the optical disc of galaxies – remain flat (Roberts 1976; Bosma 1978, 1981), the nature and distribution of non-baryonic dark matter in galaxies has been a source of debate in the scientific community. One of the main issues with determining the non-baryonic distribution of matter in the central regions of galaxies is what has become known as the ‘disc–halo degeneracy’. This degeneracy arises due to the fact that rotation curve decompositions depend critically on the mass-to-light ratio (M/L or Y) of the stellar disc. van Albada et al. (1985) showed that the rotation curve of a galaxy can be fit equally well with a variety of dark matter halo and stellar disc components, from models with barely any disc contribution, to models where the disc contributes as much as possible to the rotation curve – without creating a hole in the dark matter halo; this is also known as the ‘maximum disc hypothesis’. Sackett (1997) gave a more concrete definition of the maximum disc by stating that a disc is maximum if it contributes 85 ± 10 per cent of the rotation curve at $r = 2.2 h_r$, where $h_r$ is the scalelength of the disc. However, there is no conclusive evidence to either prove or disprove the maximum disc hypothesis and in fact there are a number of arguments both for (Athanassoula, Bosma & Papaioannou 1987; Sackett 1997; Weiner, Sellwood & Williams 2001) and

against it (Kuijken & Gilmore 1991; Courteau & Rix 1999), for the Milky Way (Bissantz & Gerhard 2002; Bissantz, Englmaier & Gerhard 2003; Hamadache et al. 2006; Tisserand et al. 2007) as well as for external galaxies (Bottema 1993; Trott & Webster 2002; Kranz, Slyz & Rix 2003; Gnedin et al. 2007). The lack of knowledge regarding the dark matter content in the central regions of galaxies is of course a major roadblock for testing theories of galaxy formation and evolution.

One way to obtain an estimate of the disc M/L, and therefore break the disc–halo degeneracy, is from stellar population synthesis (SPS) models. These however require a number of assumptions that depend on the detailed star formation history, the modelling of late phases of stellar evolution and the initial mass function (IMF), quantities that in general contain a large amount of uncertainty, and which propagate through to the determination of the physical parameters of galaxies (Bell & de Jong 2001; Conroy, Gunn & White 2009).

Another way to break the degeneracy between the disc and the dark matter halo is by using dynamical estimates of the M/L of the disc. In this study, we employ a method that uses the gas in galaxies as a tracer of the underlying gravitational potential, and which we refer to as the hydrodynamic gas response method. This method has already been used in studies in the past for both barred galaxies (Lindblad, Lindblad & Athanassoula 1996; Weiner et al. 2001; Pérez, Fux & Freeman 2004; Weiner 2004; Zánmar Sánchez et al. 2008) and spiral galaxies (Kranz, Slyz & Rix 2001, 2003).

<footnote n="1">¹GEPI, Observatoire de Paris, CNRS, Université Paris Diderot, 5 Place Jules Janssen, F-92195 Meudon, France</footnote>

<footnote n="2">²Aix Marseille Univ, CNRS, LAM, Laboratoire d'Astrophysique de Marseille, F-13013 Marseille, France</footnote>

<footnote>* E-mail: francesca.fragkoudi@obspm.fr</footnote>
This type of modelling has also been applied in order to determine the pattern speed of barred galaxies (Pérez et al. 2004; Sánchez-Menguiano et al. 2015). The method relies on simple principles: gas, being collisionless, will follow circular orbits when placed in an axisymmetric potential. However, when strong non-axisymmetries are present, the gas is forced to shock, thus creating dust lanes, as observed in many barred galaxies (Prendergast 1983; Athanas-soula 1992b). The strength and morphology of these shocks depends, among other things, on the ratio between the amount of baryonic to dark matter. Therefore, by tuning the *M/L* of the disc, such that the simulations reproduce the morphology and kinematics of the gas in the observations, we can also constrain the amount of dark matter.

In addition to its dependence on the *M/L* of the disc, the morphology and kinematics of the gas also depend on the bar pattern speed and on the height function of the disc. The effect of the latter has not been extensively studied in the past; however, we show in the following sections that it plays a crucial role on the strength and shape of the shocks created in the gas. Due to the dependence of the gas morphology on these three parameters, we need to run a large number of simulations, in order to explore the full three-dimensional parameter space, to find the model that best reproduces the observations.

This paper is structured as follows: in Section 2, we describe the observational data we have for NGC 1291, which we use to construct the models. In Section 3, we present the dynamical models and the gas response simulations. In Section 4, we describe the main results of the paper, the trends that arise from changing the three free parameters and the best-fitting models. In Section 5, we discuss the findings of this study in the context of previous results in the literature, and in Section 6, we give a summary of the main results.

## 2 OBSERVATIONS OF NGC 1291

NGC 1291 is an early-type barred galaxy (SB0/a; de Vaucouleurs et al. 1991) that, in addition to the primary bar, contains a smaller nuclear bar, as well as an outer ring, in which there is active star formation (Thilker et al. 2007). The inclination of the galaxy is at most 12° so it is practically face-on (Bosma et al. 2010) and it is at a distance of 9.7 Mpc (Tully et al. 2009), while the galaxy does not seem to be undergoing any interactions.

The galaxy contains hot gas primarily in the bulge region, as evidenced by X-ray studies (Bregman, Hogg & Roberts 1995; Hogg et al. 2001) as well as H I gas, which has a hole in the inner parts of the galaxy where the X-ray gas is found (Bosma et al. 2010). The fraction of H I gas mass compared to the stellar mass<fnref n="1" /> is 2.6 per cent (Bosma et al. 2010) and the mass of the hot gas found in the galaxy is even less (Bregman et al. 1995).

For the work presented here, we made use of the 3.6, 4.5, 8 and 24 µm images (see Fig. 1) of NGC 1291. The 3.6 and 4.5 µm images are taken from the Spitzer Survey of Stellar Structure in Galaxies (S⁴G; Sheth et al. 2010) and are used to obtain the gravitational potential of the galaxy, since these wavelengths are well suited for tracing the old stellar population, where most of the mass of the galaxy resides. The 3.6 µm image was also used to perform photometric decompositions – using the BUDDA code (Gadotti 2008) –

<figure><img src="image_4.png" /><figcaption>**Figure 1.** From left to right: the 3.6, 8 and 24 µm images for NGC 1291. In the top panels, the whole galaxy is shown, including the outer ring, while in the bottom panels, we zoom in on a region encompassing the bar. The red line outlines the bar as obtained from two-dimensional decompositions of the 3.6 µm image.</figcaption></figure>

<table><caption>**Table 1.** Some basic properties of NGC 1291.</caption><thead><tr><td>Property</td><td>Value</td></tr></thead><tbody><tr><td>Classification</td><td>SB0/a</td></tr><tr><td>Inclination (deg)</td><td>12</td></tr><tr><td>Distance</td><td>9.7 Mpc</td></tr><tr><td>1 arcsec</td><td>47 pc</td></tr><tr><td>Disc scalelength</td><td>5.8 kpc</td></tr><tr><td>Radius of primary bar</td><td>5.1 kpc</td></tr></tbody></table>

to obtain the scalelength of the disc and the bar length, which are $h_r = 5.8$ kpc and $r_B = 5.1$ kpc, respectively (Bosma et al. 2010).

To extract information about the dust lanes in the galaxy (see the middle and right panels of Fig. 1), we use the 8 and 24 µm images from the SINGS survey (Kennicutt et al. 2003), which trace the emission due to polyaromatic hydrocarbons (PAHs) and hot dust, respectively, and correlate with the presence of molecular gas (Regan et al. 2006; Bendo et al. 2007; Schinnerer et al. 2013). PAH emission is correlated to the emission from cold dust $\le 20$ K (Bendo et al. 2010) and is likely to occur in shocked regions of the interstellar medium (ISM; Tielens 2008).

Some basic properties of the galaxy are summarized in Table 1.

## 3 DYNAMICAL MODELS AND SIMULATIONS

To obtain a dynamical model of NGC 1291, we need to account for the mass distributions of the different components that make up the potential, such as the stellar (Section 3.1) and dark matter (Section 3.2) component. In general, one also needs to account for the gaseous component in the galaxy; however, as mentioned in Section 2, the H I and X-ray gas are negligible in the bar region of NGC 1291 and we therefore use the gas only as a tracer of the underlying potential.

### 3.1 Potential of the stellar component

The method for calculating the potential due to the stellar component involves a straightforward three-dimensional integration over the density distribution, which was described in more detail along with

<footnote n="1">¹ Where the stellar mass is obtained using the *M/L* from S4G (Meidt et al. 2014; Röck et al. 2015), which as we will show further down is a reasonable assumption.</footnote>
tests on the code in Fragkoudi et al. (2015). Once a mass map for the galaxy is constructed from the near infrared (NIR) image of the galaxy – by assigning an M/L – we then assign a height function and integrate over the density distribution to obtain the gravitational potential.

### 3.1.1 The mass-to-light ratio

To construct a mass map of the galaxy, we first remove contaminating sources, such as foreground stars, apply Fourier smoothing to the image, and then assign an M/L to convert the luminosity to mass. The M/L in the 3.6 µm band can be assumed to be 0.6 (Meidt et al. 2014; Röck et al. 2015), once the contamination by dust has been removed. To correct for dust contamination, we combine information from the 3.6 and 4.5 µm images, using a formula by Eskew, Zaritsky & Meidt (2012) with a new calibration by Querejeta et al. (2015), obtained using a large number of S$^4$G galaxies. The equation for transforming from the 3.6 and 4.5 µm flux densities ($F_{3.6}$ and $F_{4.5}$ in Jy) to mass (in units of M$_{⊙}$), as given by Querejeta et al. (2015) is

$$ \frac{M_{*}}{M_{\odot}} = 10^{8.35} \times \left( \frac{F_{3.6}}{Jy} \right)^{1.85} \times \left( \frac{F_{4.5}}{Jy} \right)^{-0.85} \times \left( \frac{D}{\text{Mpc}} \right)^{2}, \quad (1) $$

where the distance D is given in Mpc.

### 3.1.2 The height function

We assume a height function, which defines how the density drops off as a function of z from the equatorial plane. The height function and the scaleheight ($z_0$) smooth out the potential, and we therefore need to use the height function that best approximates that of the galaxy we are trying to model. Different functional forms can be used for the height function, and in this study, we explore two of the most commonly used ones in the literature, the isothermal-sheet model and an exponential height function (see the top panel, Fig. 2).

The isothermal-sheet model (van der Kruit & Searle 1981) is given by

$$ F(z) = \frac{1}{2z_0} \mathrm{sech}^2 \left( \frac{z}{z_0} \right), \quad (2) $$

where $1/(2z_0)$ is the normalization factor and $z_0$ is the scaleheight (also referred to as the isothermal scaleheight when a clearer distinction needs to be made from the exponential scaleheight).

Studies of edge-on galaxies have shown that the vertical profiles of galaxies can be more peaked than that expected for an isothermal height function and better matched by an exponential distribution (de Grijs, Peletier & van der Kruit 1997; de Grijs 1998). To check whether our assumptions on the form of the height function would have a significant effect on the results, we ran some additional suites of models with an exponential height function, which is given by

$$ F(z) = \frac{1}{2h_z} \exp \left( -\frac{z}{h_z} \right), \quad (3) $$

where $h_z$ is the exponential scaleheight that is related to the isothermal scaleheight by $z_0 = 2h_z$.

We also explore height functions that include the shape of the boxy/peanut (b/p) bulge, since, as was shown in Fragkoudi et al. (2015), when a b/p bulge is present in a galaxy it should be modelled, because it can have a significant effect on the model. We use a b/p

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Top panel: isothermal and exponential height functions with equivalent scaleheights ($z_0 = 2h_z = 0.8$ kpc). We see that the exponential height function is more peaked, which leads to higher mass concentrations close to the plane of the galaxy. Bottom panels: on the left we show the behaviour of $z_0$, the isothermal scaleheight for a disc with constant scaleheight, and on the right we show how $z_0$ varies for a disc with our b/p height function. The vertical dashed lines mark the end of the bar.</figcaption></figure>

height function that is a non-separable function of position, given by

$$ F(x, y, z) = \frac{1}{2z_0(x, y)} \mathrm{sech}^2 \left( \frac{z}{z_0(x, y)} \right), \quad (4) $$

where the scaleheight $z_0(x,y)$ varies like the sum of two two-dimensional Gaussians (see the bottom right panel of Fig. 2):

$$ z_0(x, y) = A_{\mathrm{sim}} \exp \left( - \left( \frac{(x - x_0)^2}{2\sigma^2} + \frac{(y - y_0)^2}{2\sigma^2} \right) \right) + A_{\mathrm{sim}} \exp \left( - \left( \frac{(x - x_1)^2}{2\sigma^2} + \frac{(y - y_1)^2}{2\sigma^2} \right) \right) + z_0^{\mathrm{disc}}, \quad (5) $$

where $A_{\mathrm{sim}}$ is the maximum scaleheight of the peanut above the disc scaleheight ($z_0^{\mathrm{disc}}$) and $\sigma^2$ is the variance of the Gaussians. Parameters ($x_0, y_0, x_1, y_1$) give the positions of the maxima of the Gaussians.

To obtain an initial set of values for the parameters of the b/p height function, we fitted the vertical mass distribution of an N-body simulation of an isolated galaxy, taken from Athanassoula, Machado & Rodionov (2013, for more details on the fitting, see Fragkoudi et al. 2015). In this study, we adopt a fiducial value for the maximum scaleheight of the peanut $A_{\mathrm{fid}}$, which is half that obtained by the fit to the simulated galaxy $A_{\mathrm{sim}}$, i.e.

$$ A_{\mathrm{fid}} = \frac{A_{\mathrm{sim}}}{2} = 0.65 \, \mathrm{kpc}. \quad (6) $$

This is done because the height of the b/p, also known as the strength of the b/p, is known to be correlated with the bar strength (Athanassoula 2008; Martínez-Valpuesta & Athanassoula 2008) – although there is a large amount of scatter in this relation.
<figure><img src="image_6.png" /><figcaption>Figure 3. Rotation curves for models 001–004. The M/L is decreasing from 1 × $\Upsilon_{3.6}$ to 0.25 × $\Upsilon_{3.6}$ in steps of 0.25.</figcaption></figure>

The simulated galaxy used to fit the b/p height function develops a strong bar, stronger than that of NGC 1291. We therefore reduce the strength of the b/p for the model of NGC 1291, since such a strong b/p is not likely to be present in more weakly barred galaxies.

## 3.2 Dark matter halo potential

The dark matter halo in our models is a spherical pseudo-isothermal halo whose density is given by

$$ \rho(r) = \rho_c \left[ 1 + \left( \frac{r}{r_c} \right)^2 \right]^{-1}, \qquad (7) $$

where $r_c$, the core radius, and $\rho_c$, the central density, are adjusted such that $V_{\text{flat}}$ is the same in all the models. The circular velocity at a given radius $r$ is given by

$$ V(r) = V_{\infty} \left[ 1 - \frac{r_c}{r} \arctan \left( \frac{r}{r_c} \right) \right]^{1/2}, \qquad (8) $$

where $V_{\infty}$ is the circular velocity at infinity. The parameters of the halo are obtained by adjusting the rotation curve of the model, such that it reproduces the rotation velocity in the outer parts ($V_{\text{flat}}$), as predicted by the Tully-Fisher (TF) relation for NGC 1291 (McGaugh et al. 2000; McGaugh 2012). For NGC 1291 $V_{\text{flat}} = 200$ km s$^{-1}$, and we use this value to construct the rotation curves (see Fig. 3).

## 3.3 The Lagrangian radius

The Lagrangian radius is the distance from the centre of the galaxy to the Lagrangian points $L_1$ (or $L_2$) that are the locations at which the forces in the rotating frame of reference are exactly zero, i.e. in the rotating frame of reference, the stars in the disc rotate at the same speed as the bar. The Lagrangian radius is often parametrized by

$$ r_L = R r_B, \qquad (9) $$

where $r_L$ is the Lagrangian radius, $r_B$ is the bar semimajor axis and $R$ is a parameter which is often used to separate bars into fast $R \le 1.4$ and slow $R > 1.4$ bars (e.g. Debattista & Sellwood 2000).

## 3.4 Hydrodynamic simulations

The gas dynamic (or gas response) simulations were run using the hydrodynamic adaptive mesh refinement grid code RAMSES (Teyssier 2002). The simulations are run in two dimensions, since, to zeroth order, gas is confined to a thin plane in disc galaxies, which allows us to therefore decrease the complexity of the problem.

The initial conditions of the simulations consist of a two-dimensional axisymmetric gaseous disc, with a homogeneous density distribution, in hydrostatic equilibrium. The disc is tapered using an exponential tapering function. The gas is modelled as an isothermal gas with adiabatic index 5/3 which corresponds to H$^1$ gas (neutral atomic hydrogen), and we give a value of 10 km s$^{-1}$ to the sound speed of the gas, which corresponds to characteristic temperatures of the ISM.

In order to avoid transients in the disc, we introduce the non-axisymmetric potential over a finite period of time, which was found empirically as ~3 bar rotations, in accordance with previous studies (Athanassoula 1992b; Patsis & Athanassoula 2000).

Sormani, Binney & Magorrian (2015) found that the gas flows in non-axisymmetric potentials are sensitive to the resolution of the grid used in these types of simulations. After carrying out a number of tests on analytic models containing a disc, bar and bulge, taken from Athanassoula (1992a), we determined that the location of the shocks remains the same (within 5 per cent) for varying resolutions, while the central region of the simulation is indeed sensitive to resolution. For more details on this, we refer the reader to the appendix. For the purposes of this study, our simulations have a box size of 60 kpc, and a maximum resolution in the shock regions of 30 pc.

## 3.5 Varying the free parameters

The goal of this work is to find the best fit between the gas response simulations and the observations, by exploring the allowed parameter space, by varying the three free parameters of the models, i.e. the disc M/L, the Lagrangian radius and the height function of the disc.

In total, we ran approximately 300 models, where we varied the M/L from 1.5 to 0.25 $\Upsilon_{3.6}$ in steps of 0.25, the Lagrangian radius $r_L = R r_B$ where $R$ was varied between 1 and 2 in steps of 0.2, and various height functions and scaleheights were tested.

**The mass-to-light ratio.** The fiducial value of the stellar M/L ($\Upsilon_{3.6}$) is calculated as explained in Section 3.1.1. We then change this value by assuming different fractions of the fiducial $\Upsilon_{3.6}$ according to

$$ \Upsilon = f_d \Upsilon_{3.6}. \qquad (10) $$

We then re-calculate the disc rotation curve of the galaxy each time and add a dark matter halo with parameters such that the rotation curve in the outer parts of the galaxy matches the value predicted by the TF relation, as shown in Fig. 3.

**The Lagrangian radius.** By varying the Lagrangian radius while keeping all the other parameters of a model constant, we effectively vary the pattern speed of the bar. There are theoretical arguments that show that self-consistent bars cannot have a corotation radius smaller than the bar length (Athanassoula 1980; Contopoulos 1980), i.e. $R$ cannot be smaller than 1. On the other hand, what controls
the upper limit of the Lagrangian radius is the angular momentum transfer in the galaxy. Therefore, the slow-down of the bar will depend on the amount of resonant material available that is able to absorb angular momentum from the inner parts of the discs (Athanassoula 2003).

Previous studies of the shape of dust lanes in gas response sim-
ulations have set limits to $\mathcal{R}$, placing the corotation radius $r_L$ =
$(1.2 \pm 0.2)r_B$ (Athanassoula 1992b); however, there is still no con-
sensus in the literature regarding the value of $\mathcal{R}$, and indeed whether
it is the same for both early- and late-type barred spirals.

*The height function.* The scaleheight of the disc plays an important role in the strength of the forces in the plane of the galaxy (Fragkoudi et al. 2015). We vary the height function of the disc by assuming two functional forms for the height functions – isothermal and exponential. Additionally, we explore the effect of adding a b/p height function on the results, i.e. changing the scaleheight at certain locations in the disc. For the models with the flat height functions, we vary the scaleheight around this value – from $z_0 = 0.5$ to $z_0 = 1.5$ kpc (and the equivalent exponential scaleheight).

4 RESULTS

To quantitatively compare the models to the data, we place pseudo-
slits perpendicular to the shock loci in both the gas flow simulations
and the observed image, where the modelled gas flow is rotated
to match the rotation of the observed image. As shown in Fig. 4,
we trace the density jumps perpendicular to the slit, where the
maxima in the density indicate the location of the shocks. The

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Comparison between the location of the shocks in the models and the observations. A pseudo-slit is placed perpendicular to the shock in the simulation (top left) and the observed, unsharp, masked 8 µm image (top right). In the bottom panel, we show the density along the slit in the two images as a function of distance *l* along the slit. The shock loci are defined by the maximum density distribution of the simulation (vertical solid line) and the 8 µm image (vertical dashed line), and Δl is defined as the distance between the two shock loci.</figcaption></figure>

distance between the location of the shocks in the simulations and
observations is denoted by Δl, which is a proxy to the goodness of
fit of the models.

4.1 Effect of M/L

Decreasing the *M/L* of the stellar disc has a significant impact on the strength and shape of the shocks formed in the gas, as can be seen in Fig. 5, where we decrease the *M/L* from 1.5 to 0.25 × γ<sub>3.6</sub> in steps of 0.25. By reducing the *M/L* and accordingly increasing the contribution of the dark matter halo the shocks become progressively weaker, and the shape of the shocks becomes rounder. Additionally, the amount of gas present in the bar region increases since the bar is not strong enough to deplete the gas there. This occurs since we are effectively reducing the non-axisymmetric perturbation in the model, and increase the axisymmetric component. The strength of the non-axisymmetric perturbation determines the strength of the shocks, since by reducing the bar potential we effectively increase the extent of the *x*₂ family of periodic orbits, while the *x*₁ orbits become less cusped; therefore, the gas tends to stay on more circular orbits.

4.2 Effect of Lagrangian radius

By changing the pattern speed of the bar, the location of the La-
grangian (or more loosely, the corotation) radius in the disc changes.
In Fig. 6, we show the effect of changing the Lagrangian radius,
while leaving all other parameters in the models the same. By
increasing the Lagrangian radius, i.e. reducing the pattern speed
of the bar, the distance of the shocks from the centre of the bar
increases – in accordance with previous results in the literature
(Athanassoula 1992b; Sánchez-Menguiano et al. 2015) – while the
shocks formed at the leading edge of the bar are also stronger. For
very high pattern speeds, the shocks move towards the inside of
where the dust lanes occur (see model 007 in Fig. 6). The strength
of the shocks also decreases for high pattern speeds and therefore
the density contrast in the gas also decreases.

We already start to see an issue between changing the Lagrangian
radius and the *M/L*, in that there are similarities in the ways changing
the *M/L* and the pattern speed affect the gas flow. In both cases,
i.e. by decreasing the *M/L* and the Lagrangian radius, the shocks
are weakened and their shape becomes rounder, which leads to
degeneracies between the different models (for example, compare
model 002 of Fig. 5 and model 007 of Fig. 6).

4.3 Effect of height function

By increasing the scaleheight of the disc, we reduce the density in
the plane of the galaxy, which effectively reduces the strength of the
forces in the plane. This results in weaker shocks in the gas, since the
ratio of the non-axisymmetric-to-total forcings is also smaller. This
can be seen in Fig. 7 where we plot the gas response simulations in
five models with increasing scaleheight (where the scaleheight of
the model is noted in the upper left corner of the plot). The shocks
become rounder, less extended, and there is more gas present in
the bar region, since the bar cannot deplete the gas. There is also
less of a central concentration in the responding gas in the models
with thicker discs, since, due to the weakening of the bar, less gas
is funnelled to the central regions.

We also explore the effects of an exponential height function
on the shocks, since it is possible that by assigning an exponential
height function to the models, the results could be altered. We show
<figure><img src="image_8.png" /><figcaption>**Figure 5.** The gas response to models with different *M/L*: the *M/L* decreases from the top left to the bottom right from 1.5 to $0.25\Upsilon_{3.6}$ (as indicated in the top left corner of the plots). The rest of the parameters are the same for each model, namely the height function is an isothermal disc with scaleheight $z_0 = 0.8$ kpc and Lagrangian radius $r_L = 1.2r_B$. The blue diagonal line shows the location of the pseudo-slit in the simulation, and the red dots show the loci of the shock in the 8 µm *Spitzer* image of NGC 1291. The insets show the density along the pseudo-slit for both the simulation (thick black line) and the 8 µm *Spitzer* image (thin red line), and are provided to guide the eye on where the strongest shocks occur with respect to the pseudo-slit. The vertical lines in the insets show the shock loci along the slit for the models (solid line) and the observed image (dash-dotted line).</figcaption></figure>

in Fig. 8 the effects of assuming an exponential height function for the models as compared to an isothermal height function, where both models have an equivalent scaleheight (the exponential scaleheight $h_z = z_0/2$, where $z_0$ is the isothermal scaleheight). As expected, the shocks in the exponential model are slightly stronger. This occurs because for the exponential height function, the density distribution is more peaked in the mid-plane of the galaxy (see Fig. 2), which results in larger forces in the plane of the galaxy, and hence stronger shocks. After running a large number of models with an exponential height function, we found that the range of acceptable values for *M/L*, $r_L$ and $h_z$ is similar for both the exponential and isothermal height functions (see Section 4.4).

When a b/p height function is added to the model, the reduction of the non-axisymmetric forces is localized around the area where the b/p is maximum; this occurs due to the thickening of the disc mainly at the locations where the b/p is maximum. We model two different b/p geometries, where the maxima occur at $r = 2$ and 1.5 kpc, respectively, along the bar. In Fig. 9, we show the effect of adding these b/p bulge models on the rotation curve of the disc (second row), where we see that the b/p bulges reduce the height of the rotation curve in the region where the b/p is maximum. The effect of the b/p bulge on the shape of the shocks is also evident in the gas response to these models, as can be seen in the bottom row of

Fig. 9, where the shocks become rounder and have discontinuities. Also the area within the bar region is not as depleted of gas as the model without the b/p, and there is less gas in the central regions of the model with the b/p bulge (for more details on the gas inflow in models with and without b/p bulges, see Fragkoudi et al. 2016).

From Figs 5 and 7, we can see how increasing the scaleheight of the disc has a similar effect on the gas flow as decreasing the *M/L*; the shocks become weaker and rounder. Therefore, we see that there is a threefold degeneracy i.e. there is a degeneracy between the height function, the *M/L* and the Lagrangian radius of the models, although the degeneracy between the *M/L* and scaleheight is strongest. This is one of the difficulties of this type of study, and one of the ways we deal with these degeneracies is by using the velocity dispersion method in order to restrict the allowed parameter space, as will be discussed in the next sub-section.

## 4.4 Best-fitting models

Out of the different suites of models that were run we construct a number of grids for which we create contour plots showing the value of the distance $\Delta l$, between the dust lane in the 8 and 24 µm images and the maximum of the density in the gas dynamic simulations. For a given scaleheight and height function, we make a grid
<figure><img src="image_9.png" /><figcaption>Figure 6. As in Fig. 5 but for different Lagrangian radii (M/L is equal to $\Upsilon_{3.6}$ in all the models).</figcaption></figure>

of models with different M/L and different Lagrangian radii (i.e. pattern speeds) as is shown in Fig. 10.

To find the best-fitting models to the observations, we followed a procedure of elimination by applying the following criteria:

(i) The location of the primary shocks should be as close as possible to the dust lanes in the image (determined as shown in Fig. 4).

(ii) The model should have a dominant strong shock, i.e. one which is stronger than any secondary shocks formed in the gas flow.

(iii) The shape of the shocks should match the dust lanes in the 8 µm image and should also reproduce features in the 24 µm image.

We choose the models with the smallest $\Delta l$, from the contour plots in Fig. 10, from the allowed parameter space of models. This reduces the number of best-fitting models to 20 models. We then examine the models visually in order to determine the goodness of fit according to the set criteria.

Additionally, we combine this method with the velocity dispersion method in order to help break any remaining degeneracies. We use the relation between the scaleheight, velocity dispersion and surface density of an isothermal sheet in order to constrain the relation between the scaleheight and the M/L. This has been used in a number of studies in the past (e.g. Bottema 1993; Bershady et al. 2010; Martinsson et al. 2013). Assuming the disc is isothermal, the relation is given by

$$z_0 = \frac{\sigma_z^2}{\pi G \Upsilon I}, \quad (11)$$

where I is the surface intensity in units of flux per unit area and $\Upsilon$ is the M/L (Binney & Tremaine 2008; Bershady et al. 2010).

From Bosma et al. (2010), we take a value for the stellar velocity dispersion corresponding to the location of the pseudo-slit. This allows to constrain the parameter space, and excludes some of the models that have an M/L incompatible with the scaleheight predicted by the relation. It is possible that in the bar region the disc is not exactly isothermal, which is why we allow for a range of values (±25 per cent) for the M/L given a certain scaleheight. Thus, this relation can exclude some extreme cases in which for an adopted M/L the scaleheight is very far off from that given by the theoretical value of the above relation. This allowed region is given by the grey-shaded region in the plots in Fig. 10.

The four best-fitting models are shown in Fig. 11 together with the images of NGC 1291 in 8 and 24 µm for comparison. They all have a scaleheight of 0.8 kpc (equivalent to 0.4 kpc for the models with an exponential height function), while models that are even thinner than this tend to give a worse fit. In general, we expect the thickness of the disc to be between 0.8 and 1.2 kpc, since also for 1.2 kpc the models give worse fits than for 0.8 kpc. In all the models, fast bars are preferred with a value of $R \le 1.4$, since for slower bars the shocks tend to be generally too offset from the bar to correspond to the locations of the dust lanes.

The best-fitting models are summarized in Table 2. They are 001 and 107 that correspond to $z_0 = 0.8$ kpc without a b/p bulge and have M/L = 1 and $1.25\Upsilon_{3.6}$, respectively, and Lagrangian radius $r_L = 1.2$. Model 139 is also a relatively good fit to the dust lanes and has a b/p bulge with a length of 1.5 kpc along the bar, with an
<figure><img src="image_10.png" /><figcaption>**Figure 7.** As in Fig. 5 but changing the scaleheight of the models. The other parameters of the models are the same, namely $M/L = \Upsilon_{3.6}$ and $R = 1.2$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparing the gas morphology of models with an isothermal (left) and exponential (right) height function.</figcaption></figure>

$M/L = 1.25\Upsilon_{3.6}$ and $r_L = 1.4$. Additionally, model 909 with an exponential height function and $h_z = 0.4$ kpc (equivalent to $z_0 = 0.8$ kpc), $M/L = 0.75$ and $r_L = 1.2$ also gives a good fit to the observations.

It should be noted that none of the models give a *perfect* fit to the morphology of the dust lanes. This could be due to a number of reasons, such as the fact that we symmetrize the model although the galaxy is not perfectly symmetric, or due to small effects from inclination (i.e. even though we assume the galaxy is face-on, we know that it has a small, but perhaps non-negligible inclination). Additionally, although a large number of simulations were run, the grid we explore is still quite ‘coarse’ and there are regions of the parameter space which have not been explored in detail. It is

possible also, for example, that there is a boxy, rather than a strongly peanut-shaped, bulge.

# 5 DISCUSSION

## 5.1 The height function of NGC 1291

The best-fitting models predict a disc with scaleheight $\approx 0.8$ kpc for NGC 1291. The thicker the scaleheight of the model, the larger the $M/L$ that needs to be adopted in order to be able to reproduce the shocks in the gas. However, due to the fact that there is an inverse proportional relation between the scaleheight and the $M/L$ via the vertical velocity dispersion, the $M/L$ cannot be arbitrarily increased for a thick disc, since this would lead to a higher stellar velocity dispersion than that which is observed. Additionally, by increasing the scaleheight of the disc, the shape of the rotation curve of the model changes; therefore, even by increasing the $M/L$, the shock loci of the gas do not reproduce well the shapes of the dust lanes, which become too round for thicker discs.

An interesting point to note is that for borderline maximal disc cases such as NGC 1291, the functional form of the height function (i.e. whether we use an exponential or an isothermal) can have a significant effect on the results. Using an exponential can lead to a lower $M/L$, while using an isothermal can lead to a higher one. This, combined with the fact that the scaleheight can affect the results even more than the functional form, shows that a lot of attention must be paid to the vertical structure of the disc when carrying out this kind of study.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Effect of changing the length of the b/p bulge on the disc rotation curves and gas response of the models: Top row: the variation of the scaleheight of the model along the x-axis for model 001 that does not contain a b/p bulge, model 034 that contains a b/p bulge with length 2 kpc and model 066 that contains a b/p bulge with length 1.5 kpc (left, middle and right columns, respectively). Second row: the rotation curve of the three models. Bottom row: the gas surface density of the models.</figcaption></figure>

NGC 1291 seems to have a barlens, which hints at the existence of a b/p bulge (Athanassoula et al. 2013; Laurikainen et al. 2014). However, out of the best-fitting models shown in Fig. 11, only one contains a b/p bulge. There is some evidence that the shocks in the *central-most* part of the galaxy are rounder than the shape of the gas shocks for models without a b/p bulge. This hints at the possibility that there is a thicker structure in the central region, one which weakens the shocks thus making them rounder.<fnref n="2" /> It is therefore possible that a thicker region, for example due to a boxy-shaped bulge, rather than a strong peanut- or X-shaped bulge, exists in NGC 1291.

Indeed, since we have only explored a small fraction of the parameter space for parameters of the b/p bulge, such as the length, width and height of the b/p, it definitely cannot be excluded that NGC 1291 could contain a weak b/p bulge. However, the results seem to suggest that the galaxy does not contain a *strong* b/p bulge, i.e. one with a large double peaked maximum, however a weaker peanut or boxy-shaped bulge, or indeed a peanut bulge with smaller length and width, cannot be excluded.

## 5.2 Dynamical determination of the 3.6 µm M/L

The best-fitting models obtained in this study have an *M/L* that falls within the range predicted by SPS models and other approaches, which predict a value of $\Upsilon_{3.6} = 0.6 M_{\odot} L_{\odot}^{-1}$ (assuming a Chabrier IMF) with small variations for age and metallicity (Meidt et al. 2014; Norris et al. 2014; Röck et al. 2015).

SPS methods are highly dependent on the assumed IMF – indeed the choice of IMF can change the *M/L* by up to a factor of 10 (Bell & de Jong 2001) and in more recent works by up to a factor of 4 (Röck et al. 2015). On the other hand, since our technique provides a dynamical determination of the *M/L*, and as such is independent of IMF, it provides an independent confirmation of the robustness of these methods. Our best-fitting models have an *M/L* that lies in the range of $M/L_{3.6} = 0.6–0.75$, in line with the uncertainty (i.e. about 0.1dex) present in the above-mentioned methods (Meidt et al. 2014).

## 5.3 Maximal or sub-maximal disc?

We examine the rotation curves of all the best-fitting models to explore whether the disc of NGC 1291 is maximal. Model 001 has a disc that contributes 74 per cent of the rotation curve of the galaxy at $r = 2.2 h_r$. Therefore, according to the definition by Sackett (1997) the disc of this galaxy is just short of maximum. Both models 107 and 139 have a disc that contributes 80 per cent of the total circular velocity at $r = 2.2 h_r$, which therefore makes

<footnote n="2">² As shown in the appendix, one should be cautious when examining the central region of the simulations, since this region is dependent on the grid resolution of the simulations.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Contour plots of $\Delta l$ (i.e. the goodness of fit) for suites of models with different scaleheights (the scaleheight is denoted on the top of each plot). Note that the smaller the $\Delta l$ (i.e. lighter colours), the better the goodness of fit. The bottom right contour plot corresponds to an exponential height function with $h_z = 0.4$ kpc (for the sake of brevity, we do not include all the suites of models run with an exponential height function). The x-axis of the plots corresponds to M/L (in units of $\Upsilon_{3.6}$). The y-axis gives $\mathcal{R}$, i.e. the ratio of the Lagrangian radius to bar semimajor axis, $r_L/r_B$. The black dots correspond to points on the grids where we have models. The shaded area of the plots shows the constrained parameter space after taking into account the vertical velocity dispersion of the galaxy (see the text).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Two-dimensional gas surface density for the best-fitting models. In the left column, the top row shows model 001 and the bottom row 107 that have $z_0 = 0.8$ kpc, $M/L = 1.0$ and $1.25\Upsilon_{3.6}$, respectively, and $r_L = 1.2$. The second column shows the 8 and 24 µm images of NGC 1291 for comparison. In the left column, the top row shows model 139 that has $z_0 = 0.8$ kpc with a b/p with length = 1.5 kpc along the bar, $M/L = 1.25$ and $r_L = 1.4$, while the bottom row shows model 969 that has an exponential height function with $h_z = 0.4$ kpc, $M/L = 0.75$ and $r_L = 1.2$ ($r_L$ is given in units of $r_B$).</figcaption></figure>

<table><caption>**Table 2.** Best-fitting models: from left to right: the model number, the $M/L$ (in units of $\Upsilon_{3.6}$), the ratio of the Lagrangian to bar radius, the type of height function used (either isothermal or exponential), the scaleheight $z_0$ in kpc (or the equivalent $h_z$ in the case of an exponential profile), the contribution of the disc circular velocity to the total circular velocity at $r = 2.2r_h$ (above 75 per cent is maximal), and whether or not the model contains a b/p bulge.</caption><thead><tr><th>Model</th><th><i>M/L</i> (Υ<sub>3.6</sub>)</th><th><i>R</i></th><th>HF</th><th><i>z</i><sub>0</sub> (kpc)</th><th>Max (per cent)</th><th>b/p</th></tr></thead><tbody><tr><td>001</td><td>1</td><td>1.2</td><td>iso</td><td>0.8</td><td>74</td><td>no</td></tr><tr><td>107</td><td>1.25</td><td>1.2</td><td>iso</td><td>0.8</td><td>80</td><td>no</td></tr><tr><td>139</td><td>1.25</td><td>1.4</td><td>iso</td><td>0.8</td><td>80</td><td>yes</td></tr><tr><td>909</td><td>0.75</td><td>1.2</td><td>exp</td><td>0.8</td><td>63</td><td>no</td></tr></tbody></table>

them both maximal. The disc of model 909 contributes 63 per cent of the total rotation velocity at $r = 2.2h_r$, and thus would be dubbed sub-maximal. However, it is worth noting that the scalelength of NGC 1291 is unusually large (5.8 kpc); thus, even for model 909, the disc contributes most of the circular velocity in the inner regions of the galaxy and at least until the end of the bar (which is at $r = 5.1$ kpc). Thus, all our best-fitting models seem to be consistent with a maximal disc in the inner regions of the galaxy, and at least up to the end of the bar, although model 909 is slightly sub-maximal according to the definition by Sackett (1997).

These results are consistent with most previous results in the literature using a similar approach, i.e. gas dynamical modelling. In the work by Weiner et al. (2001), they show that the galaxy they study, NGC 4123 ($V_{max} \approx 140$ km s$^{-1}$), has a maximum disc. Zánmar Sánchez et al. (2008) find that the disc of NGC 1365 ($V_{max} \approx 300$ km s$^{-1}$) is massive although not quite maximal. The work of Kranz et al. (2003) showed that the two most massive

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Plot of the disc mass fraction as a function of maximum rotation velocity. We translate our value of $f_d$ to the equivalent of that used in Kranz et al. (2003), and add it to the plot (red circle). We see that NGC 1291 falls on the upper right side of the plot, consistent with galaxies with $V_{max} \ge 200$ having maximal discs. The blue squares correspond to previous works on the $M/L$ of barred galaxies using gas response simulations, i.e. Lindblad et al. (1996), Weiner et al. (2001), Weiner et al. (2004) and Zánmar Sánchez et al. (2008). Adapted from Kranz et al. (2003).</figcaption></figure>

galaxies in their sample (those with $V_{max} \ge 200$ km s$^{-1}$) had maximal discs, while the less massive galaxies were sub-maximal (see Fig. 12). It seems therefore that there is some tension between the results found in Weiner et al. (2001) and those in Kranz et al. (2003), since NGC 4123 should have a sub-maximal disc according to the findings of Kranz et al. (2003). We are consistent with these findings, i.e. the disc of NGC 1291, which is a galaxy with
$V_{\text{max}} \ge 200 \text{ km s}^{-1}$, is maximal. This seems to argue that massive barred galaxies tend to have massive discs as most of the aforementioned studies also find (e.g. Athanassoula et al. 1987). Whether or not this also extends down to lower masses remains an open question.

Our results are inconsistent with some other methods for dynamically determining the *M/L*, such as for example the DiskMass survey (Martinsson et al. 2013), which make use of the stellar vertical velocity dispersion. In their study, they found that all of the 30 galaxies used for the study have sub-maximal discs, even though about a third of their sample are spiral galaxies with $V_{\text{max}} \ge 200 \text{ km s}^{-1}$. This inconsistency could be due to the fact that the velocity dispersion probed in these studies corresponds to a younger (and therefore colder) stellar population, which leads to a smaller velocity dispersion and subsequently a smaller *M/L* (Aniyan et al. 2016). On the other hand, we are consistent with the findings of Bosma et al. (2010) for NGC 1291 who found, by exploring the vertical velocity dispersion as well as the photometric properties of the bar, that the disc of the galaxy is maximal.

## 5.4 Fast- or slow-rotating bar?

The three best-fitting models for NGC 1291 have fast-rotating bars with $R \le 1.4$. This is consistent with observational estimates of the bar pattern speed for early-type galaxies (Elmegreen 1996; Corsini 2011), although there are indications that bars in late-type galaxies might rotate slowly (Rautiainen, Salo & Laurikainen 2008). The results are also in agreement with theoretical values predicted by Athanassoula (1992b), who used the shape of the dust lanes to constrain $R$ and found that it should be of the order of $1.2 \pm 0.2$.

Additionally, this is consistent with previous gas-dynamical studies of other barred galaxies, which are also found to have fast-rotating bars, such as NGC 4123 (Weiner et al. 2001) and NGC 1365, which was found to have a fast-rotating bar by two independent studies ($r_L \approx 1.2r_B$, Lindblad et al. 1996; Zánmar Sánchez et al. 2008).

This result is also consistent with the fact that NGC 1291 is compatible with the maximum disc hypothesis, since there are theoretical arguments that suggest that galaxies embedded in concentrated haloes will lose angular momentum and thus slow down (Athanassoula 2003).

## 5.5 Additional physical processes

The simulations used in this study are kept as 'stripped-down' as possible in order to avoid significant computational costs due to the large number of simulations that need to be run in order to explore the three-dimensional parameter space as finely as possible. We therefore trade detailed physics (such as gas cooling, self-gravity, star formation and feedback), three dimensions and very high resolution in order to be able to run a larger number of simulations. Here, we briefly discuss the effects of leaving out these physical processes on our results.

As previously mentioned, and as is shown in our resolution study in the appendix, the resolution should not play a major role on the shape and loci of the shocks (to within 5 per cent accuracy). However, as was already shown by Sormani et al. (2015), the central regions are affected quite severely by the grid resolution, which is why we exclude these regions from the analysis in this work.

Fujimoto et al. (2014) show that in three dimensions, cold gas in the bar region can have significant vertical velocity dispersions (see their fig. 4). This vertical displacement could indeed have an effect on the shape of the shocks; however, we expect this effect

to be minimal, since the radial and tangential velocities in the disc are higher than the vertical velocity dispersion in the gas (see for example figs 10 and 11 in Athanassoula 1992b).

By including more detailed physics in our simulations, such as gas cooling and self-gravity, the gas would fragment and develop clumps (Nimori et al. 2013; Fujimoto et al. 2014), while star formation and feedback from supernovae can then disperse these dense molecular clouds (Dobbs, Pringle & Burkert 2012; Renaud et al. 2013). While including these recipes into our simulations would change the detailed structure of the shocks in the dust lane, we do not expect the shock loci to be affected globally. Additionally, the correct implementation of star formation and feedback recipes is a complex issue, well beyond the scope of this paper.

Therefore, considering the large amount of simulations needed to explore the parameter space, we think that the reduced realism of our simulations will not compromise the conclusions of this work; however, we plan to further explore the effects of more detailed physics and three-dimensional gas flows in future work.

# 6 SUMMARY

We explore the dynamical structure of NGC 1291, by using the gas in the galaxy as a tracer of the underlying gravitational potential, in order to constrain the disc *M/L* and therefore the amount of dark matter in its central regions. The potential of the galaxy is obtained from the 3.6 and 4.5 µm *Spitzer* images, and has three free parameters, namely the height function of the galaxy, the *M/L* and the bar pattern speed (or equivalently, the Lagrangian radius). We run hydrodynamic gas response simulations in these potentials, and explore the three-dimensional parameter space, to obtain the best-fitting models, by matching the shock loci in the gas response simulations to the loci of the observed dust lanes in the 8 and 24 µm images.

There are a number of interesting trends that emerge from this study: the effects of the scaleheight and functional form of the height function on the gas flows in hydrodynamic gas response simulations have not been extensively studied in the past. We show that by increasing the scaleheight of the disc, or by employing a less peaked height function, the forces in the plane of the galaxy are reduced, which leads to weaker shocks with a rounder morphology. This effect is similar to that which occurs when the *M/L* of the model is decreased, while increasing the contribution of the dark matter halo to maintain the outer part of the rotation curve at a value predicted by the TF relation. In both cases, i.e. when we increase the scaleheight or decrease the *M/L*, we are decreasing the strength of the non-axisymmetric forcings in the plane of the galaxy – albeit in a slightly different way in each case. Furthermore, in accordance with previous results in the literature (Athanassoula 1992b; Sánchez-Menguiano et al. 2015), we show that by decreasing the pattern speed of the bar, the shocks in the gas are displaced from the bar semimajor axis. Conversely, by increasing the pattern speed of the bar the morphology of the shocks becomes rounder and, for very fast bars, the shocks are weaker.

We find that there are degeneracies between the three main free parameters that make it difficult to unambiguously find the best-fitting model for the galaxy. We therefore combine the gas response method with the velocity dispersion method in order to further constrain the parameter space.

The best-fitting models obtained for NGC 1291 suggest the following:

(i) The *M/L* of the disc at 3.6 µm is within the range given by SPS and other methods, i.e. $\Upsilon_{3.6} = 0.6 M_{\odot}/L_{\odot}$ (Meidt et al. 2014;
Norris et al. 2014; Röck et al. 2015). Since our method does not rely on any assumptions about the IMF, it provides a separate independent measurement of $Υ_{3.6}$.

(ii) This *M/L* suggests that NGC 1291 has a maximal disc (three out of four best-fitting models have maximal discs). This confirms previous findings using similar methods, i.e. that massive barred galaxies have maximal discs (Lindblad et al. 1996; Weiner et al. 2001; Kranz et al. 2003).

(iii) NGC 1291 contains a fast bar, with $R \le 1.4$, also in agreement with previous results in the literature (Elmegreen 1996; Rautiainen et al. 2008; Corsini 2011).

(iv) The scaleheight of NGC 1291 is around the range of 0.8–1.2 kpc (for a sech² height function).

The results show that a wealth of information can be extracted from the dynamical modelling of barred galaxies via gas response simulations. We plan to expand this work to include other galaxies from the $S^4G$ survey in order to obtain a larger statistical sample and probe a broader range of masses. This should help shed some light on conflicting results in the literature, especially in the low-mass regime, where there remains considerable debate as to whether or not discs are maximal.

## ACKNOWLEDGEMENTS

FF would like to thank Miguel Querejeta for helpful discussions on obtaining mass maps from NIR images of galaxies. The authors would also like to thank Dimitri Gadotti for providing us with the structural parameters of NGC 1291. We acknowledge financial support to the Detailed Anatomy of Galaxies network from the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme FP7/2007-2013/ under REA grant agreement number PITN-GA-2011-289313. We also acknowledge financial support from the CNES (Centre National d'Etudes Spatiales – France). This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the programme 'Investissements d'Avenir' supervised by the Agence Nationale pour la Recherche.

## REFERENCES

* Aniyan S., Freeman K. C., Gerhard O. E., Arnaboldi M., Flynn C., 2016, MNRAS, 456, 1484

* Athanassoula E., 1980, A&A, 88, 184

* Athanassoula E., 1992a, MNRAS, 259, 328

* Athanassoula E., 1992b, MNRAS, 259, 345

* Athanassoula E., 2003, MNRAS, 341, 1179

* Athanassoula E., 2008, IAU Symp. 245, Formation and Evolution of Galaxy Bulges. p. 93

* Athanassoula E., Bosma A., Papaioannou S., 1987, A&A, 179, 23

* Athanassoula E., Machado R. E. G., Rodionov S. A., 2013, MNRAS, 429, 1949

* Bell E. F., de Jong R. S., 2001, ApJ, 550, 212

* Bendo G. J. et al., 2007, MNRAS, 380, 1313

* Bendo G. J. et al., 2010, MNRAS, 402, 1409

* Bershady M. A., Verheijen M. A. W., Swaters R. A., Andersen D. R., Westfall K. B., Martinsson T., 2010, ApJ, 716, 198

* Binney J., Tremaine S., 2008, Galactic Dynamics, 2nd edn. Princeton University Press, Princeton, NY

* Bissantz N., Gerhard O., 2002, MNRAS, 330, 591

* Bissantz N., Englmaier P., Gerhard O., 2003, MNRAS, 340, 949

* Bosma A., 1978, PhD thesis, Groningen Univ.

* Bosma A., 1981, AJ, 86, 1825

* Bosma A., Gadotti D. A., de Blok W. J. G., Athanassoula E., 2010, in Verdes-Montenegro L., Del Olmo A., Sulentic J., eds, ASP Conf. Ser. Vol. 421, Galaxies in Isolation: Exploring Nature Versus Nurture, Kinematics of the Barred Spiral Galaxy NGC 1291. Astron. Soc. Pac., San Francisco, CA, p. 53

* Bottema R., 1993, A&A, 275, 16

* Bregman J. N., Hogg D. E., Roberts M. S., 1995, ApJ, 441, 561

* Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486

* Contopoulos G., 1980, A&A, 81, 198

* Corsini E. M., 2011, Mem. Soc. Astron. Ital. Suppl., 18, 23

* Courteau S., Rix H.-W., 1999, ApJ, 513, 561

* de Grijs R., 1998, MNRAS, 299, 595

* de Grijs R., Peletier R. F., van der Kruit P. C., 1997, A&A, 327, 966

* de Vaucouleurs G., de Vaucouleurs A., Corwin H. G., Jr, Buta R. J., Paturel G., Fouqué P., 1991, Third Reference Catalogue of Bright Galaxies. Volume I: Explanations and references. Volume II: Data for galaxies between 0$^h$ and 12$^h$. Volume III: Data for galaxies between 12$^h$ and 24$^h$. Springer, New York, NY

* Debattista V. P., Sellwood J. A., 2000, ApJ, 543, 704

* Dobbs C. L., Pringle J. E., Burkert A., 2012, MNRAS, 425, 2157

* Elmegreen B., 1996, in Buta R., Crocker D. A., Elmegreen B. G., eds, ASP Conf. Ser. Vol. 91, IAU Colloq. 157: Barred Galaxies, Pattern Speeds in Barred Galaxies. Astron. Soc. Pac., San Francisco, p. 197

* Eskew M., Zaritsky D., Meidt S., 2012, AJ, 143, 139

* Fragkoudi F., Athanassoula E., Bosma A., Iannuzzi F., 2015, MNRAS, 450, 229

* Fragkoudi F., Athanassoula E., Bosma A., 2016, MNRAS, 462, L41

* Fujimoto Y., Tasker E. J., Wakayama M., Habe A., 2014, MNRAS, 439, 936

* Gadotti D. A., 2008, MNRAS, 384, 420

* Gnedin O. Y., Weinberg D. H., Pizagno J., Prada F., Rix H.-W., 2007, ApJ, 671, 1115

* Hamadache C. et al., 2006, A&A, 454, 185

* Hogg D. E., Roberts M. S., Bregman J. N., Haynes M. P., 2001, AJ, 121, 1336

* Kennicutt Jr, R. C. et al., 2003, PASP, 115, 928

* Kim W.-T., Seo W.-Y., Stone J. M., Yoon D., Teuben P. J., 2012, ApJ, 747, 60

* Kranz T., Slyz A., Rix H.-W., 2001, ApJ, 562, 164

* Kranz T., Slyz A., Rix H.-W., 2003, ApJ, 586, 143

* Kuijken K., Gilmore G., 1991, ApJ, 367, L9

* Laurikainen E., Salo H., Athanassoula E., Bosma A., Herrera-Endoqui M., 2014, MNRAS, 444, L80

* Lindblad P. A. B., Lindblad P. O., Athanassoula E., 1996, A&A, 313, 65

* McGaugh S. S., 2012, AJ, 143, 40

* McGaugh S. S., Schombert J. M., Bothun G. D., de Blok W. J. G., 2000, ApJ, 533, L99

* Martínez-Valpuesta I., Athanassoula E., 2008, in Knapen J. H., Mahoney T. J., Vazdekis A., eds, ASP Conf. Ser. Vol. 390, Pathways Through an Eclectic Universe, Boxy/Peanut Bulges and Stellar Bars. Astron. Soc. Pac., San Francisco, CA, p. 463

* Martinsson T. P. K., Verheijen M. A. W., Westfall K. B., Bershady M. A., Andersen D. R., Swaters R. A., 2013, A&A, 557, A131

* Meidt S. E. et al., 2014, ApJ, 788, 144

* Nimori M., Habe A., Sorai K., Watanabe Y., Hirota A., Namekata D., 2013, MNRAS, 429, 2175

* Norris M. A., Meidt S., Van de Ven G., Schinnerer E., Groves B., Querejeta M., 2014, ApJ, 797, 55

* Patsis P. A., Athanassoula E., 2000, A&A, 358, 45

* Pérez I., Fux R., Freeman K., 2004, A&A, 424, 799

* Prendergast K. H., 1983, in Athanassoula E., ed., Proc. IAU Symp. 100, Internal Kinematics and Dynamics of Galaxies, Theoretical Studies of Gas Flow in Barred Spirals Galaxies. D. Reidel Publishing Co., Dordrecht, p. 215

* Querejeta M. et al., 2015, ApJS, 219, 5

* Rautiainen P., Salo H., Laurikainen E., 2008, MNRAS, 388, 1803

* Regan M. W. et al., 2006, ApJ, 652, 1112

* Renaud F. et al., 2013, MNRAS, 436, 1836

* Roberts M. S., 1976, Comments Astrophys., 6, 105

* Röck B., Vazdekis A., Peletier R. F., Knapen J. H., Falcón-Barroso J., 2015, MNRAS, 449, 2853

* Sackett P. D., 1997, ApJ, 483, 103
Sánchez-Menguiano L., Pérez I., Zurita A., Martínez-Valpuesta I., Aguerri
J. A. L., Sánchez S. F., Comerón S., Díaz-García S., 2015, MNRAS,
450, 2670

Schinnerer E. et al., 2013, ApJ, 779, 42

Sheth K. et al., 2010, PASP, 122, 1397

Sormani M. C., Binney J., Magorrian J., 2015, MNRAS, 449, 2421

Teyssier R., 2002, A&A, 385, 337

Thilker D. A. et al., 2007, ApJS, 173, 538

Tielens A. G. G. M., 2008, ARA&A, 46, 289

Tisserand P. et al., 2007, A&A, 469, 387

Trott C. M., Webster R. L., 2002, MNRAS, 334, 621

Tully R. B., Rizzi L., Shaya E. J., Courtois H. M., Makarov D. I., Jacobs B. A., 2009, AJ, 138, 323

van Albada T. S., Bahcall J. N., Begeman K., Sancisi R., 1985, ApJ, 295, 305

van der Kruit P. C., Searle L., 1981, A&A, 95, 105

Weiner B. J., 2004, in Ryder S., Pisano D., Walker M., Freeman K., eds,
Proc. IAU Symp. 220, Dark Matter in Galaxies, The Dark Matter density
Problem in Massive Disk Galaxies. Astron. Soc. Pac., San Francisco,
p. 265

Weiner B. J., Sellwood J. A., Williams T. B., 2001, ApJ, 546, 931

Zánmar Sánchez R., Sellwood J. A., Weiner B. J., Williams T. B., 2008, ApJ, 674, 797

central regions of the gas response simulations is dependent on
the adopted grid resolution. We wanted to explore whether this
would affect our ability to draw robust conclusions from the type
of analysis carried out in this paper. To do this, we ran tests using a
number of analytic models containing a disc, bar and bulge, taken
from Athanassoula (1992a). Our aim was to determine the effects
of resolution on the shocks induced in the gas flows.

We found, in agreement with Sormani et al. (2015), and as can be
seen in Fig. A1, that the size of the central disc (sometimes referred
to as the ‘x₂’ disc) is reduced for higher resolutions, presumably due
to the gas being able to better resolve the x₁ orbits. We measured
the extent of the central disc for the different resolutions, and found
that the disc size converges for higher resolutions (see Fig. A2),
which is also in agreement with what was found by Sormani et al.
(2015).

We evaluated the effect of resolution on the shock loci, by plac-
ing pseudo-slits perpendicular to the shocks, similarly to Fig. 4.
We found that although the central regions of the models are sig-
nificantly affected, the shock loci are the same for the different
resolutions, to within 5 per cent. This was found for a large number
of models.

As can be seen in Fig. A1, we also found that at high resolution,
the gas starts developing instabilities, again in agreement with the
findings in Sormani et al. (2015). There is no consensus as to why
the flow becomes unstable at high resolutions, but it is a feature
common to all hydrodynamic codes, and could be connected with
the onset of turbulence due to high shear (Kim et al. 2012).

APPENDIX: TESTS ON SIMULATIONS: THE
EFFECT OF GRID RESOLUTION

Recent work by Sormani et al. (2015) using gas response simula-
tions, signalled a worrying trend, i.e. that the morphology of the

<figure><img src="image_16.png" /></figure>

**Figure A1.** Gas surface density with *x*₁ and *x*₂ orbits overplotted. From top left to bottom right we increase the resolution, as noted at the top of each panel. The first panel shows the gas surface density for a Cartesian grid with a resolution 40 pc, where we have a refinement of 2¹⁰ cells per side. For each subsequent panel, we increase the resolution by a factor of 2 and employ an adaptive mesh refinement (AMR) scheme. For the plots titled AMR, the resolution indicated is the maximum resolution reached by the AMR grid.
<figure><img src="image_17.png" /></figure>

**Figure A2.** Size of the $x_2$ disc (along the y-axis), as a function of grid resolution for two models, 001 and 052 from Athanassoula (1992a). Model 052 has a slightly stronger bar. The size of the $x_2$ disc seems to converge for higher resolutions. For model 001, the size of the disc seems to increase for higher resolutions but this could be also due to instabilities in the gas.

We found therefore that for the purposes of this study, the shock loci are a more robust diagnostic than the central $x_2$ disc, since it is not dependent on the resolution of the simulations. For high enough resolutions, however, the central gaseous disc seems to also converge and therefore, for high enough resolutions, this could also be used as another constraint when comparing gas-dynamical models with observations.

This paper has been typeset from a TeX/LATEX file prepared by the author.