# The effects of Boxy/Peanut bulges on galaxy models

F. Fragkoudi, E. Athanassoula, A. Bosma, F. Iannuzzi

► To cite this version:

F. Fragkoudi, E. Athanassoula, A. Bosma, F. Iannuzzi. The effects of Boxy/Peanut bulges on galaxy models. Monthly Notices of the Royal Astronomical Society, 2015, 450 (1), pp.229–245. 10.1093/mn-ras/stv537 . hal-01439786

HAL Id: hal-01439786

https://hal.science/hal-01439786

Submitted on 15 Sep 2021

**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_1.png" /></figure>

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

# The effects of Boxy/Peanut bulges on galaxy models

F. Fragkoudi,\* E. Athanassoula, A. Bosma and F. Iannuzzi

Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388 Marseille, France

Accepted 2015 March 9. Received 2015 March 2; in original form 2014 September 8

## ABSTRACT

We examine the effects that the modelling of a Boxy/Peanut (B/P) bulge will have on the estimates of the stellar gravitational potential, forces, orbital structure and bar strength of barred galaxies. We present a method for obtaining the potential of disc galaxies from surface density images, assuming a vertical density distribution (height function), which is let to vary with position, thus enabling it to represent the geometry of a B/P. We construct a B/P height function after the results from a high-resolution, N-body+SPH simulation of an isolated galaxy and compare the resulting dynamical model to those obtained with the commonly used, position-independent ‘flat’ height functions. We show that methods that do not allow for a B/P can induce errors in the forces in the bar region of up to 40 per cent and demonstrate that this has a significant impact on the orbital structure of the model, which in turn determines its kinematics and morphology. Furthermore, we show that the bar strength is reduced in the presence of a B/P. We conclude that neglecting the vertical extent of a B/P can introduce considerable errors in the dynamical modelling. We also examine the errors introduced in the model due to uncertainties in the parameters of the B/P and show that even for generous but realistic values of the uncertainties, the error will be noticeably less than that of not modelling a B/P bulge at all.

**Key words:** galaxies: bulges – galaxies: kinematics and dynamics – galaxies: structure.

## 1 INTRODUCTION

Many edge-on disc galaxies can be seen to contain boxy-, peanut- or X-shaped isophotes, which are usually grouped together into the category of Boxy/Peanut (hereafter B/P) bulges. As a result of a number of theoretical studies (see Athanassoula 2008b, 2015 for general reviews of this subject), including orbital structure and stability analysis (Binney 1981; Pfenniger 1984, 1985; Patsis, Skokos & Athanassoula 2002; Skokos, Patsis & Athanassoula 2002a,b; Martinez-Valpuesta, Shlosman & Heller 2006; Harsoula & Kalapotharakos 2009; Contopoulos & Harsoula 2013), and numerical simulations (Combes & Sanders 1981; Combes et al. 1990; Raha et al. 1991; Mihos et al. 1995; Athanassoula & Misiriotis 2002; Athanassoula 2003, 2005; Bureau & Athanassoula 2005; Martinez-Valpuesta et al. 2006), these structures are now known to be due to vertical instabilities in the bar, which cause it to ‘puff up’, giving rise to boxy or peanut-like shapes. These studies also show that once a bar forms, a B/P bulge will form soon after.

Observational studies have further confirmed the link between B/P bulges and bars (see Kormendy & Kennicutt 2004 and Kormendy 2015 for reviews on the subject), by showing that the fraction of edge-on disc galaxies with B/P bulges is comparable to the fraction of disc galaxies containing bars (Lütticke, Dettmar

& Pohlen 2000). Kinematic studies of edge-on barred galaxies and B/P bulges also confirm the connection between the two structures (Athanassoula & Bureau 1999; Bureau & Athanassoula 1999, 2005; Chung & Bureau 2004, and references therein). Therefore, barred galaxies at present and past epochs will contain B/P bulges, and in fact, one is also believed to be present in our own Galaxy (Weiland et al. 1994; Howard et al. 2009; McWilliam & Zoccali 2010; Shen et al. 2010; Ness et al. 2012, Ness et al. 2013a,b; Nataf et al. 2013, 2015; Vásquez et al. 2013; Wegg & Gerhard 2013; Gardner et al. 2014; Nataf, Cassisi & Athanassoula 2014).

Bars are found in about two thirds of disc galaxies in the local Universe, with variable strengths (Eskridge et al. 2000; Menéndez-Delmestre et al. 2007; Barazza, Jogee & Marinova 2008; Aguerri, Méndez-Abreu & Corsini 2009; Gadotti 2009), and are known to be the main drivers of the secular phase of galaxy evolution. The torque they induce into the disc causes outward angular momentum transfer, which in turn will cause a redistribution of matter in the disc. They are thus responsible for driving gas to the centre of their host galaxy (Athanassoula 1992b), forming discy pseudo-bulges (Kormendy & Kennicutt 2004; Athanassoula 2005), redistributing stars in the galactic disc (Sellwood & Binney 2002; Roškar et al. 2008; Minchev & Famaey 2010), and possibly creating a fuel reservoir for AGN activity (Shlosman, Begelman & Frank 1990; Coelho & Gadotti 2011; Emsellem et al. 2015, but see also Lee et al. 2012, for a review see Combes 2001). However, even though the effect of bars on all these processes has been thoroughly

<footnote>\* E-mail: francesca.fragkoudi@lam.fr</footnote>
examined, a study of the effects of B/P bulges on all these processes has not, until present, been carried out.

As a first step towards understanding the effect B/P bulge geometry may have on the aforementioned processes, in this paper we focus on their influence on dynamical models of their host galaxies. These models are obtained directly from images of the galaxies' surface brightness, by first assuming a vertical density distribution, or height function, and subsequently deriving the potential of the galaxy. They have been used extensively in the literature, with one of their most important implementations being in simulations that study the response of gas in a fixed potential. These response simulations are used to study the dark matter content of galaxies and to test the maximum disc hypothesis (Kranz, Slyz & Rix 2001, 2003; Weiner, Sellwood & Williams 2001; Slyz, Kranz & Rix 2003; Pérez, Fux & Freeman 2004), the bar pattern speed (Lindblad, Lindblad & Athanassoula 1996; Kalapotharakos, Patsis & Grosbøl 2010b) as well as the kinematical and morphological properties of gas in galaxies (Lin et al. 2011, 2013). Dynamical models have also been used in studies determining the bar strength (Buta & Block 2001; Laurikainen & Salo 2002) and the orbital structure of galaxies (Quillen, Frogel & Gonzalez 1994; Patsis, Athanassoula & Quillen 1997; Kalapotharakos, Patsis & Grosbøl 2010a; Patsis, Kalapotharakos & Grosbøl 2010). Furthermore, they have been used to study gravitational torques in barred and spiral galaxies in order to establish the amount of gas inflow and by extension determine the importance of secular evolution (Zaritsky & Lo 1986; Haan et al. 2009; Foyle, Rix & Zibetti 2010). In all of these aforementioned studies, the geometry of the B/P bulge is not taken into account when constructing the height function, and instead a position independent, 'flat' height function is assumed. This is partly due to the lack of an analytical model for a B/P bulge, as well as to the inherent difficulty of detecting these bulges in face-on or intermediate inclination galaxies, which are the galaxies generally used in these studies.

Various methods, however, have been proposed over the past few years, which allow either for the detection of B/P bulges, or at least for an educated guess at their existence. By viewing a large number of N-body+SPH simulations, and covering a wide range of viewing angles, Athanassoula et al. (2014) have shown that B/P bulges manifest themselves in face-on projections as the so-called barlens (Laurikainen et al. 2011), which renders their detection fairly easy. Strong observational arguments for this have been presented in Laurikainen et al. (2014). Another method proposed by Debattista et al. (2005), uses signatures in the stellar kinematics of face-on or almost face-on galaxies and was implemented by Méndez-Abreu et al. (2008) who confirmed the existence of a B/P bulge in NGC 98. Furthermore, it is possible to detect signatures of B/P bulges by examining the morphological features of inclined galaxies (Athanassoula & Beaton 2006; Erwin & Debattista 2013).

We therefore believe that a study of the effect of B/P bulges on models of their host galaxies, and by extension of the necessity of including the geometry of B/P bulges in the height function of these models, is called for. To this aim, we first introduce in Section 2 a straightforward and reliable method for calculating the potential, forces and derivatives of forces of a general density distribution $\rho(r, \phi, z)$. We present some tests which demonstrate that the method can give highly accurate results, while also allowing the flexibility to choose an arbitrary height function, without being restricted to one which is constant with position.

We then used our code on an image of an N-body+SPH simulated galaxy, which is presented in Section 2.3, thus obtaining a realistic potential for a barred galaxy. In order to create this model, we

assign a thickness and a height function to the galaxy. These height functions are introduced in Section 3 and include two 'flat' height functions, a function which describes peanut bulges (from which we construct our fiducial B/P model), and another which describes boxy bulges.

The main results are presented in Section 4, where we examine the effect B/P bulges have on the potential and forces (4.1), on the periodic orbits (4.2) and on the bar strength (4.3). We find that B/P bulges indeed have a significant effect on the results and therefore conclude that they should be included when modelling their host galaxy.

In Section 5, we explore the errors which will be induced in the results by using a B/P model which is not exactly the 'correct' one. This is necessary since it is not trivial to observationally obtain the exact parameters of B/P bulges, and this can introduce errors in the model. We show that for a range of uncertainties, the errors induced in the results are less than those induced by not modelling the B/P at all. We also introduce a new method for calculating the bar strength, $Q_{int}^{int}$, which takes into account the variation of the non-axisymmetric forcings along the whole extent of the bar.

In Section 6, we give a summary and list the main conclusions of our work.

## 2 METHOD AND TESTS

To create a dynamical model of a galaxy, we first need its density distribution. The two-dimensional surface density can be obtained from surface brightness images of a face-on disc galaxy by assuming an M/L ratio (see for example Querejeta et al. 2014). It is important that these images are taken in a wavelength range which minimizes the effect of dust extinction and traces the old stellar population (for example, the Spitzer 3.6 µm band). In this work, we use an image of a face-on simulated galaxy, and thus we do not need to account for dust extinction, nor assign an M/L ratio, as our two-dimensional image gives the surface density directly. Once we have the surface density, we also need to assign a height function, and together these give us the three-dimensional density distribution, from which we can calculate the potential of the galaxy due to the stellar component. Our method for calculating the potential involves a straightforward three-dimensional integration over the density distribution and we refer to it throughout the paper as the 3DF method. We calculate the potential in Cartesian coordinates by

$$ \Phi(x, y, z) = -G \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{\rho(x', y', z')}{\sqrt{\sum_{j=1}^{3}(x_j - x')^2 + \epsilon^2}} dx'dy'dz' \quad (1) $$

where G is Newton's gravitational constant, $\rho$ is the density and $\epsilon$ is the softening length which is necessary to eliminate the noise at the expense of a small bias (Merritt 1996; Athanassoula et al. 2000). We can differentiate the expression in equation (1) analytically with respect to x, y and z, to obtain expressions for the force and its derivatives. We thus rely heavily on an adequate integration algorithm, specifically one which can deal with singularities. To tackle this we use CQUAD, a doubly adaptive integration algorithm (Galassi et al. 2003), which requires more function evaluations than other integration routines, but is more successful in dealing with difficult integrands. It computes the integral within the desired relative error limits (or precision), which the user can set. Since we mainly work in the $z=0$ plane, we focus in what follows on the non-zero quantities in the plane: the potential $\Phi$ and the two non-zero
components of the force $F_x$ and $F_y$. As mentioned, the above, as well as what follows, concerns the potential and forces of the stellar component of the galaxy.

## 2.1 Tests on the method: relative errors

In order to test the accuracy of our method, we create a model of a barred galaxy containing a disc, a bar and a classical bulge, using density distributions which have analytic solutions for the potential and forces. We then calculate the potential and forces for this model using the 3DF method, and compare the results against the analytic solutions. The general results of these tests are very positive, which demonstrates the ability of our code to deal with difficulties such as cuspy and/or non-axisymmetric density configurations. For more details on the model of the galaxy, we refer the reader to Appendix A.

To calculate the relative errors of the potential and the derivatives of the force, we use the relation

$$ Error = 0.5 \frac{|R_1 - R_2|}{|R_1| + |R_2|}, \qquad (2) $$

where $R_1$ and $R_2$ are, respectively, the analytic and 3DF solutions.

To calculate the relative errors for the force, we use the relation

$$ Error = \frac{|F_{i1} - F_{i2}|}{\sqrt{F_{i1}^2 + F_{j1}^2}}, \qquad (3) $$

where i and j can be either the x or the y component of the force, and the subscripts 1 and 2 stand for the analytic and 3DF solutions, respectively. We therefore normalize the error of each component of the force by the total force at each point. This is done because our main interest in the forces is for the calculation of orbits and because on the symmetry axes of the x and y components of the force (in the static frame of reference), the analytic estimates of $F_x$ and $F_y$ will be equal to zero.

We stress that the precision with which the code calculates the results is an input parameter to the code. The accuracy can be as high as the user wants it to be (within the limits of machine precision), at the expense, of course, of computation time. We require a three dimensional integration and, due to the propagation of error at each integration, the relative precision we ask of the CQUAD algorithm for each integration has to be larger than that which we wish to achieve. Practically this means that if we ask for a relative precision of $10^{-3}$, we will obtain a relative precision of approximately $10^{-1}$. This is sufficient for our purposes as the error is less than 1 per cent for all variables. The softening is set to 10 pc throughout the paper. For this precision and softening, the maximum error of the potential is 0.3 per cent, of $F_x$ 0.6 per cent and of $F_y$ 0.7 per cent.

## 2.2 Tests on the method: orbits

Even though from the results of the relative errors we see that the 3DF method gives highly accurate results, we would like to confirm that the noise in the force field does not prevent orbits from running smoothly, as they would in an analytic potential. To do this, we calculate a number of periodic orbits in the analytic potential and in the potential derived using the 3DF method for the galaxy in the frame corotating with the bar, in the model described in Appendix A. The grid used for the orbits and throughout the paper is 200 × 200, and the orbits are calculated using a Kick-Drift-Kick leapfrog algorithm (Hockney & Eastwood 1988; Quinn et al. 1997; Springel 2005).

<figure><img src="image_3.png" /><figcaption>Figure 1. Comparison of orbits in the analytic and 3DF potentials of our model galaxy. The thin black line gives the outline of the bar. The orbits calculated in the 3DF potential are given in solid thick black lines, and the orbits calculated in the analytic potential are given in dashed red lines. We plot some $x_1$ orbits (along the bar), some $x_2$ orbits (perpendicular to the bar), and an almost circular orbit outside the bar region. We see that the orbits in the two potentials almost completely overlap, so that the red and black lines are practically indistinguishable.</figcaption></figure>

In Fig. 1, we plot a few of these orbits. In this figure and all throughout the paper the bar major axis is along the x-axis. In the figure, we show some $x_1$ orbits, which are elongated along the bar, some $x_2$ orbits which are perpendicular to the bar, as well as some nearly circular orbits outside the bar region. We see that the orbits calculated in the 3DF potential are a very good approximation of those calculated in the analytic potential, as the two practically coincide. Thus the error which is introduced in the potential from our 3DF method and the adopted value of precision is sufficient for our purposes.

## 2.3 The image

We use the 3DF method on the density distribution derived from a face-on image of a simulated isolated galaxy and different height functions (which are described in Section 3). The initial conditions of the simulation from which the image was constructed, include a live spherical dark matter halo, an exponential stellar disc and 75 per cent gas (for more information on the simulation the reader is referred to run gtr116 in Athanassoula, Machado & Rodionov 2013). The snapshot we use is taken well into the secular evolution phase of the galaxy, specifically at 8 Gyr after the start of the simulation, by which point a strong bar and B/P bulge have formed. The image we use is constructed from the 'stars' component of the snapshot and has a morphology reminiscent of that of many strongly barred galaxies, such as IC 4290 (see Fig. 2).

In order to decrease the noise of the image, we require a snapshot with a large number of particles. We create a snapshot with 40 times the particles of the original snapshot, following the procedure described in Athanassoula (2005). To further reduce the noise
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Visual comparison between IC 4290 and gtr116. The two galaxies have striking morphological similarities and are classified as having the same bar strength (for more details see Section 4.3).</figcaption></figure>

in the image we apply some smoothing, by Fourier decomposing and recomposing it.

The Fourier components are calculated as follows.

$$a_n(r) = \frac{1}{\pi} \int_{-\pi}^{\pi} \Sigma(r, \theta) \cos(n\theta) \, d\theta, \qquad (4)$$

$$b_n(r) = \frac{1}{\pi} \int_{-\pi}^{\pi} \Sigma(r, \theta) \sin(n\theta) \, d\theta, \qquad (5)$$

where $a_n$ and $b_n$ are the even and odd Fourier components, $\theta$ is the azimuthal angle, $r$ the radius and $\Sigma$ gives the surface density. We then reduce the high frequency noise by recomposing the image as

$$\Sigma(r, \theta) = \frac{a_0}{2} + \sum_{m=2}^{m=n_F} (a_m(r)\cos(m\theta) + b_m(r)\sin(m\theta)), \qquad (6)$$

using only a limited number of even Fourier components (in our case $n_F = 26$). We show in Figs 3(a) and (b) the surface density of the original image and of the Fourier recomposed image, respectively, both in arbitrary units, and in Fig. 3(c) we show the residual image of the two. As the images of surface density are in arbitrary units, the density, as well as the potential and its derivatives will also be in arbitrary units in what follows.

## 3 HEIGHT FUNCTIONS USED

In order to obtain the three-dimensional density of a galaxy disc from a two-dimensional image, we need to assume a height function, which defines how the density drops off as a function of $z$ from the

equatorial plane $z = 0$. The height function and the scaleheight ($z_0$) will of course affect the results, and we therefore need to use the height function which best approximates that of the galaxy we are trying to model.

The height function can be either constant or can change with position. In the case where it is constant with respect to position we assume, for simplicity, that the density distribution can be written as

$$\rho(x, y, z) = \Sigma(x, y)F(z), \qquad (7)$$

where $\rho$ is the three-dimensional density distribution, $\Sigma$ is the two-dimensional surface density, and $F$ is the height function. In the more general case, where the height function depends on position in the galaxy, as would be for example the height function describing a B/P bulge, the scaleheight changes as a function of position. In this case, the density distribution would be given by

$$\rho(x, y, z) = \Sigma(x, y)F(x, y, z), \qquad (8)$$

where the normalization of the height function is

$$\int_{-\infty}^{\infty} F(x, y, z) \, dz = 1. \qquad (9)$$

It is worth noting that the mass of the model is always the same; the height function simply determines the volume density of the galaxy, by setting the thickness of the disc.

### 3.1 Flat height functions

Up to now in the literature, position-independent or ‘flat’ height functions have been used when modelling barred disc galaxies. We therefore also use two flat height functions in this paper, to check the discrepancy which will be created in the model by (a) using a flat height function and a B/P height function, and (b) using two different flat height functions. We adopt two commonly used functions, the isothermal-sheet model (van der Kruit & Searle 1981):

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

and the sech-law model (van der Kruit 1988):

$$F(z) = \frac{1}{\pi z_0} \mathrm{sech} \left( \frac{z}{z_0} \right), \qquad (11)$$

where $1/(2z_0)$ and $1/(\pi z_0)$ are the respective normalization factors.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left: original image showing the surface density of the stellar component of the gtr116 simulation. Middle: model from the Fourier recomposition, using up to $n_F = 26$ even Fourier components. Right: residual image after subtracting the model from the original image.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Left: side-on image of the surface stellar density of the simulated galaxy gtr116. Middle: the scaleheight of the simulation (red crosses) is plotted along the x-axis (for y = 0, i.e. the side-on projection). The solid black line shows the fit of z₀(x, 0) to the data, which gives the scaleheight of the fiducial peanut height function. Right: plot of the scaleheight of the simulation (red crosses) along x = 3 which is where the maximum of the peanut occurs (end-on projection). The solid black line shows the values of z₀(x, y) along x = 3.</figcaption></figure>

## 3.2 Peanut height function

To obtain a height function for the peanut, we examined the particle distribution along different cuts in x and y from the simulation introduced in the previous section. We found that the sum of two two-dimensional Gaussians for the scaleheight can provide a reasonable approximation to the B/P shape. As can be seen in Fig. 4 and as commented below, this choice may fail at certain points, but provides an overall fair representation of the structure.

The resulting B/P height function is a non-separable function of position and is given by

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

The scaleheight $z_0(x, y)$ varies like the sum of two two-dimensional Gaussians:

$$z_0(x, y) = A \exp \left( - \left( \frac{(x - x_0)^2}{2\sigma^2} + \frac{(y - y_0)^2}{2\sigma^2} \right) \right) + A \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 (13)$$

where A is the maximum scaleheight of the peanut above the disc scaleheight, $z_0^{\mathrm{disc}}$. The variance of the Gaussians is given by $\sigma^2$, ($x_0$, $y_0$) is the position of the maximum of the first Gaussian and ($x_1$, $y_1$) the position of the maximum of the second Gaussian. We fit these two two-dimensional Gaussians to values of the scaleheight obtained from the simulation along y = 0 and x = 3 (which is where the maximum of the scaleheight occurs). In the remainder of the paper, we refer to this as our fiducial peanut (or B/P) model.

To obtain the scaleheights, we take cuts along the x- and y-axes and fit the vertical particle distribution with a sech² function. We thereby determine the variation of $z_0$ from bin to bin along the cut. The results can be seen in Fig. 4. In the side-on view (panel b), we see that the scaleheight along y = 0 behaves approximately like the combination of two Gaussians, except in the central region where the scaleheight drops below that of the outer disc. For a cut along x = 3, where the peanut is maximum (end-on view, panel c), the behaviour of $z_0$ is still well approximated by a Gaussian, although our fit slightly under predicts the value of the scaleheight.

Along some cuts at x values intermediate between the centre and the peanut maximum, the Gaussian approximation fails to represent the behaviour of the scaleheight with y. In fact the behaviour of the scaleheight in the presence of a B/P bulge is quite complex, and cannot be grasped entirely by a simple analytic function. However,

as it turns out, the fitted function shown in Fig. 4 underestimates the value of $z_0$ at these points. This directly translates into an underestimation of the effect of the peanut in those regions. In summary, our fiducial model for the peanut height function shown in Fig. 4 will result into a conservative estimate of the effect of the real peanut present in the image we adopt as our starting point. Given the scope of this paper, which is to demonstrate the generic effect of a peanut bulge on a galaxy model, we find this approximation more than satisfactory.

## 3.3 Boxy height function

The B/P bulge might at times have rather boxy isophotes. This could be due to projection effects, whereby the peanut is projected at such an angle that the isophotes appear boxy (Athanassoula & Misiriotis 2002). However, boxy isophotes might be present even when the bar is seen side-on, i.e. they might be the real shape of the B/P bulge (see Patsis et al. 2002 for a discussion based on orbits). This tends to be the case for galaxies with weak bars, where instead of a strong x-shape or peanut forming, boxy isophotes are seen (Athanassoula 2008a).

To model a boxy bulge, we use a height function which drops off as sech² with height from the z = 0 plane,

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

where the scaleheight is a top-hat function,

$$z_0(x, y) = \begin{cases} z_0^{\mathrm{bulge}} & |x| \le x_{\mathrm{max}} \ \& \ |y| \le y_{\mathrm{max}} \\ z_0^{\mathrm{disc}} & \text{otherwise.} \end{cases}, \qquad (15)$$

and where $z_0^{\mathrm{disc}}$ gives the scaleheight of the disc and $z_0^{\mathrm{bulge}}$ gives the scaleheight, or strength, of the boxy bulge. This is quite a simplified model of the boxy bulge, with only two free parameters, its strength and length (which is set by $L = 2x_{\mathrm{max}}$). The thickness of the box, i.e. $y_{\mathrm{max}}$, is set by the width of the bar.

We create the fiducial boxy height function such that it best approximates the fiducial peanut height function, so such that we can examine whether the former can be used as an approximation for the latter, as there is one less parameter to model. The fiducial boxy height function therefore has a height equal to the height of the fiducial peanut model and its length is such that the boxy bulge finishes approximately where the peanut scaleheight is in between its maximum and minimum (see top-right panel of Fig. 5).
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Errors from not taking into account the proper B/P geometry: the top row gives the scaleheight (in red and black lines) along the bar major axis for the setups we are comparing in the plots. The second, third and fourth rows give the relative difference between the two setups being compared for the potential, Fₓ and Fᵧ, respectively (see the colourbar for values of the relative difference). The dark green line represents the ellipse fitted to the outer isophote of the bar. First column: difference between the sech and sech² setup. Second column: difference between the fiducial peanut height function and a sech² height function. Third column: difference between a boxy height function and a sech² height function. Fourth column: difference between a boxy height function and our fiducial peanut height function. We see that not including a peanut or a boxy bulge where there is one will induce large errors in the potential and forces and also that a boxy height function is not a good approximation for a peanut height function. For details see the text (Section 4.1).</figcaption></figure>
# 4 BOXY/PEANUT OR NO BOXY/PEANUT?

We wish to investigate whether accounting for the geometry of the B/P in the height function will significantly change the model of its host galaxy, and therefore whether we should include it in the modelling when it is present. Thus, in the next three subsections we investigate the effect B/Ps will have on the potential and forces (Section 4.1), on the periodic orbits (Section 4.2) and on the bar strength (Section 4.3) of the model.

## 4.1 B/P effect on potential and forces

We calculate the potential and the forces for the density distribution given by the image described in Section 2.3, and the different height functions described in Section 3. The results of this subsection are shown in Fig. 5. In the top row, we plot the scaleheight along the bar major axis for the two setups we are comparing in the plots. We show two-dimensional plots of the relative difference for the potential in the second row, of the *x* component of the force in the third row, and of the *y* component of the force in the fourth row. The green line represents an ellipse fitted to the outer isophote of the bar. Each column gives one of the following comparisons (from left to right).

(i) *Two flat height functions: sech and sech²*

We compare the models obtained by implementing two flat height functions, sech and sech², with equivalent scaleheights, in order to demonstrate that different flat height functions do not significantly affect the results. In the very centre, the difference for the potential is only around 1 per cent and for the forces it is 5 per cent, while for the rest of the grid the difference between the two setups is well below 1 per cent in all cases. If we decrease the value of the scaleheight the two height functions produce even more similar results. This happens because as the disc tends to become infinitesimally thin, the shape of the height function becomes less important. Equivalently, if we increase the value of the scaleheight, and hence the thickness of the disc, then the difference in the results obtained with the two height functions increases.

We see therefore that the scaleheight, and not the vertical profile of the height function, is primarily responsible for creating differences in the models.

(ii) *Flat and peanut height functions*

We compare a flat height function and the height function of our fiducial peanut bulge, i.e. a peanut with parameters fitted to our simulated galaxy. The differences that arise from using these two setups are significant for the potential and forces, as can be seen in the second column of Fig. 5. This is especially true near and around the region of the maximum height of the peanut, and in general in and around the bar. The force can be different in the two cases by up to 40 per cent, which is not surprising since around the maximum of the peanut the scaleheight is more than three times the value of the scaleheight of the disc. This can be seen in the top row of the figure, which demonstrates how the scaleheight varies along *x* for the two height functions. The larger scaleheight reduces the forces in the plane of the disc, due to a reduction of the density in the plane. Therefore, we see that by not taking into account the geometry of the B/P bulge, we induce significant errors in the model, i.e. in the potential and its derivatives.

(iii) *Flat and boxy height functions*

In the third column of Fig. 5, we compare a flat height function and the height function of the fiducial boxy bulge. We see that a boxy height function will also induce large differences compared to the flat height function, and in fact in the central regions our fiducial

boxy height function has an even larger effect than the peanut. Boxy bulges are usually associated with weaker bars in simulations and are therefore typically less strong than peanut bulges (although at early times boxy bulges can be as strong as peanut bulges – see figs 2 and 3 in Athanassoula & Misiriotis 2002). In observations, boxy bulges can appear as strong as peanut bulges (see for example Chung & Bureau 2004) although it is hard to distinguish whether these are truly boxy bulges, or simply peanut bulges seen at an angle. It is therefore reasonable to assume that a substantial amount of boxy bulges will be somewhat less strong than our fiducial boxy height function. However, for the sake of simplicity and to be able to compare our fiducial boxy height function with our fiducial peanut height function, we give the former the same strength as the latter. Therefore, our fiducial boxy height function can be thought of as an upper limit for the effect of a boxy bulge on the model of its galaxy.

(iv) *Peanut and boxy height functions*

We compare our fiducial peanut height function to our fiducial boxy height function in order to see to what extent the boxy height function approximates the peanut height function as it has one less free parameter than the peanut height function. These results can be seen in the fourth column of Fig. 5. For the potential and forces the match between the boxy height function and the peanut height function is quite poor, especially in the central region where the scaleheights of the two height functions are very different. Therefore, the boxy height function is not a good approximation to a peanut bulge.

## 4.2 B/P effect on periodic orbits

In this section, we examine how some of the most important families of periodic orbits will be influenced by taking into account the geometry of a B/P. To do this we study two models: the model with the fiducial peanut bulge height function, and the model with the sech² height function *without* a B/P bulge. We set the pattern speed of both models to be such that corotation occurs just outside the bar radius, within the range $1.4 > R_{CR}/R_{bar} > 1$, where $R_{CR}$ and $R_{bar}$ are the corotation and bar radius, respectively (e.g. Athanassoula 1992b). The orbits are calculated in a frame of reference corotating with the pattern speed of the bar.

In Fig. 6(a), we show a few typical orbits in the bar region for the potentials we are examining, overplotted on the image of our simulated galaxy, gtr116, shown face-on. The three most important families of orbits in the bar region are shown, i.e. the $x_1$ (red lines, extended along the bar major axis), $x_2$ (cyan lines, perpendicular to the bar major axis) and 3/1 (green lines, asymmetric with respect to the y-axis) families, which are stable along most of their extent.

In Fig. 6(b), we plot the characteristic diagram of periodic orbits, for the two cases with and without a peanut. The characteristic diagram gives the value at which the orbit intersects the y-axis ($y_0$) as a function of its Jacobian energy ($E_j$, i.e. energy in the rotating frame of reference; Binney & Tremaine 2008). The Jacobian energy is in arbitrary units, since, as already mentioned, the mass is also in arbitrary units. We see that the characteristic diagram of the two models differs significantly. The most noticeable effect due to the presence of a B/P is the change in the bifurcation loci of the upper and lower branch of the 3/1 family. This indicates that taking into account the geometry of a B/P in the model changes the location of the 3:1 resonance, and therefore the 3/1 family of periodic orbits appears at higher energies. Thus orbits of the 3/1 family will differ in the two cases, as can be seen in Fig. 7(b), where we plot two 3/1 orbits in the two models for the same cut along the y-axis. The extent of the 3/1 family of orbits is also significantly increased for the case with a B/P bulge, surpassing the extent of the $x_1$ family of
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Left: some typical orbits in the bar region for the model with a sech² height function, overplotted on the image of the simulated galaxy gtr116: In red the $x_1$ bar supporting orbits, in cyan an $x_2$ orbit perpendicular to the bar and in green the 3/1 orbits (asymmetric with respect to the y-axis). Right: characteristic diagram (intersection of each orbit with the y-axis as a function of the Jacobi energy) for the models created from the image of the simulated galaxy gtr116 and the two height functions. The solid black line gives the characteristic diagram for the model with the fiducial peanut bulge and the dashed red line the characteristic diagram for a model with a flat sech² height function. The dotted blue line shows the zero velocity curve (ZVC) for the sech² model (the ZVC of the two models are very similar).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Left: two $x_1$ orbits with the same Jacobian energy ($E_J = -0.9$) calculated in the two potentials: with (solid black line) and without (dashed red line) a B/P bulge. We see that in the B/P model the $x_1$ orbit's length and height are reduced (its extent along the x-axis is reduced by ~12 per cent and along the y-axis by ~46 per cent – measured, respectively, at $y = 0$ and $x = 0$). Right: two 3/1 periodic orbits for $y_0 = 1.2$ in the two potentials (colours as before). Again the orbits differ significantly.</figcaption></figure>

periodic orbits, which is in fact shorter compared to the $x_1$ family in the model without a B/P bulge.

The $x_1$ family also suffers changes, in the $E_J$ region between -1.1 and -0.8. In this area of the diagram, the maximum extent of the orbits along the x-axis reaches the region where the effect of the B/P is maximum; therefore, for these energies the orbits of the two models differ. In Fig. 7(a), we show an $x_1$ orbit of the same energy ($E_J = -0.9$) plotted in the two models. When a B/P is present, the maximum extent of the orbit along the x-axis is reduced by 12 per

cent, while its maximum extent along the y-axis is reduced by 46 per cent (measured at $x = 0$).

For the $x_2$ family, the highest (lowest) value of $y_0$ increases (decreases) for the model with a B/P bulge (see Fig. 6b), i.e. the entire extent of the $x_2$ family is increased by about 43 per cent. As the extent of the $x_2$ family is related to the distance between the two Inner Lindblad Resonances (ILRs; Athanassoula 1992a), the distance between these two resonances will therefore also increase. This increase is due to a weakening of the non-axisymmetric
<table><caption>**Table 1.** Errors in bar strength.</caption><thead><tr><th>Comparison</th><th>⟨Error Q<sub>T</sub>⟩</th><th>MAX(Error Q<sub>T</sub>)</th><th>Q<sub>b</sub></th></tr></thead><tbody><tr><td>sech–sech<sup>2</sup></td><td>1 per cent</td><td>1.6 per cent</td><td>1.5 per cent</td></tr><tr><td>peanut–no peanut</td><td>27 per cent</td><td>74 per cent</td><td>4 per cent</td></tr><tr><td>peanut–boxy</td><td>14 per cent</td><td>42 per cent</td><td>16 per cent</td></tr></tbody></table>

Note. We show the average and maximum of the relative errors of $Q_T$, as well as the relative error of $Q_b$, for three different comparisons of the setups.

perturbation: when the geometry of B/Ps is taken into account in the model, the scaleheight of the galaxy is increased, where the B/P is maximum, and therefore – since the amount of mass is the same – the volume density in the plane of the galaxy is decreased. This leads to a decrease in the radial and tangential forces in such a way that the non-axisymmetric perturbation is diminished, thus changing the distance between the two ILRs. This is in accordance with results from both Contopoulos & Grosbol (1989) and Athanasoula (1992a), the latter of which showed that there are a number of model parameters which can affect this distance. In particular, in figs 6 and 7 of Athanassoula (1992a), we see that the distance between the ILRs can increase due to a decrease in the bar mass or the pattern speed, or due to an increase in the central mass concentration of the galaxy or the axial ratio of the bar.

The differences between the orbital families of the two models will have effects on their stellar, as well as their gaseous dynamics. The extent of the $x_2$ orbits plays a crucial role on the shape of the shock loci in the gas (Athanassoula 1992a, 1992b), and therefore the shape of the gas shocks in the two models should differ significantly; conversely, the shape and strength of the shocks influence the amount of gas inflow towards the centre of the galaxy, and it is therefore likely that there will be a measurable difference in the amount of gas inflow in models with and without a B/P bulge. This is further supported by the results in Section 4.3, where we show that B/P bulges reduce the strength of the bar; we plan to address in upcoming work the extent to which the gas flows will be affected.

## 4.3 B/P effect on bar strength

We study the effect of the B/P bulge on one of the measures of bar strength, which involves calculating the non-axisymmetric forcings on the disc due to the bar, i.e. the bar-induced torque (Combes & Sanders 1981; Buta & Block 2001). The magnitude of this non-axisymmetric perturbation is given by

$$Q_T(r) = \frac{F_T^{\max}(r)}{\langle F_R(r) \rangle}, \qquad (16)$$

where $F_T$ is the tangential force $F_T(r) = (1/r) (\partial\Phi/\partial\phi)$, and $\langle F_R(r) \rangle$ is the average over azimuth of the radial force $F_R(r) = \partial\Phi/\partial r$. The forces are calculated directly from the image, as described in Section 2. In order to obtain a single measure of the bar strength for a galaxy, the quantity $Q_b$, which is the maximum of $Q_T$ in the bar region, is commonly defined as the bar strength. In what follows we investigate the effect that a B/P height function will have on both $Q_b$ and $Q_T$.

The results of this study are discussed in paragraphs (i), (ii) and (iii), where we examine models with the flat, the fiducial peanut and the fiducial boxy height function, respectively. In Table 1, we show the maximum and average relative errors of $Q_T$, denoted MAX(Error $Q_T$) and $\langle Error Q_T \rangle$, respectively, as well as the relative error of $Q_b$,

when comparing two models with different height functions. The average relative error of $Q_T$ over radius is calculated according to

$$\langle \text{Error } Q_T \rangle = \frac{1}{n} \times \sum_{i=1}^{n} \left( \left| \text{abs} \left( \frac{Q_{T_1}(r_i) - Q_{T_2}(r_i)}{Q_{T_1}(r_i)} \right) \right| \times 100 \right). \quad (17)$$

Plots of these results can also be seen in Fig. 8, where it is worth noting that for the flat and peanut height functions, the strength of the bar is $Q_b \approx 0.55–0.6$. According to Buta & Block (2001), this represents a strong bar case (between bar class 5 and 6), which corresponds to approximately 20 per cent of their sample of SB galaxies. We have already shown in Fig. 2 the striking morphological similarity between IC 4290 and our galaxy, and we note here that IC 4290 is also classified by Buta & Block (2001) as a class 6 barred galaxy, with $Q_b = 0.56$. Therefore the results presented in this section, as well as in previous and subsequent sections, correspond straightforwardly and quantitatively to strongly barred galaxies. However, even weakly barred galaxies will have B/P bulges, albeit weaker ones, and therefore the results will also apply to these galaxies although to a lesser extent. We intend to carry out a full statistical study of the effects of different strength B/P bulges on the models of their host galaxy, together with a full comparison to observations, elsewhere.

### (i) A model with a sech and sech² height function

In previous work by Laurikainen & Salo (2002), the effect of position-independent height functions and height functions which only vary as a function of radius was examined, and these were found not to change $Q_b$ in a significant way. For realistic height functions such as the exponential, sech, or sech² models, they found that $Q_b$ was affected by less than 5 per cent, which is consistent with our own results. This is confirmed in Fig. 8 and Table 1, where we see that for an equivalent scaleheight, the height functions of sech and sech² will produce very similar bar strengths, which will tend to become even more similar the thinner the disc.

### (ii) A model with the fiducial peanut height function

We plot $Q_T$ for our fiducial B/P height function. In the region around the scaleheight maxima (at $r = 3$), $Q_T$ is significantly flatter than the model without a peanut, due to the reduction of the strength of the tangential and radial forces. The value of $Q_b$ will not be significantly different from the case with the flat height function, due to the fact that the maxima of the peanut and the maximum of $Q_T$ are at a relatively large distance from each other. The torque induced by the bar in the two cases however is significantly different, as can be seen both in the plot and in the second row of Table 1. By using the $Q_b$ method of measuring bar strength, the bars will be judged as having the same strength (class 6), and hence the same effect on the disc, even though the forces in the plane of the galaxy are significantly reduced in the presence of a peanut. Therefore, in Section 5.3 we introduce another measure of bar strength, which can capture the reduction in bar strength when a B/P bulge is present.

### (iii) A model with the fiducial boxy height function

We also plot $Q_T$ for our fiducial boxy height function. We see again that where the boxy bulge is maximum, $Q_T$ is flattened due to the decrease in the strength of the bar forces in the plane. We also see that where the top-hat boxy function ends, $Q_T$ exceeds the values of the sech and sech² curves. This is due to the effect of the boxy bulge on the tangential and radial forces, with the former increasing just outside the boxy bulge while the latter is decreased in the whole disc due to the overall decrease in mass–density in the plane of the galaxy. The combination of these two effects results in the torque becoming large in the region just outside the boxy bulge. As a
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Strength of non-axisymmetric forcings ($Q_T$) as a function of radius, for models with different height functions: sech (solid red line), sech² (dashed green line), the fiducial peanut setup (thick black solid line) and the fiducial boxy setup (dash-dotted blue line). The vertical solid black line indicates the radius at which the scaleheight of the fiducial peanut is maximum and the vertical dashed black line indicates the end of the bar.</figcaption></figure>

consequence, $Q_b$ is overestimated by 16 per cent compared to the fiducial peanut case (third row of Table 1). We investigate the boxy height function as an alternative to using the peanut height function – since there is one parameter less to model – and conclude that even if one is merely interested in $Q_T$, a simple boxy height function is not a good approximation to a peanut function.

# 5 ERRORS DUE TO BOXY/PEANUT MODELLING

In this section, we investigate how much error will be induced if we include a B/P bulge in the model, but with the main peanut parameters differing from that of our fiducial model. This type of error is induced due to observational uncertainties, as it is not trivial to observationally obtain the correct parameters for the B/P bulge we want to model. This is due to the physics of the problem, not the numerical part of the calculation (as is the error referred to in Section 2.1 which can be made arbitrarily small) and is therefore practically an unavoidable source of uncertainty. Nevertheless, as we will discuss below, there do exist empirical and theoretical arguments which can constrain the parameter space of a B/P bulge.

The height function we have chosen for our fiducial B/P bulge, the sum of two two-dimensional Gaussians, has three degrees of freedom. Thus inaccuracies in the modelling of the B/P bulge can also be introduced in three ways: by estimating wrongly the height of the Gaussians (which corresponds to a change in peanut strength), or the distance between the maxima of the Gaussians (which corresponds to a change in the peanut length), or the widths of the Gaussians (which corresponds to a change in peanut ‘width’, i.e. how peaked or thin the peanut is at its maximum).

## 5.1 Potential and forces

In this subsection, we investigate how much error is introduced in the potential and forces by incorrectly modelling the B/P bulge.

### 5.1.1 Peanut strength uncertainties

The maximum value of the scaleheight of the peanut, also called the peanut strength, is a value which is not trivial to find observationally. Numerical studies have shown that the strength of the peanut correlates with the bar strength (Athanassoula 2008a). However this relation has a considerable scatter, and can merely give an approximate estimate. Debattista et al. (2005) showed that for face-on, or nearly face-on N-body simulated galaxies, an indicator of the presence and strength of a B/P bulge is the fourth order Gauss-Hermite moment of the line-of-sight velocity dispersion, $h_4$. However, this relation has not been quantified in such a way which would allow a direct measurement of peanut strength from $h_4$. Studies of orbital structure (Patsis et al. 2002; Skokos et al. 2002a) have also suggested specific families of periodic orbits which are responsible for giving the B/P bulge its height, but no direct measurement of the strength of the B/P bulge is available from orbital structure either.

In order to measure the error due to peanut strength uncertainties, we use a grid of models with varying peanut strength and compare them to our fiducial peanut setup. These results can be seen in Table 2 below. In this and in all subsequent tables, the term ‘Average Error’, indicates the average error within the outer isophote of the bar and ‘Maximum Error’ corresponds to the maximum error found in the grid, excluding the central-most point. We see that an over- or underestimation of the error by the same amount will produce similar errors in the model. In all the cases studied, the error induced by an incorrect peanut strength is always less than that induced by not modelling a B/P bulge at all.

### 5.1.2 Peanut width uncertainties

In Table 3, we show the errors for a grid of models with different peanut width errors. For all cases considered, the error induced due to a miscalculation of the peanut width is less than that induced by not modelling a peanut at all, apart from the maximum error induced in the potential when the peanut width is 50 per cent larger than in the fiducial scenario. This is due to a sharp increase in scaleheight
<table><caption>**Table 2.** Percentage error due to peanut strength uncertainty.</caption><thead><tr><th rowspan="2">Peanut strength</th><th colspan="3">Average error</th><th colspan="3">Maximum error</th></tr><tr><th>Φ</th><th>F<sub>x</sub></th><th>F<sub>y</sub></th><th>Φ</th><th>F<sub>x</sub></th><th>F<sub>y</sub></th></tr></thead><tbody><tr><td>+50 per cent</td><td>1.3 per cent</td><td>2 per cent</td><td>3 per cent</td><td>3 per cent</td><td>10 per cent</td><td>6 per cent</td></tr><tr><td>+25 per cent</td><td>0.7 per cent</td><td>1 per cent</td><td>2 per cent</td><td>1.5 per cent</td><td>5 per cent</td><td>3 per cent</td></tr><tr><td>-25 per cent</td><td>0.7 per cent</td><td>1 per cent</td><td>2 per cent</td><td>1.5 per cent</td><td>6 per cent</td><td>4 per cent</td></tr><tr><td>-50 per cent</td><td>1.4 per cent</td><td>3 per cent</td><td>4 per cent</td><td>3 per cent</td><td>14 per cent</td><td>10 per cent</td></tr><tr><td>no peanut</td><td>3 per cent</td><td>8 per cent</td><td>10 per cent</td><td>7 per cent</td><td>37 per cent</td><td>28 per cent</td></tr></tbody></table>

*Note.* Average and maximum errors of the potential and forces for setups with different peanut strength error, within the area enclosed by the outer isophote of the bar. The last row gives the error induced by not including a B/P bulge in the model at all.

<table><caption>**Table 3.** Percentage error due to peanut width uncertainty.</caption><thead><tr><th rowspan="2">Peanut width</th><th colspan="3">Average error</th><th colspan="3">Maximum error</th></tr><tr><th>Φ</th><th>F<sub>x</sub></th><th>F<sub>y</sub></th><th>Φ</th><th>F<sub>x</sub></th><th>F<sub>y</sub></th></tr></thead><tbody><tr><td>+50 per cent</td><td>3 per cent</td><td>5 per cent</td><td>6 per cent</td><td>15 per cent</td><td>26 per cent</td><td>27 per cent</td></tr><tr><td>+25 per cent</td><td>1 per cent</td><td>2 per cent</td><td>3 per cent</td><td>4 per cent</td><td>10 per cent</td><td>12 per cent</td></tr><tr><td>-25 per cent</td><td>1 per cent</td><td>2 per cent</td><td>3 per cent</td><td>3 per cent</td><td>11 per cent</td><td>6 per cent</td></tr><tr><td>-50 per cent</td><td>2 per cent</td><td>5 per cent</td><td>6 per cent</td><td>5 per cent</td><td>11 per cent</td><td>15 per cent</td></tr><tr><td>no peanut</td><td>3 per cent</td><td>8 per cent</td><td>10 per cent</td><td>7 per cent</td><td>37 per cent</td><td>28 per cent</td></tr></tbody></table>

*Note.* Average and maximum errors of the potential and forces for setups with different peanut width errors within the area enclosed by the outer isophote of the bar. The last row gives the error that would be present if we do not model a B/P bulge at all.

in the central region for large peanut widths (see solid red line in Fig. 10c), which is where the potential is most affected. This error however is confined only to the potential and to the central most grid points, and should not have significant effects on orbital calculations in most of the galaxy.

### 5.1.3 Peanut length uncertainties

Of the three parameters – length, strength and width – length has the least uncertainty, due to a method proposed by Athanassoula et al. (2014). The method determines the length of B/Ps for face-on and moderately inclined galaxies, which uses the shape of the projected isophotes in the bar region. They demonstrated that the barlens and the peanut are the same component and therefore that the size of the former can be used to estimate the length of the latter. For galaxies with larger inclinations the length can be estimated from other morphological features in the isophotes created by the B/P bulge (Athanassoula & Beaton 2006; Erwin & Debattista 2013), while orbital structure studies confirm the aforementioned results and also give clues as to the length of the peanut (Patsis et al. 2002; Patsis, Skokos & Athanassoula 2003). Due to all this, the uncertainties of the length estimates are rather small, certainly smaller than the corresponding ones for strength and width, which is why we use a smaller range of uncertainties for the peanut length.

We carry out comparisons for a grid of models with different peanut length errors and give the results in Table 4. As expected, the more we change the length of the peanut away from the fiducial value, the larger the errors will be, although there is an asymmetry in the error induced with respect to over- and underestimating the length; by underestimating the peanut length by a certain amount, we induce more error than by overestimating it by the same amount. By decreasing the length of the peanut we induce more error in the central regions of the galaxy, which is where the potential is most affected, due to the two Gaussians overlapping in the centre and thus increasing the scaleheight (see dotted magenta line, Fig. 10e).

This can be seen in Table 4, for the case of −30 per cent peanut length where, for the potential, the maximum and average errors are larger than that of the +30 per cent peanut length case.

For all the cases considered, the average error induced is smaller or equal to that of not modelling the B/P bulge. Given that the length is a fairly well-constrained quantity, large errors are not expected to be present due to the peanut length in the modelling of the B/P bulge.

### 5.1.4 Combinations of uncertainties

It is likely that a combination of different kinds of error will contribute to the total error budget of a model of the B/P. It is not in the scope of this paper to explore the full parameter space of the possible error combinations, instead we choose a few cases in order to get a feel of the amount of error that can be induced. By 'combination of error', we refer to a combination of all the different sources of error. By '+50 per cent' ('−50 per cent') we refer to a setup with peanut strength and width which are 50 per cent larger (smaller) than the fiducial value, and a peanut length which is 30 per cent larger (smaller) than the fiducial. We do not find it necessary to further increase the error in peanut length, since it is the best constrained quantity out of the three parameters. The '+25 per cent' ('−25 per cent') setup corresponds to one with peanut strength, width and length which is 25 per cent larger (smaller) than the fiducial value. In Table 5, we see that all the combinations of uncertainties will introduce less error in the model than not including a B/P bulge at all.

## 5.2 Periodic orbits

As shown in Section 4.2, the presence of a B/P bulge will affect the extent and shape of the different families of periodic orbits which make up the bar. In this section, we qualitatively explore the errors introduced in the calculation of periodic orbits due to incorrect modelling of a B/P bulge. To do this, we examine the characteristic
<table><caption>**Table 4.** Percentage error due to peanut length uncertainty.</caption><thead><tr><th rowspan="2">Peanut length</th><th colspan="3">Average error</th><th colspan="3">Maximum error</th></tr><tr><th>Φ</th><th>F<sub>x</sub></th><th>F<sub>y</sub></th><th>Φ</th><th>F<sub>x</sub></th><th>F<sub>y</sub></th></tr></thead><tbody><tr><td>+30 per cent</td><td>1 per cent</td><td>6 per cent</td><td>5 per cent</td><td>4 per cent</td><td>24 per cent</td><td>13 per cent</td></tr><tr><td>+16 per cent</td><td>0.9 per cent</td><td>3 per cent</td><td>3 per cent</td><td>3 per cent</td><td>12 per cent</td><td>7 per cent</td></tr><tr><td>-16 per cent</td><td>1 per cent</td><td>4 per cent</td><td>4 per cent</td><td>5 per cent</td><td>11 per cent</td><td>10 per cent</td></tr><tr><td>-30 per cent</td><td>3 per cent</td><td>8 per cent</td><td>9 per cent</td><td>13 per cent</td><td>25 per cent</td><td>28 per cent</td></tr><tr><td>no peanut</td><td>3 per cent</td><td>8 per cent</td><td>10 per cent</td><td>7 per cent</td><td>37 per cent</td><td>28 per cent</td></tr></tbody></table>

*Note.* Average and maximum errors of the potential and forces for setups with different peanut length errors within the area enclosed by the outer isophote of the bar. The last row gives the error induced by not modelling a B/P bulge at all.

<table><caption>**Table 5.** Percentage error due to combination of uncertainties.</caption><thead><tr><th rowspan="2">Combination</th><th colspan="3">Average error</th><th colspan="3">Maximum error</th></tr><tr><th>Φ</th><th>F<sub>x</sub></th><th>F<sub>y</sub></th><th>Φ</th><th>F<sub>x</sub></th><th>F<sub>y</sub></th></tr></thead><tbody><tr><td>+50 per cent</td><td>2 per cent</td><td>3 per cent</td><td>5 per cent</td><td>6 per cent</td><td>23 per cent</td><td>22 per cent</td></tr><tr><td>+25 per cent</td><td>1 per cent</td><td>4 per cent</td><td>4 per cent</td><td>6 per cent</td><td>19 per cent</td><td>18 per cent</td></tr><tr><td>-25 per cent</td><td>1 per cent</td><td>6 per cent</td><td>6 per cent</td><td>5 per cent</td><td>29 per cent</td><td>22 per cent</td></tr><tr><td>-50 per cent</td><td>2 per cent</td><td>7 per cent</td><td>8 per cent</td><td>6 per cent</td><td>35 per cent</td><td>26 per cent</td></tr><tr><td>no peanut</td><td>3 per cent</td><td>8 per cent</td><td>10 per cent</td><td>7 per cent</td><td>37 per cent</td><td>28 per cent</td></tr></tbody></table>

*Note.* Average and maximum errors of the potential and forces for setups with different combinations of errors within the area enclosed by the outer isophote of the bar. The last row gives the error induced by not modelling a B/P bulge at all.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** (a) Characteristic diagram for models with different B/P setups. See Fig. 7 and the text in Section 4.2 for more details on the interpretation of the characteristic diagram. The dotted blue line gives the ZVC; the characteristic diagram for the model with 0, 50 and 100 per cent the fiducial peanut strength is given by the dashed red line, the magenta dash-dotted line and the solid black line, respectively. (b) The same 3/1 orbit, which cuts the y-axis at y = 1, in three models: without a B/P bulge (dashed red line), with 50 per cent peanut strength (magenta dash-dotted line) and with 100 per cent peanut strength (black solid line).</figcaption></figure>

diagram of the most relevant families of periodic orbits for three models with different peanut strengths (100, 50 and 0 per cent of the fiducial strength), shown in Fig. 9(a).

For the model with 50 per cent the fiducial strength the extent of the $x_2$ orbits is reduced by ~19 per cent, while if we do not add a B/P at all, the extent of the $x_2$ family is reduced by ~43 per cent (more than double the error for the 50 per cent peanut strength case).

The bifurcation locus of the 3/1 family for the model with 50 per cent the peanut strength occurs about halfway between the locus of the models with and without a B/P bulge. Additionally, the extent

of the 3/1 family in the characteristic diagram for this model is almost the same as for the model with the fiducial peanut, while the extent of the 3/1 family without a B/P is significantly shortened. In Fig. 9(b), we can see how the 3/1 orbits are affected by the incorrect modelling of the B/P bulge: for the same cut along the y-axis the 3/1 orbits are more elongated in the case with the fiducial peanut model, while they become less elongated and more concave with respect to the bar as the peanut strength is reduced. However, as expected, the orbits in the model with 50 per cent the peanut bulge better match the orbit of the fiducial B/P setup than the model without the B/P bulge.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Top row – left: values of $z_0$ for the fiducial peanut strength (solid thick black line), for 50 per cent larger (solid thin red line), 25 per cent larger (dashed green line), 25 per cent less (dotted magenta line), 50 per cent less (dot-dashed cyan line) and 0 per cent peanut strength, i.e. an isothermal sheet (thin solid black line). Right: bar-induced torque $Q_T$ as a function of radius for models with the aforementioned height functions (respective colours). Middle row – left: the values of $z_0$ for the fiducial peanut width (solid thick black line), 50 per cent larger (solid thin red line), 25 per cent larger (dashed green line), 25 per cent less (dotted magenta line) and 50 per cent less peanut width (dot-dashed cyan line). Right: bar-induced torque $Q_T$, for aforementioned models (respective colours). Bottom row – left: the values of $z_0$ for models with the fiducial peanut length (thick solid black line), 30 per cent longer peanut (dashed red line), 30 per cent shorter peanut (dash-dotted green line) and a model without a peanut (i.e. 0 per cent peanut strength-thin solid black line). Right: bar-induced torque $Q_T$ for the aforementioned models (respective colours). In all plots, the vertical lines correspond to the positions of the peanut maxima for each respective height function.</figcaption></figure>

As already noted in Section 4.2, the $x_1$ family of orbits is also affected by the presence of the B/P bulge, when the maximum extent of the orbits reach the region where the effect of the B/P is maximum, i.e. around $(x, y) = (\pm 3 \text{ kpc}, 0 \text{ kpc})$. On the characteristic

diagram this occurs in the region around $(E_1, y_0) = (0.9, 0.5 \text{ kpc})$. However, even by underestimating the strength of the B/P by 50 per cent, the $x_1$ family is quite similar to the $x_1$ family of the fiducial B/P case.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Left: values of $z_0$ for different combinations of uncertainties. Right: strength of non-axisymmetric forcings $Q_T$ as a function of radius, for the different combinations of uncertainties. In order to not overclutter the plot, the positions of the peanut maxima are given by the vertical arrows, from left to right, for the -50 per cent, -25 per cent, fiducial, +25 per cent and +50 per cent cases.</figcaption></figure>

We see that in general, the characteristic diagram of the model with 50 per cent the fiducial strength has features which are more towards the fiducial peanut model and therefore even with such a large error in peanut strength, the characteristic diagram of this model reproduces relatively well the characteristic diagram of the fiducial B/P model and certainly better than the model without a B/P bulge. Similar results are found when considering errors in peanut width and length, and we therefore conclude that it is preferable to include a B/P in the model; the orbital structure of the model is significantly affected when a B/P bulge is present, and by adding a B/P, even with large errors in its parameters, the periodic orbits reproduce the correct structure more closely than when not including a B/P at all.

## 5.3 Bar strength

In this section, we examine how both the relative errors of $Q_T$ and those of its maximum value $Q_b$ will be affected by uncertainties in the different parameters of the peanut model. We also introduce a new measure of bar strength, $Q_T^{int}$, which takes into consideration the integrated bar-induced torque, along the entire range of the bar. We do so because even though $Q_b$ remains relatively unchanged when adding a B/P bulge to the model, $Q_T$ over its whole range is significantly affected (see for example Fig. 10f and Table 1), and we wish to have a measure of this change with a single number. The bar strength as defined by $Q_T^{int}$ is given by

$$Q_T^{int} = \frac{1}{r_{\text{disc}}} \int_0^{r_{\text{bar}}} Q_T \, dr, \qquad (18)$$

where $r_{\text{disc}}$ is the disc scalelength.

To get a good estimate of the difference of $Q_T$ between two models over the entire radial range, it is best to carry out a point by point comparison, and then consider the radially averaged relative $Q_T$ error. The relative error of $Q_T^{int}$ is a better proxy for this error than $Q_b$, although there are cases where the relative error of $Q_T^{int}$ is small, while the average relative error of $Q_T$ is much more significant (such as the first row of Table 8) or vice versa (first row of Table 9). Therefore, it is possible to have two cases with identical $Q_T^{int}$, but locally different $Q_T$.

### 5.3.1 Peanut strength uncertainties

We see in Figs 10(a) and (b) that $Q_T$ increases as we reduce the strength of the peanut, and it reaches its maximum value when the

<table><caption>**Table 6.** Percentage error of bar strength due to peanut strength uncertainty.</caption><thead><tr><th>Peanut strength</th><th>⟨Error Q<sub>T</sub>⟩</th><th>MAX(Error Q<sub>T</sub>)</th><th>Q<sub>b</sub></th><th>Q<sub>T</sub><sup>int</sup></th></tr></thead><tbody><tr><td>+50 per cent</td><td>8 per cent</td><td>17 per cent</td><td>1.4 per cent</td><td>6 per cent</td></tr><tr><td>+25 per cent</td><td>4 per cent</td><td>9 per cent</td><td>0.7 per cent</td><td>3 per cent</td></tr><tr><td>-25 per cent</td><td>5 per cent</td><td>11 per cent</td><td>0.5 per cent</td><td>4 per cent</td></tr><tr><td>-50 per cent</td><td>11 per cent</td><td>26 per cent</td><td>1.6 per cent</td><td>7 per cent</td></tr><tr><td>no peanut</td><td>27 per cent</td><td>74 per cent</td><td>4 per cent</td><td>20 per cent</td></tr></tbody></table>

Note. The error induced in the bar strength due to different amount of error in the peanut strength. We see the effect of these uncertainties on the average and maximum relative error in $Q_T$ (⟨Error $Q_T$⟩ and MAX(Error $Q_T$), respectively), as well as on the relative errors of $Q_b$ and $Q_T^{int}$.

peanut strength is zero, which corresponds to the height function of a flat isothermal sheet.

The values of the average and maximum relative errors of $Q_T$ (⟨Error $Q_T$⟩ and MAX(Error $Q_T$), respectively), as well as the relative error of $Q_b$ and $Q_T^{int}$ can be seen in Table 6 (and in all subsequent tables in the following subsections). We see that when we compare an isothermal sheet to the fiducial peanut model the error in $Q_b$ is of the order of 4 per cent. This is not representative of the large change that the average relative error of $Q_T$ undergoes (27 per cent). This is due to the fact that the maximum of $Q_T$ does not change much, even though $Q_T$ itself is affected by a significant amount over its entire range (see Fig. 10b). On the other hand, the change in $Q_T^{int}$, which takes into account the whole bar region, is more representative of the change in the average relative error of $Q_T$ (20 per cent).

In all the cases and for all the measurements of bar strength, the error introduced in the model due to uncertainty in peanut strength is not as large as the error introduced when not including a B/P bulge in the model.

### 5.3.2 Peanut width uncertainties

We compare setups with varying peanut widths to our fiducial model. Comparisons for $Q_T$ can be seen in Figs 10(c) and (d) and as previously mentioned, the mismatch between the different models is found when the scaleheights of the models are different. The extent of the region of $Q_T$ which is flattened is reduced when the width of the B/P is reduced, as expected. Conversely, when we increase the width of the B/P bulge, the area of $Q_T$ which is flattened is increased.
<table><caption>**Table 7.** Percentage error of bar strength due to peanut width uncertainty.</caption><thead><tr><th>Peanut width</th><th>⟨Error Q<sub>T</sub>⟩</th><th>MAX(Error Q<sub>T</sub>)</th><th>Q<sub>b</sub></th><th>Q<sub>T</sub><sup>int</sup></th></tr></thead><tbody><tr><td>+50 per cent</td><td>6 per cent</td><td>12 per cent</td><td>9 per cent</td><td>6 per cent</td></tr><tr><td>+25 per cent</td><td>4 per cent</td><td>11 per cent</td><td>5 per cent</td><td>4 per cent</td></tr><tr><td>-25 per cent</td><td>7 per cent</td><td>20 per cent</td><td>3 per cent</td><td>6 per cent</td></tr><tr><td>-50 per cent</td><td>16 per cent</td><td>39 per cent</td><td>5 per cent</td><td>13 per cent</td></tr><tr><td>no peanut</td><td>27 per cent</td><td>74 per cent</td><td>4 per cent</td><td>20 per cent</td></tr></tbody></table>

Note. As in Table 6 but for errors in peanut width.

<table><caption>**Table 8.** Percentage error of bar strength due to peanut length uncertainty.</caption><thead><tr><th>Peanut length</th><th>⟨Error Q<sub>T</sub>⟩</th><th>MAX(Error Q<sub>T</sub>)</th><th>Q<sub>b</sub></th><th>Q<sub>T</sub><sup>int</sup></th></tr></thead><tbody><tr><td>+30 per cent</td><td>16 per cent</td><td>35 per cent</td><td>17 per cent</td><td>5 per cent</td></tr><tr><td>+16 per cent</td><td>9 per cent</td><td>19 per cent</td><td>7 per cent</td><td>3 per cent</td></tr><tr><td>+8 per cent</td><td>5 per cent</td><td>10 per cent</td><td>3 per cent</td><td>1 per cent</td></tr><tr><td>-8 per cent</td><td>4 per cent</td><td>10 per cent</td><td>2 per cent</td><td>1 per cent</td></tr><tr><td>-16 per cent</td><td>9 per cent</td><td>19 per cent</td><td>4 per cent</td><td>3 per cent</td></tr><tr><td>-30 per cent</td><td>16 per cent</td><td>39 per cent</td><td>8 per cent</td><td>8 per cent</td></tr><tr><td>no peanut</td><td>27 per cent</td><td>74 per cent</td><td>4 per cent</td><td>20 per cent</td></tr></tbody></table>

Note. As in Table 6 but for errors in peanut length.

Values for the errors in the bar region are given in Table 7. We see that errors in peanut width do not induce very large errors in the average relative error of $Q_T$, compared to the errors induced when not including a B/P bulge. The errors induced in $Q_T^{int}$ are not very large either, although $Q_b$, in the case of +50 per cent peanut width, has a relative error larger than that of not including a B/P bulge. This again shows the importance of carrying out a point by point comparison, and a comparison of $Q_T^{int}$, in order to determine the errors induced in bar strength due to uncertainties.

### 5.3.3 Peanut length uncertainties

The results of this study are shown in Figs 10(e) and (f) and in Table 8. Something worth noting in the Fig. 10(f) is that the flattening of $Q_T$ occurs at the positions where the maxima of the peanut are found (which are indicated by the corresponding vertical lines).

In Table 8, we see the errors induced in the different measurements of bar strength due to uncertainties in peanut length. For the case where the peanut length is 30 per cent larger than the fiducial value, $Q_b$ has a relative error of 17 per cent compared to the error of 4 per cent in $Q_b$ when we do not add a peanut. This seems to suggest that it can be counter-productive to include a B/P in the model when there exist uncertainties in peanut length. However, if we examine the average relative error in $Q_T$, we see that the error induced in $Q_T$ is in fact larger when we do not model a B/P than when we miscalculate its length by +30 per cent. This points once again to the need for examining the average errors of $Q_T$ and not just $Q_b$, as the errors induced in $Q_b$ are not representative of the error induced in $Q_T$.

We also see in Fig. 10(f), that even though $Q_T$ of the two cases (of fiducial peanut and +30 per cent peanut length) differs significantly point by point, the area under the curve for the two cases is quite similar. This is reflected in the value of the relative error of $Q_T^{int}$, which only suffers a change of around 5 per cent while the average error of $Q_T$ suffers a change of 16 per cent. We see therefore that $Q_T^{int}$ is not always a good approximation for the average relative error of $Q_T$.

<table><caption>**Table 9.** Percentage error of bar strength due to a combination of uncertainties.</caption><thead><tr><th>All Errors</th><th>⟨Error Q<sub>T</sub>⟩</th><th>MAX(Error Q<sub>T</sub>)</th><th>Q<sub>b</sub></th><th>Q<sub>T</sub><sup>int</sup></th></tr></thead><tbody><tr><td>+50 per cent</td><td>19 per cent</td><td>42 per cent</td><td>47 per cent</td><td>25 per cent</td></tr><tr><td>+25 per cent</td><td>16 per cent</td><td>35 per cent</td><td>28 per cent</td><td>14 per cent</td></tr><tr><td>-25 per cent</td><td>16 per cent</td><td>50 per cent</td><td>7 per cent</td><td>11 per cent</td></tr><tr><td>-50 per cent</td><td>20 per cent</td><td>63 per cent</td><td>7 per cent</td><td>17 per cent</td></tr><tr><td>no peanut</td><td>27 per cent</td><td>74 per cent</td><td>4 per cent</td><td>20 per cent</td></tr></tbody></table>

Note. As in Table 6 but for different combinations of uncertainties.

The important thing to note is that all the cases considered induce less error in the average error of $Q_T$ than not modelling the B/P at all.

### 5.3.4 Combinations of uncertainties

As has already been discussed, the most likely scenario is that of a combination of different sources of error affecting our model. The combinations of errors shown in Table 9 and Fig. 11, are as in Section 5.1.4. We see that the average and maximum relative error in $Q_T$ for all the combinations is less than that of not modelling the B/P at all. The scaleheights and bar strength for these models can be seen in Figs 10(a) and (b), respectively.

## 6 SUMMARY AND CONCLUSIONS

In this paper, we present the effects of a B/P height function on the potential, forces, periodic orbits and bar strength of a barred galaxy. We show that such height functions significantly affect the results, which consequently hints to the effects that a B/P bulge will have on its host galaxy.

We present a method for calculating the potential and forces due to the stellar component of a disc galaxy, based on a three-dimensional integration of the stellar density distribution, which can be obtained from images of not too inclined galaxies combined with a given height function. The method gives robust results for different test cases, as well as allowing for any general height function to be used, thus allowing for complex density distributions to be modelled.

We used our code on an image extracted from a N-body+SPH simulation of an isolated galaxy, together with two flat, position-independent height functions, and two position-dependent height functions. Of the two position-dependent height functions, one models a peanut bulge and one models a boxy bulge. To create an accurate and physically motivated fiducial height function for the peanut, we shaped and fitted our peanut height function to the B/P bulge of the simulated galaxy.

We found, in accordance with previous results in the literature (Laurikainen & Salo 2002), that for the two flat height functions the potential and forces do not vary much, provided the setups have equivalent scaleheights. This also holds true for the bar strength $Q_T$, which does not change much for different flat height functions.

However, we found that for boxy or peanut height functions the potential and forces vary significantly with respect to the case in which a flat height function is used (see Fig. 5). For the potential, the difference can be up to 7 per cent for an extended region within the bar. For $F_x$ the difference can be as large as 37 per cent, while for $F_y$ this difference can be as large as 28 per cent. We therefore concluded that if a B/P bulge is present, one should include it when creating a dynamical model of the galaxy.
To further confirm this result, we examined the effect of the B/P bulge on the morphology of the most important families of periodic orbits found in barred galaxies. We see that by taking into account the B/P geometry (i.e. by using our fiducial peanut height function) for a given energy, the elongation of the $x_1$ orbits – the bar-supporting orbits elongated parallel to the bar – is decreased; this effect is most noticeable for orbits in the region where the B/P is maximum (around $\pm 3$ kpc, see Fig. 7a) as expected. By adding a B/Pt to the model the extent of the $x_2$ family – the family of orbits perpendicular to the bar – is increased by ~43 per cent (see Fig. 6b), as is the distance between the two ILRs. Additionally, the position of the 3:1 resonance is changed; the 3/1 family – elongated along the bar and asymmetric with respect to the y-axis – appears at larger energies and is much more extended in the characteristic diagram (see Fig. 6b). All the aforementioned effects will have an impact on the stellar as well as the gaseous kinematics of the galaxy. The shape and strength of the shocks in the gas will be affected, which in turn affects the amount of gas inflow to the central parts of the galaxy. This could have an impact on the formation of discy bulges and possibly on the fuel reservoir for AGN activity. We plan to investigate in future work the extent of the effects of B/P bulges on gas flows in galaxies.

We also studied the effect of the B/P bulge on the bar strength, as given by the non-axisymmetric forcings due to the bar, $Q_T$. The shape as well as the maximum of $Q_T$ are significantly affected by taking into account the geometry of a B/P bulge. We found it useful to define a new quantity for measuring bar strength, $Q_{\text{bar}}^n$, which allows us to extract information about the strength of the bar by using its whole extent. The presence of a Boxy/Peanut bulge, especially at the points where its scaleheight is maximum, reduces the bar strength (see Fig. 8) which confirms that the presence of a B/P bulge reduces the bar induced torques.

Even though taking into account the geometry of B/P bulges will affect the model, it is not trivial to obtain their parameters observationally. We therefore examined how much error would be introduced in the results by introducing uncertainties in the B/P parameters. Each source of error individually (peanut strength, peanut length and peanut width), as well as combinations of the different sources of error, induce errors in the results which in general are considerably less than those induced by not modelling the peanut at all. So, for realistic values of uncertainties in the peanut parameters, the error in including a peanut will be less than the error induced by not including a peanut in the model.

The simulated galaxy, we chose for this study contains a strong bar, corresponding to bar classes 5 and 6 from the Buta & Block (2001) classification. Therefore, the results of this study can be straightforwardly and quantitatively applied to real galaxies with similar bar and peanut strength, which account for approximately 20 per cent of SB galaxies in the local Universe. Our results are also qualitatively relevant to all barred galaxies in the secular evolution phase, although for reduced bar and peanut strength the effect of the B/P bulge on the model is also reduced. In this work, we have presented an in depth study of the effects of a B/P bulge on its galaxy model, focusing on a particular test case; we plan to present a quantitative statistical study of the effects of these bulges on their host galaxies elsewhere.

## ACKNOWLEDGEMENTS

We acknowledge financial support to the DAGAL 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. EA and AB also acknowledge financial support from the CNES (Centre National d'Etudes Spatiales - France). The simulation analysed here was made using HPC resources from GENCI-TGCC/CINES (Grants 2013 - x2013047098 and 2014 - x2014047098). We would also like to thank Dimitri Gadotti for useful discussions on the properties of barred galaxies and for helpful suggestions on the manuscript and Sergey Rodionov for software assistance.

## REFERENCES

* Aguerri J. A. L., Méndez-Abreu J., Corsini E. M., 2009, A&A, 495, 491

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

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

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

* Athanassoula E., 2005, MNRAS, 358, 1477

* Athanassoula E., 2008a, in Wada K., Combes F., eds, Mapping the Galaxy and Nearby Galaxies, Astrophys Space Sci. Proc., Springer-Verlag, New York, p. 47

* Athanassoula E., 2008b, in Bureau M., Athanassoula E., Barbuy B., eds, Proc. IAU Symp. 245, Formation and Evolution of Galactic Bulges, Cambridge Univ. Press, Cambridge, p. 93

* Athanassoula E., 2015, in Laurikainen E., Peletier R., Gadotti D., eds, Galactic Bulges, Boxy/peanut/X bulges, barlenless and the thick part of galactic bars: What are they and how did they form?, Springer, preprint (arXiv:1503.04804)

* Athanassoula E., Beaton R. L., 2006, MNRAS, 370, 1499

* Athanassoula E., Bureau M., 1999, ApJ, 522, 699

* Athanassoula E., Misiriotis A., 2002, MNRAS, 330, 35

* Athanassoula E., Fady E., Lambert J. C., Bosma A., 2000, MNRAS, 314, 475

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

* Athanassoula E., Laurikainen E., Salo H., Bosma A., 2014, preprint (arXiv:1405.6726)

* Barazza F. D., Jogee S., Marinova I., 2008, ApJ, 675, 1194

* Binney J., 1981, MNRAS, 196, 455

* Binney J., Tremaine S., 2008, Galactic Dynamics, 2nd edn. Princeton Univ. Press, Princeton, NJ

* Bureau M., Athanassoula E., 1999, ApJ, 522, 686

* Bureau M., Athanassoula E., 2005, ApJ, 626, 159

* Buta R., Block D. L., 2001, ApJ, 550, 243

* Chung A., Bureau M., 2004, AJ, 127, 3192

* Coelho P., Gadotti D. A., 2011, ApJ, 743, L13

* Combes F., 2001, in Aretxaga I., Kunth D., Mújica R., eds, Advanced Lectures on the Starburst-AGN Connection. World Scientific, Singapore, p. 223

* Combes F., Sanders R. H., 1981, A&A, 96, 164

* Combes F., Debbasch F., Friedli D., Pfenniger D., 1990, A&A, 233, 82

* Contopoulos G., Grosbol P., 1989, A&AR, 1, 261

* Contopoulos G., Harsoula M., 2013, MNRAS, 436, 1201

* Debattista V. P., Carollo C. M., Mayer L., Moore B., 2005, ApJ, 628, 678

* Dehnen W., 1993, MNRAS, 265, 250

* Emsellem E., Renaud F., Bournaud F., Elmegreen B., Combes F., Gabor J. M., 2015, MNRAS, 446, 2468

* Erwin P., Debattista V. P., 2013, MNRAS, 431, 3060

* Eskridge P. B. et al., 2000, AJ, 119, 536

* Foyle K., Rix H.-W., Zibetti S., 2010, MNRAS, 407, 163

* Gadotti D. A., 2009, MNRAS, 393, 1531

* Gadotti D. A., 2011, MNRAS, 415, 3308

* Galassi M., Davies J., Theiler J., Gough B., Jungman G., Booth M., Rossi F., 2003, Gnu Scientific Library: Reference Manual. Network Theory Ltd

* Gardner E., Debattista V. P., Robin A. C., Vásquez S., Zoccali M., 2014, MNRAS, 438, 3275

* Haan S., Schinnerer E., Emsellem E., García-Burillo S., Combes F., Mundell C. G., Rix H.-W., 2009, ApJ, 692, 1623

* Harsoula M., Kalapotharakos C., 2009, MNRAS, 394, 1605
Hockney R. W., Eastwood J. W., 1988, Computer simulation using particles. Hilger, Bristol

Howard C. D. et al., 2009, ApJ, 702, L153

Kalapotharakos C., Patsis P. A., Grosbøl P., 2010a, MNRAS, 403, 83

Kalapotharakos C., Patsis P. A., Grosbøl P., 2010b, MNRAS, 408, 9

Kormendy J., 2015, in Laurikainen E., Peletier R., Gadotti D. A., eds, Galactic Bulges, Section 7. Springer Int. Publ. AG

Kormendy J., Kennicutt R. C., Jr, 2004, ARA&A, 42, 603

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

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

Laurikainen E., Salo H., 2002, MNRAS, 337, 1118

Laurikainen E., Salo H., Buta R., Knapen J. H., 2011, MNRAS, 418, 1452

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

Lee G.-H., Woo J.-H., Lee M. G., Hwang H. S., Lee J. C., Sohn J., Lee J. H., 2012, ApJ, 750, 141

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

Lin L.-H., Taam R. E., Yen D. C. C., Muller S., Lim J., 2011, ApJ, 731, 15

Lin L.-H., Wang H.-H., Hsieh P.-Y., Taam R. E., Yang C.-C., Yen D. C. C., 2013, ApJ, 771, 8

Lütticke R., Dettmar R.-J., Pohlen M., 2000, A&AS, 145, 405

McWilliam A., Zoccali M., 2010, ApJ, 724, 1491

Martinez-Valpuesta I., Shlosman I., Heller C., 2006, ApJ, 637, 214

Méndez-Abreu J., Corsini E. M., Debattista V. P., De Rijcke S., Aguerri J. A. L., Pizzella A., 2008, ApJ, 679, L73

Menéndez-Delmestre K., Sheth K., Schinnerer E., Jarrett T. H., Scoville N. Z., 2007, ApJ, 657, 790

Merritt D., 1996, AJ, 111, 2462

Mihos J. C., Walker I. R., Hernquist L., Mendes de Oliveira C., Bolte M., 1995, ApJ, 447, L87

Minchev I., Famaey B., 2010, ApJ, 722, 112

Miyamoto M., Nagai R., 1975, PASJ, 27, 533

Nataf D. M. et al., 2013, ApJ, 769, 88

Nataf D. M., Cassisi S., Athanassoula E., 2014, MNRAS, 442, 2075

Nataf D. M. et al., 2015, MNRAS, 447, 1535

Ness M. et al., 2012, ApJ, 756, 22

Ness M. et al., 2013a, MNRAS, 430, 836

Ness M. et al., 2013b, MNRAS, 432, 2092

Patsis P. A., Athanassoula E., Quillen A. C., 1997, ApJ, 483, 731

Patsis P. A., Skokos C., Athanassoula E., 2002, MNRAS, 337, 578

Patsis P. A., Skokos C., Athanassoula E., 2003, MNRAS, 342, 69

Patsis P. A., Kalapotharakos C., Grosbøl P., 2010, MNRAS, 408, 22

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

Pfenniger D., 1984, A&A, 134, 373

Pfenniger D., 1985, A&A, 150, 112

Querejeta et al., 2014, preprint (arXiv:1410.0009)

Quillen A. C., Frogel J. A., Gonzalez R. A., 1994, ApJ, 437, 162

Quinn T., Katz N., Stadel J., Lake G., 1997, preprint (arXiv:e-prints)

Raha N., Sellwood J. A., James R. A., Kahn F. D., 1991, Nature, 352, 411

Roškar R., Debattista V. P., Stinson G. S., Quinn T. R., Kaufmann T., Wadsley J., 2008, ApJ, 675, L65

Sellwood J. A., Binney J. J., 2002, MNRAS, 336, 785

Shen J., Rich R. M., Kormendy J., Howard C. D., De Propris R., Kunder A., 2010, ApJ, 720, L72

Shlosman I., Begelman M. C., Frank J., 1990, Nature, 345, 679

Skokos C., Patsis P. A., Athanassoula E., 2002a, MNRAS, 333, 847

Skokos C., Patsis P. A., Athanassoula E., 2002b, MNRAS, 333, 861

Slyz A. D., Kranz T., Rix H.-W., 2003, MNRAS, 346, 1162

Springel V., 2005, MNRAS, 364, 1105

van der Kruit P. C., 1988, A&A, 192, 117

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

Vásquez S. et al., 2013, A&A, 555, A91

Wegg C., Gerhard O., 2013, MNRAS, 435, 1874

Weiland J. L. et al., 1994, ApJ, 425, L81

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

Zaritsky D., Lo K. Y., 1986, ApJ, 303, 66

## APPENDIX A: THE ANALYTIC GALAXY MODEL

To model the disc, we use a Miyamoto-Nagai density (Miyamoto & Nagai 1975) which is defined by the potential-density pair,

$$ \Phi_{\text{MN}}(R, z) = - \frac{GM_D}{\sqrt{R^2 + (a + \sqrt{z^2 + b^2})^2}}, \quad (A1) $$

$$ \begin{aligned} \rho_{\text{MN}}(R, z) = & \left( \frac{b^2 M_D}{4\pi} \right) \frac{a R^2 + (a + 3\sqrt{z^2 + b^2})(a + \sqrt{z^2 + b^2})^2}{[R^2 + (a + \sqrt{z^2 + b^2})^2]^{5/2}(z^2 + b^2)^{3/2}}, \quad (A2) \end{aligned} $$

where $R$ and $z$ are the cylindrical coordinates, $G$ is Newton's gravitational constant, $M$ the total mass of the system and $a$ and $b$ are its characteristic lengths. We set the parameters $a$ and $b$ to 9 and 1.8, respectively, such that we obtain a realistic exponential disc with a scalelength of about 3 kpc with its mass set to 0.56 times the total mass of the system (Gadotti 2011).

The bulge is modelled using a Dehnen sphere (Dehnen 1993), where we set $\gamma = 0.5$ in order to obtain a cuspy density distribution. The potential density pair is given by

$$ \rho(r) = \frac{5}{8\pi} \frac{r_B M_B}{\sqrt{r(r_B + r)^{7/2}}}, \quad (A3) $$

and

$$ \Phi(r) = \frac{-2GM_B}{3r_B} \left( 1 - \left( \frac{r}{r + r_B} \right)^{3/2} \right), \quad (A4) $$

where $r_B$ is a characteristic radius of the system. The mass of the bulge is set to 0.34 the total mass of the model which is a typical value for the bulge mass (Gadotti 2011).

The bar is modelled using a Ferrers ellipsoid (Ferrers 1877), whose density is given by

$$ \rho = \begin{cases} \rho_0(1 - m^2)^\mu & m \le 1 \\ 0 & m \ge 1 \end{cases}, \quad (A5) $$

where $m^2$ is

$$ m^2 = \frac{x^2}{\alpha^2} + \frac{y^2}{\beta^2} + \frac{z^2}{\gamma^2}. \quad (A6) $$

The central density of the bar is given by $\rho_0$, while $n$ sets the decrease in bar density as a function of position and $\alpha$, $\beta$ and $\gamma$ give the sizes of the three semiprincipal axes. The mass of the bar is 0.1 times the total mass, which is again a typical value for real galaxies (Gadotti 2011). The bar's semimajor axis is set to $a = 5$ kpc, with an axial ratio of $a/b = 2.5$, and we use the inhomogeneous $n = 1$ case.

When integrating orbits in these potentials, we do so in the rotating frame of reference of the bar, where the bar potential rotates with a pattern speed which is set such that corotation occurs just outside the end of the bar within the range $1.4 > R_{CR}/R_{bar} > 1$, where $R_{CR}$ and $R_{bar}$ are the corotation and bar radius, respectively (e.g. Athanassoula 1992b). The equations of motion in a rotating frame of reference are taken from chapter 3 of Binney & Tremaine (2008), where the fictitious forces due to rotation are taken into account.

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