# Deprojection of external barred galaxies from photometry

Behzad Tahmasebzadeh,1,2 ID Ling Zhu,1 ID Juntai Shen,3,4,1 ID Ortwin
Gerhard,5 ID Yujing Qin6 ID

1 Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan
Road, Shanghai 200030, China
2 Department of Astronomy and Space Sciences, University of Chinese Academy of
Sciences, 19A Yuquan Road, Beijing 100049, China
3 Department of Astronomy, School of Physics and Astronomy, Shanghai Jiao Tong
University, 800 Dongchuan Road, Shanghai 200240, China
4 Key Laboratory for Particle Astrophysics and Cosmology (MOE) / Shanghai Key
Laboratory for Particle Physics and Cosmology, Shanghai 200240, China
5 Max-Planck-Institut für Extraterrestrische Physik, Gießenbachstraße 1, 85748
Garching, Germany
6 Department of Astronomy/Steward Observatory, 933 North Cherry Avenue,
University of Arizona, Tucson, AZ 85721, USA behzadtahmaseb@gmail.com,Corr
author: lzhu<EMAIL_ADDRESS>
###### Abstract

The observations of external galaxies are projected to the 2D sky plane.
Reconstructing the 3D intrinsic density distribution of a galaxy from the 2D
image is challenging, especially for barred galaxies, but is a critical step
for constructing galactic dynamical models. Here we present a method for
deprojecting barred galaxies and we validate the method by testing against
mock images created from an N-body simulation with a peanut-shaped bar. We
decompose a galaxy image into a bulge (including a bar) and a disk. By
subtracting the disk from the original image a barred bulge remains. We
perform multi-Gaussian expansion (MGE) fit to each component, then we
deproject them separately by considering the barred bulge is triaxial while
the disk is axisymmetric. We restrict the barred bulge to be aligned in the
disk plane and has a similar thickness to the disk in the outer regions. The
3D density distribution is thus constructed by combining the barred bulge and
the disk. Our model can generally recover the 3D density distribution of disk
and inner barred bulge regions, although not a perfect match to the peanut-
shaped structure. By using the same initial conditions, we integrate the
orbits in our model-inferred potential and the true potential by freezing the
N-body simulation. We find that $85\%$ of all these orbits have similar
morphologies in these two potentials, and our model supports the orbits that
generate a boxy/peanut-shaped structure and an elongated bar similar to these
in the true potential.

###### keywords:

galaxies: bulges - galaxies: fundamental parameters - galaxies: photometry -
galaxies: structure.

††pubyear: 2021††pagerange: Deprojection of external barred galaxies from
photometry–15

##  1 Introduction

A large fraction of disk galaxies ($30\%$ to $65\%$) have bars with various
strengths in the central regions Eskridge et al. (2000); Menéndez-Delmestre et
al. (2007); Barazza et al. (2008); Aguerri et al. (2009); Gadotti (2009);
Erwin (2018). In the case of face-on and moderately inclined galaxies, bars
appear as non-axisymmetric perturbations in the surface density maps. For the
edge-on or highly inclined galaxies, bars are detectable by particular
kinematic signatures, e.g., a positive correlation between mean velocity and
the third Gauss-Hermite moment $h_{3}$ Bureau & Athanassoula (2005); Li et al.
(2018). Bars can redistribute the angular momentum and energy of the disk
material, so they drive the morphological evolution of disk galaxies Weinberg
(1985); Debattista & Sellwood (1998); Athanassoula (2003); Kormendy &
Kennicutt (2004). They are well known to play a major role in the secular
galaxy evolution Friedli & Benz (1993); Sheth et al. (2005); Masters et al.
(2011); Gadotti (2011).

A bar has been included in the dynamical models for Milky Way by the
Schwarzschild orbit-superposition method Zhao (1996); Wang et al. (2013) and
the made-to-measure (M2M) method Long et al. (2013); Hunt et al. (2013);
Portail et al. (2015, 2017). There are several commonly used implementations
of the Schwarzschild’s orbit-superposition method for external galaxies by
taking different assumptions of their geometries: they could be spherical
Richstone & Tremaine (1985); Breddels et al. (2013); Kowalczyk et al. (2017),
axisymmetric Cretton et al. (1999); Gebhardt et al. (2000); Valluri et al.
(2004) or triaxial Cappellari et al. (2006); Van den Bosch et al. (2008); Zhu
et al. (2018b). However, bar shapes and the figure rotation are not included
explicitly in these models.

A major goal of dynamical modeling is to obtain the underlying mass
distribution, including the central supermassive black hole mass, the stellar
mass profile, and the dark matter halo profile. However, without including the
bar properly, these results could be significantly biased Brown et al. (2013).
Dynamical modeling of external barred galaxies are still in the early stage.
Blaña Díaz et al. (2018) made a triaxial bulge/bar/disk M2M model for M31
taking an N-body model which generally matches the bulge properties of M31 as
an initial condition of the M2M algorithm. A bar has been included in the
recently developed Schwarzschild code, SMILE Vasiliev (2013) and FORSTAND
Vasiliev & Valluri (2020), which are, however, only applied to mock data
created from a simulation by using its real 3D density distribution.
Estimating the 3D density distribution of a real barred galaxy is a key step
still missing before we can create proper dynamical models to a real barred
galaxy from observations in a general sense.

It is non-trivial to obtain the intrinsic 3D density distributions from the 2D
images on the observational plane for triaxial systems, and even for
axisymmetric systems when not observed edge-on Rybicki (1987). Several
approaches are developed to make the 3D deprojection feasible. One of such
approaches is to consider ellipsoid luminosity density profiles for stellar
systems with ellipsoidal radius of
$m=\sqrt{x^{2}+\frac{y^{2}}{p^{2}}+\frac{z^{2}}{q^{2}}}$ which requires the
axial ratios $q\leq p\leq 1$. The parameters could then be derived by fitting
the projected model to the observed image Contopoulos (1956); Stark (1977);
Binney (1985). However, these methods are not able to reproduce isophotal
twists and ellipticity variations along the radius. Non-parametric
deprojections are also attempted Magorrian (1999); Bissantz & Gerhard (2002),
but they are complicated and require considerable time to converge. Multi-
Gaussian expansion (MGE) Bendinelli (1991); Monnet et al. (1992); Emsellem et
al. (1994a); Cappellari (2002) is an efficient method to describe the surface
brightness of a 2D image and deproject it to the 3D luminosity density
distribution. It can reproduce isophotal twists and ellipticity variations
along the radius for a triaxial rigid body. In addition, the gravitational
potential and force components can be computed based on the 3D density
distribution derived from this method Van den Bosch et al. (2008). MGE has
been widely used in dynamical models for different types of galaxies (e.g.,
Emsellem et al., 1994b; Cretton & van den Bosch, 1999; Cappellari, 2008; Zhu
et al., 2018a), but not yet including bars explicitly.

A popular method of deprojecting barred galaxies is to describe their face-on
view by analytical, image stretching, or Fourier-based methods Gadotti et al.
(2007); Noordermeer & van der Hulst (2007); Li et al. (2011), then adopting a
scale height for the vertical profile to construct the 3D density distribution
(e.g., Weiner et al. (2001)). This approach works for galaxies with moderate
inclination angles ( $<\hskip 0.56917pt60^{\circ}$) and the vertical thickness
of the bar is a significant source of uncertainties Zou et al. (2014);
Fragkoudi et al. (2015). The 3D intrinsic shape of bars can also be obtained
with some statistical approaches. Méndez-Abreu et al. (2010); Méndez-Abreu et
al. (2018) developed an approach to derive the intrinsic shape of bulges/bar
based on the geometric information extracted from 2D photometric
decomposition. They assume a bulge/bar is a triaxial ellipsoid that shares the
same equatorial plane as an oblate disk. They use this method to estimate the
intrinsic axes ratios of $83$ barred galaxies from CALIFA survey, and they
find that $68\%$ of bars in their sample are prolate-triaxial ellipsoids and
$32\%$ are oblate-triaxial ellipsoids.

In this paper, we present an efficient way of deprojecting barred galaxies
based on the MGE algorithm, and validate the method by applying it to mock
images created from a simulation. We first create mock images with different
projection angles as described in §2. We introduce our deprojecting approach
for barred galaxies in §3. For verification of our model, in §4 we calculate
the potential, force and analyze the orbital structures according to the 3D
density distribution from our model, and compare with those from the original
simulation. In §5, we summarize and list the main conclusions of our work.

##  2 the mock data

We use an influential N-body simulation of a Milky Way-like galaxy from Shen
et al. (2010). It has a central bar-shaped structure that matches many
observed properties of the Milky Way Qin et al. (2015). The total stellar mass
of the simulation is $4.25\times 10^{10}M_{\odot}$ with $10^{6}$ equal-mass
particles. A rigid pseudo-isothermal DM halo potential is adopted
$\Phi=\frac{1}{2}V_{\mathrm{0}}^{2}\ln\left(1+\frac{r^{2}}{R_{\mathrm{c}}^{2}}\right)$,
in which the scale velocity and scale radius are $V_{\mathrm{0}}=250\hskip
1.42271pt\mathrm{km}\mathrm{s}^{-1}$ and $R_{c}=15\hskip
1.42271pt\mathrm{kpc}$, respectively. The bar has a half-length of $\sim\hskip
0.85355pt4\hskip 1.42271pt\mathrm{kpc}$ and rotates with a pattern speed of
$\Omega_{\mathrm{p}}\sim 39\hskip 2.27626pt\mathrm{km\hskip
1.13791pts^{-1}\hskip 1.13791ptkpc^{-1}}$ (corotation radius $\sim\hskip
0.85355pt4.7\hskip 1.42271pt\mathrm{kpc}$). The end-to-end separation between
the outer two edges of the $\mathrm{X}$-shaped structures is $\sim\hskip
0.85355pt4\hskip 1.42271pt\mathrm{kpc}$ along the major axis and $\sim\hskip
0.85355pt2.4\hskip 1.42271pt\mathrm{kpc}$ along the vertical minor axis Li &
Shen (2012).

Throughout the paper, we use the coordinate $(x,y,z)$ to describe the
intrinsic 3D structure, where $x$, $y$, $z$ are aligned with the long,
intermediate and short axes of the galaxy. While we use the coordinate
$(x^{\prime},y^{\prime})$ to describe the projected structure to the 2D
observational plane, and $z^{\prime}$ is along the line-of-sight. The
orientation of a projection is specified by the viewing angles
$(\theta,\varphi,\psi)$. $\theta$ and $\varphi$ give the orientation of the
line-of-sight with respect to the principal axes of the object. For instance,
projections along the intrinsic major, intermediate, and minor axes correspond
respectively to $(\theta=90^{\circ},\varphi=0^{\circ})$,
$(\theta=90^{\circ},\varphi=90^{\circ})$ and
$(\theta=0^{\circ},\varphi=0^{\circ},...,90^{\circ})$. The angle $\psi$ is
required to determine the rotation of the object around the line of sight (see
Fig. 2 in de Zeeuw & Franx (1989)).

The two coordinate systems are related as Binney (1985):

$\left(\begin{array}[]{l}{x^{\prime}}\\\ {y^{\prime}}\\\
{z^{\prime}}\end{array}\right)=\mathbf{R}\cdot\mathbf{P}\cdot\left(\begin{array}[]{l}{x}\\\
{y}\\\ {z}\end{array}\right),$ (1)

where matrix $\mathrm{P}$ is responsible for the projection onto the sky plane
defined as:

${\mathbf{P}=\left(\begin{array}[]{ccc}{-\sin\varphi}&{\cos\varphi}&{0}\\\
{-\cos\theta\cos\varphi}&{-\cos\theta\sin\varphi}&{\sin\theta}\\\
{\sin\theta\cos\varphi}&{\sin\theta\sin\varphi}&{\cos\theta}\end{array}\right)},$
(2)

and matrix $\mathrm{R}$ expresses the rotation around the line-of-sight by
angle $\psi$:

${\mathbf{R}=\left(\begin{array}[]{ccc}{\sin\psi}&{-\cos\psi}&{0}\\\
{\cos\psi}&{\sin\psi}&{0}\\\ {0}&{0}&{1}\end{array}\right)}.$ (3)

To simplify the description of multiple components of a galaxy, we always
align the major axis of the disk with $x^{\prime}$ axis and minor with
$y^{\prime}$. Note that this could be different from the natural coordinate,
i.e. $y^{\prime}$ aligning with north, when dealing with real observational
data. Here $\psi$ is defined from the $y^{\prime}$ axis, thus we always have
$\psi_{\mathrm{disk}}=90^{\circ}$ for the disk aligned in this way.

To create a mock image, we project the simulation snapshot to the
observational plane with certain viewing angles. We put it at a distance of
$41\hskip 1.42271pt\mathrm{Mpc}$ and then create a surface brightness map with
a spatial resolution of $1\hskip 1.42271pt\mathrm{arcsec}\hskip
2.27626pt\mathrm{pixel}^{-1}$. At this distance, $1\hskip
1.42271pt\mathrm{kpc}$ equals $5\hskip 1.42271pt\mathrm{arcsec}$. We produce
several mock images from the simulation adopting different viewing angles of
$\theta$ and $\varphi$ as listed in Table 1, which are labeled as true values.
To keep the disk major axis aligned with $x^{\prime}$-axis, all of our mock
galaxies have $\psi^{\mathrm{true}}=90^{\circ}$. Each mock image will be taken
as an independent galaxy from observation. In the following sections, we will
illustrate our deprojection model with the mock galaxy $I_{1}$.

Name |  $\theta^{\mathrm{true}}$ $(^{\circ})$ |  $\varphi^{\mathrm{true}}$ (∘)
---|---|---
$I_{1}$ |  60 |  -45
$I_{2}$ |  60 |  -90
$I_{3}$ |  60 |  0
$I_{4}$ |  80 |  -45
$I_{5}$ |  80 |  -90
$I_{6}$ |  80 |  0

Table 1: Mock data sets with different viewing angles of
$\theta^{\mathrm{true}}$ and $\varphi^{\mathrm{true}}$. Inclination angle
$\theta^{\mathrm{true}}$ is the orientation of the disk,
$\theta^{\mathrm{true}}=90^{\circ}$ means edge-on and
$\theta^{\mathrm{true}}=0^{\circ}$ means face-on. $\varphi^{\mathrm{true}}$
describes the orientation of the bar, $\varphi^{\mathrm{true}}=90^{\circ}$
means side-on and $\varphi^{\mathrm{true}}=0^{\circ}$ means end-on. To keep
the disk major axis aligned with $x^{\prime}$-axis, we have
$\psi^{\mathrm{true}}=90^{\circ}$ for all cases.

##  3 Derojection

Here, we give a step-by-step description of our method. We first decompose the
galaxy to a bulge and a disk by using GALFIT Peng et al. (2010). Secondly, we
fit the bulge and the disk using 2D MGE separately. Then, we deproject each
component individually from 2D MGE to 3D MGE. Finally, the deprojected galaxy
is simply the sum of the axisymmetric disk and the triaxial bulge. Meanwhile,
the bulge major axis is restricted to be aligned in the disk plane.

This method allows different intrinsic shapes for the two components. And we
have the freedom to align the bulge in the disk plane with different position
angles.

###  3.1 2D bulge-to-disk decomposition

We use GALFIT 3.0.5
111https://users.obs.carnegiescience.edu/peng/work/galfit/galfit.html to
decompose the surface brightness map of our mock galaxy into 2D elliptical
bulge and disk. Poisson noise is taken as the uncertainty of the image, which
is used to weight the data points in GALFIT fitting.

We use a standard exponential profile to describe the surface brightness of a
disk:

$\Sigma(r^{\prime})=\Sigma_{0}\exp\left(-\frac{r^{\prime}}{r_{s}}\right),$ (4)

where $r_{s}$ is the scale length of disk and $r^{\prime}=\left(x^{\prime
2}+\frac{y^{\prime 2}}{q_{d}^{2}}\right)^{\frac{1}{2}}$,
$(x^{\prime},y^{\prime})$ are aligned with the major and minor axis of the
disk, and $q_{d}$ is the axial ratio.

We use a Sérsic function to describe the surface brightness of an elliptical
bulge:

$\Sigma(r^{\prime})=\Sigma_{e}\exp\left[-\kappa\left(\left(\frac{r^{\prime}}{r_{e}}\right)^{\frac{1}{n}}-1\right)\right],$
(5)

Where $r_{e}$ is the effective radius (or half-light radius), $\Sigma_{e}$ is
the surface brightness at the effective radius $r_{e}$, and $n$ is called the
"Sérsic index", which controls the concentration of the profile. When $n$ is
small, it has a shallow inner profile and a steep truncation at a large
radius. Inversely, When $n$ is large, it has a steep inner profile and a
highly extended outer wing. $\kappa$ is not a free parameter, it coupled to
$n$ by $\gamma\left(2n;\kappa\right)=\frac{1}{2}\Gamma(2n)$, $\Gamma$ and
$\gamma$ are the Gamma function and lower incomplete Gamma function,
respectively. Here $r^{\prime}$ is defined as:

$r^{\prime}=\left(x_{p}^{2}+\frac{y_{p}^{2}}{q_{b}^{2}}\right)^{\frac{1}{2}}.$
(6)

$q_{b}=b/a$, where $b$ and $a$ are the minor and major axis of the elliptical
bulge in the observational plane. $x_{p}$ and $y_{p}$ are aligned with the
major and minor axis of the elliptical bulge:

$\begin{split}&x_{p}=\hskip
4.55254ptx^{\prime}\cos(90-\mathrm{\psi}^{\mathrm{proj}}_{\mathrm{bar}})+y^{\prime}\sin(90-\mathrm{\psi}^{\mathrm{proj}}_{\mathrm{bar}})\\\
&y_{p}=-x^{\prime}\sin(90-\mathrm{\psi}^{\mathrm{proj}}_{\mathrm{bar}})+y^{\prime}\cos(90-\mathrm{\psi}^{\mathrm{proj}}_{\mathrm{bar}})\end{split}$
(7)

where $\mathrm{\psi}^{\mathrm{proj}}_{\mathrm{bar}}$ is the position angle of
the elliptical bulge measured counter-clockwise from the $y^{\prime}$-axis.

Fig. 1 presents the 2D surface brightness profile of (a) the mock galaxy
$I_{1}$, (b) GALFIT exponential disk, (c) barred bulge, which is the residual
of subtracting GALFIT disk from the original image, (d) GALFIT elliptical
bulge. The radial profiles of the surface brightness along the major axis is
presented in Fig. 2. The sum of an exponential disk and an elliptical bulge
gives a reasonable fit to the image. although the model is slightly dimmer in
the regions of $40\lesssim R\lesssim 70$ $\mathrm{arcsec}$, and brighter at
$R\gtrsim 70$ $\mathrm{arcsec}$. Our best-fit parameters are $r_{s}=20.90$
$\mathrm{arcsec}$, $\Sigma_{0}=11.63$ $\mathrm{mag}$ $\mathrm{arcsec}^{-2}$,
$q_{d}=0.51$ for the disk component and $n=1.093$, $r_{e}=5.88$
$\mathrm{arcsec}$, $\Sigma_{e}=11.69$ $\mathrm{mag}$ $\mathrm{arcsec}^{-2}$,
$q_{b}=0.46$, and $\mathrm{\psi}^{\mathrm{proj}}_{\mathrm{bar}}=71^{\circ}$
for the elliptical bulge. Note that we have $\psi_{\mathrm{disk}}=90^{\circ}$,
thus the position angle of the bar is $-19^{\circ}$ different from the disk.

Figure 1: 2D surface brightness profiles of: (a) original mock galaxy $I_{1}$,
(b) GALFIT exponential disk of the best-fit model, (c) barred bulge which is
obtained by subtracting the GALFIT disk from the original image, (d) the
GALFIT elliptical bulge (Sérsic profile). Figure 2:  Surface brightness
profile along the major axis. The blue and green dashed curves are the
exponential disk and the elliptical bulge from GALFIT, while the black solid
curve indicates the sum of them. The red dots are the data of the mock galaxy
$I_{1}$.

###  3.2 2D MGE fit to the bulge and disk

We then separately fit 2D MGEs to (i) the 2D surface decomposed disk and (ii)
the residual barred bulge mock image obtained by subtracting the disk model
from the original mock image. This is done with the $\mathrm{MgeFit}$ software
from Cappellari (2002) 222http://www-astro.physics.ox.ac.uk/~mxc/software/.
The MGE describes the surface brightness (in the unit of
$\mathrm{L_{sun}}\mathrm{pc}^{-2}$ converted from surface brightness in the
unit of $\mathrm{mag}$ $\mathrm{arcsec}^{-2}$) written as Cappellari (2002):

$\mathbf{\Sigma}\left(R^{\prime},\theta^{\prime}\right)=\sum\limits_{\begin{subarray}{c}i=0\end{subarray}}^{N}\frac{L_{j}}{2\pi\sigma^{\prime
2}_{j}q^{\prime}_{j}}\exp\left[-\frac{1}{2\sigma^{\prime
2}_{j}}\left(x_{j}^{\prime 2}+\frac{y_{j}^{\prime 2}}{q^{\prime
2}_{j}}\right)\right],$ (8)

with ($j$ refers to each Gaussian).

${x_{j}^{\prime}=R^{\prime}\sin\left(\theta^{\prime}-\psi_{j}^{\prime}\right)},\hskip
14.22636pt{y_{j}^{\prime}=R^{\prime}\cos\left(\theta^{\prime}-\psi_{j}^{\prime}\right)},$
(9)

where $(R^{\prime},\theta^{\prime})$ are polar coordinates in the sky plane.
$L_{j}$ indicates the observed total luminosity, $q_{j}^{\prime}$ is the
projected flattening and we assume $0\leqslant q_{j}^{\prime}\leqslant 1$,
$\sigma_{j}^{\prime}$ is the scale length along the projected major axis, and
$\psi_{j}^{\prime}$ is the position angle measured counter-clockwise from the
$y^{\prime}$-axis to the major axis of each Gaussian component. We denote:

$\psi_{j}^{\prime}=\psi+\Delta\psi_{j}^{\prime}$ (10)

where $\Delta\psi_{j}^{\prime}$ is the isophotal twist of each Gaussian that
can be measured directly.

For the disk component, we always have
$\psi_{j}^{\prime}=\psi_{\mathrm{disk}}=90^{\circ}$ by aligning its major axis
with $x^{\prime}$, and $\Delta\psi_{j}^{\prime}=0^{\circ}$ for all Gaussian
components by assuming an axisymmetric oblate shape. While for the barred
bulge component which may be triaxial, the isophotal twist of each Gaussian
component is allowed to be different. Here we measure twists with respect to
the disk position angle
($\Delta\psi_{j}^{\prime}=\psi_{j}^{\prime}-\psi_{\mathrm{disk}}$). Details of
fitting are presented in Table 2 and 3 respectively. Results are shown in Fig.
3. There are weak signals of possible spurs which are offset from the major
axis of the inner region of the barred bulge in Fig. 3 (black solid line
contours). The offset-spurs are isophotal signatures corresponding to the
vertically boxy/peanut part of the bar in moderately inclined barred galaxies
Erwin & Debattista (2016).

Figure 3:  (3) The contours of the disk image (black solid line), over-plotted
with contours of the best-fit MGE without twist between different Gaussians
(red dashed line). (3) The MGE fit including twist for the barred bulge.

$j$ |  $L_{j}$ $(\mathrm{L_{\odot}pc^{-2}})$ |  $\sigma^{\prime}_{j}$ $(\mathrm{arcsec})$ |  $q^{\prime}_{j}$ |  $\Delta\psi_{j}^{\prime}$ $(^{\circ})$
---|---|---|---|---
$1$ |  39.380 |  2.400 |  0.547 |  0
$2$ |  42.381 |  5.671 |  0.558 |  0
$3$ |  48.220 |  9.592 |  0.554 |  0
$4$ |  54.451 |  14.300 |  0.555 |  0
$5$ |  74.615 |  21.589 |  0.557 |  0
$6$ |  55.593 |  34.948 |  0.551 |  0
$7$ |  10.840 |  54.956 |  0.568 |  0

Table 2: Details of MGE fit for the disk component in Fig. 3. $j$ is the
number of each individual Gaussian for which, $L_{j}$ is the central flux in
the unit of $(\mathrm{L_{\odot}pc^{-2}})$, $\sigma^{\prime}_{j}$ presents the
size in the unit of $(\mathrm{arcsec})$, $q^{\prime}_{j}$ indicates the
flattening. Position angles of all disk Gaussians components are fixed to be
$90^{\circ}$.

$j$ |  $L_{j}$ $(\mathrm{L_{\odot}pc^{-2}})$ |  $\sigma^{\prime}_{j}$ $(\mathrm{arcsec})$ |  $q^{\prime}_{j}$ |  $\Delta\psi_{j}^{\prime}$ $(^{\circ})$
---|---|---|---|---
$1$ |  1666.431 |  0.864 |  0.544 |  -27.796
$2$ |  7924.190 |  2.136 |  0.593 |  -17.000
$3$ |  1256.558 |  4.848 |  0.480 |  -17.000
$4$ |  1907.776 |  6.499 |  0.480 |  -28.000
$5$ |  117.221 |  10.212 |  0.600 |  -28.000

Table 3: Details of MGE fit for the barred bulge component in Fig. 3.
$\Delta\psi_{j}^{\prime}$ shows the isophotal twist with respect to the disk.

###  3.3 Construction of 3D MGE density

We first deproject each component from 2D MGE surface brightness to 3D
intrinsic density distribution separately, then add the 3D density
distributions of the two components together for a whole galaxy. An oblate
shape and a triaxial ellipsoid are used to describe disk and barred bulge,
respectively. The intrinsic coordinate system denoted by $(x,y,z)$ is aligned
with the galaxy’s principal axes.

The triaxial MGE luminosity density is defined as Cappellari (2002):

$\rho(x,y,z)=\sum\limits_{\begin{subarray}{c}j=0\end{subarray}}^{N}\frac{L_{j}}{(\sigma_{j}\sqrt{2\pi})^{3}q_{j}p_{j}}\exp\left[-\frac{1}{2\sigma^{2}_{j}}\left(x^{2}+\frac{y^{2}}{p^{2}_{j}}+\frac{z^{2}}{q_{j}^{2}}\right)\right],$
(11)

where $N$ is the number of Gaussian functions. $p_{j}=B_{j}/A_{j}$ and
$q_{j}=C_{j}/A_{j}$ are the axial ratios. $A_{j}$, $B_{j}$, and $C_{j}$ are
the intrinsic major, intermediate, and minor axes of the Gaussians.

####  3.3.1 Deprojection of barred bulge: triaxial case

A deprojection could be understood as a transformation between two coordinates
as following:

$\left(\begin{array}[]{l}{x}\\\ {y}\\\
{z}\end{array}\right)=\mathbf{R^{-1}}\cdot\mathbf{P^{-1}}\cdot\left(\begin{array}[]{l}{x^{\prime}}\\\
{y^{\prime}}\\\ {z^{\prime}}\end{array}\right).$ (12)

From parametric Eq. (8), we have the parameters
$(L_{j},q^{\prime}_{j},\sigma^{\prime}_{j},\Delta\psi_{j}^{\prime})$ for each
Gaussian from observations. Given a set of viewing angles
$(\theta,\varphi,\psi)$, we can infer the intrinsic quantities
$(\sigma_{j},p_{j},q_{j})$ through the following equations Van den Bosch et
al. (2008):

$1-q_{j}^{2}=\frac{\delta_{j}^{\prime}\left[2\cos 2\psi_{j}^{\prime}+\sin
2\psi_{j}^{\prime}(\sec\theta\cot\varphi-\cos\theta\tan\varphi)\right]}{2\sin^{2}\theta\left[\delta_{j}^{\prime}\cos\psi_{j}^{\prime}\left(\cos\psi_{j}^{\prime}+\cot\varphi\sec\theta\sin\psi_{j}^{\prime}\right)-1\right]},$
(13) $p_{j}^{2}-q_{j}^{2}=\frac{\delta_{j}^{\prime}\left[2\cos
2\psi_{j}^{\prime}+\sin
2\psi_{j}^{\prime}(\cos\theta\cot\varphi-\sec\theta\tan\varphi)\right]}{2\sin^{2}\theta\left[\delta_{j}^{\prime}\cos\psi_{j}^{\prime}\left(\cos\psi_{j}^{\prime}+\cot\varphi\sec\theta\sin\psi_{j}^{\prime}\right)-1\right]},$
(14)
$u_{j}=\frac{1}{q_{j}^{\prime}}\sqrt{p_{j}^{2}\cos^{2}\theta+q_{j}^{2}\sin^{2}\theta\left(p_{j}^{2}\cos^{2}\varphi+\sin^{2}\varphi\right)}.$
(15)

Where $u_{j}=\sigma_{j}^{\prime}/\sigma_{j}$, $\delta^{\prime
2}_{j}=1-q^{\prime 2}_{j}$, and
$\psi_{j}^{\prime}=\psi+\Delta\psi_{j}^{\prime}$. Note that, in principle, the
absolute value of the position angle $\psi$ of a triaxial structure could be a
free parameter, however, the relative difference of position angle of the
Gaussians $\Delta\psi_{j}^{\prime}$ are measured directly from the MGE
fitting.

All Gaussians are supposed to have the same viewing angles
$(\theta,\varphi,\psi)$, the allowed viewing angles is thus the intersection
of these allowed for each Gaussian component. The allowed range of viewing
angles are restricted by the Gaussian with the minimum of flattening
$q^{\prime}_{j}$ and the maximum difference of twist $\Delta\psi_{j}^{\prime}$
among the Gaussians. Cappellari (2002); Van den Bosch et al. (2008). To avoid
non-physical restriction on the allowed viewing angles caused by a particular
Gaussian, we avoid the Gaussians with too small $q^{\prime}_{j}$ and too large
$\Delta\psi_{j}^{\prime}$ difference, when the error does not change
significantly during the MGE fitting (see Fig.1 in Van den Bosch et al.
(2008)).

We emphasize that if we consider the barred bulge and disk as one rigid body,
deprojection is impossible due to the large difference between
$\Delta\psi_{j}^{\prime}$ of the barred bulge Gaussians and disk Gaussians.
Therefore we need to deproject the barred bulge and disk separately.

####  3.3.2 Deprojection of disk: axisymmetric case

We consider the axisymmetric oblate shape for the disk. We always align the
disk major axis with $x^{\prime}$-axis, thus we have
$\psi_{\mathrm{disk}}=90^{\circ}$ and $\Delta\psi_{j}^{\prime}=0$ for all
Gaussians of the disk. For an axisymmetric oblate system $\varphi$ is
irrelative, therefore Eqs. 13 and 14 are simplified to:

$q^{2}_{j}=\frac{q^{\prime 2}_{j}-\cos^{2}\theta}{\sin^{2}\theta},\hskip
8.5359pt\hskip 8.5359ptp_{j}=1.$ (16)

and

$\sigma_{j}=\sigma^{\prime}_{j},$ (17)

where $(\theta>0^{\circ})$ is the inclination angle of the disk, $q_{j}$ is
the intrinsic flattening and $q^{\prime}_{j}$ is the observed flattening of
each Gaussian component in the disk. Axisymmetric MGE deprojection above is
only valid up until $(\cos^{2}\theta<q^{\prime 2}_{j})$ for all Gaussian
components. It means the minimum inclination is imposed by the flattest
Gaussian in an axisymmetric MGE fit.

In addition, the intrinsic flattening $q$ of disks in late-type galaxies have
a narrow distribution centered at $q\sim 0.26$ Rodríguez & Padilla (2013). We
can roughly derive the inclination of the disk with its observed flattening
$q_{\mathrm{obs}}$, by assuming an intrinsical flattening of $0.26$ from
Eq.16. We denote the disk inclination angle derived in this way as
$\theta_{\mathrm{disk}}^{\mathrm{derive}}$.

###  3.4 Allowed viewing angles

Fig. 4 shows the allowed viewing angles ($\theta$ vs. $\varphi$) of the galaxy
by combining the disk and barred bulge. In Fig. 4a, we show the allowed
regions of viewing angles for the barred bulge, which are intersections of the
allowed viewing angles of each individual Gaussian in the barred bulge.

Then, we consider the major axis of the barred bulge is aligned in the disk
plane. Thus, the inclination angle of the barred bulge should be the same as
the disk. This combination narrows down the allowed inclination angle $\theta$
as shown in Fig. 4b.

The intrinsic position angle $\psi$ of an isolated triaxial system is unknown
as discussed in Van den Bosch et al. (2008). However, in our case, we have a
reference disk which is aligned as $\psi_{\mathrm{disk}}=90^{\circ}$. For the
barred bulge fixed in the disk plane, we restricted it to be:

$|\psi_{\mathrm{bar}}-90^{\circ}|\leq 5^{\circ}.$ (18)

Note that it is
$\psi_{j}^{\prime}=\psi_{\mathrm{bar}}+\Delta\psi_{j}^{\prime}$ that goes into
Eqs. (13-15), the large isophotal twist between the barred bulge and disk is
considered in $\Delta\psi_{j}^{\prime}$ (see Table 3).

We can further constrain the inclination angle $\theta$ with the disk as shown
in Fig. 4d, based on what we roughly derived in Section 3.3.2:

$|\theta-\theta_{\mathrm{disk}}^{\mathrm{derive}}|\leq 10^{\circ},$ (19)

The last constraint we consider refers to the scale height of the disk and
barred bulge. We consider similar thickness of the disk and the barred bulge
at outer regions of the barred bulge:

$|(\sigma q)_{\mathrm{disk}}-(\sigma q)_{\mathrm{bar}}|\leq 10\hskip
0.56917pt\%\hskip 2.84544pt(\sigma q)_{\mathrm{bar}}.$ (20)

Where $(\sigma q)_{\mathrm{bar}}$, and $(\sigma q)_{\mathrm{disk}}$ indicate
size and flattening of the outermost barred bulge Gaussian and a Gaussian of
the disk that has a roughly similar size. We consider $(\sigma
q)_{\mathrm{bar}}$ of outer most Gaussian of the barred bulge (5th Gaussian in
Table. 3) to be close to $(\sigma q)_{\mathrm{disk}}$ of the 4th Gaussian of
the disk (in Table. 2). The allowed regions of viewing angles after imposing
all the constraints are shown in Fig. 4e.

Figure 4:  The parameters space of allowed viewing angles $\theta$ vs.
$\varphi$ for the deprojection of barred bulge + disk. Plots from left- to
right-hand side: $\mathrm{a)}$ all the allowed viewing angles of the barred
bulge, $\mathrm{b)}$ orientations which are also allowed for the disk,
$\mathrm{c)}$ the additional constraint of
$|\psi_{\mathrm{bar}}-90^{\circ}|\leq 5^{\circ}$, $\mathrm{d)}$ the additional
constraint of $|\theta-\theta_{\mathrm{disk}}^{\mathrm{derive}}|\leq
10^{\circ}$, $\mathrm{e)}$ add the last constraint that outer barred bulge and
the inner disk share the similar scale height. The gray lines show the true
value of $\theta^{\mathrm{true}}=60^{\circ}$ and
$\varphi^{\mathrm{true}}=-45^{\circ}$ for the mock galaxy $I_{1}$. Red dot in
last panel indicates a model with ($\theta=58^{\circ}$,
$\varphi_{\mathrm{bar}}=-47^{\circ}$, $\psi_{\mathrm{bar}}=92^{\circ}$) from
the final allowed angles.

###  3.5 3D density distribution

We compare the model inferred density distribution with the true density
distribution of the simulation. In Fig. 5, we show the surface density
distribution of the simulation projected on the $x-y$, $x-z$, and $y-z$
planes. In Fig. 5, we show one model chosen from the allowed regions of
viewing angles with
$(\theta,\varphi_{\mathrm{bar}},\psi_{\mathrm{bar}})=(58^{\circ},-47^{\circ},92^{\circ})$.
We emphasize that $\theta$ is the same for the disk and the barred bulge as we
align the barred bulge in the disk plane. Any $\varphi$ is allowed for disk
because it is assumed to be an oblate axisymmetric structure. We always have
$\psi_{\mathrm{disk}}=90^{\circ}$ because the disk major axis is aligned with
the $x^{\prime}$-axis in the observational plane. The residual between the
true and model inferred density distribution is shown in Fig. 5.

The 3D density distribution from our model generally matches the original
simulation, with a triaxial barred bulge located at the center of a disk. We
obtain the bulge-dominated area along each principal axis
$(x_{\mathrm{bulge}},y_{\mathrm{bulge}},z_{\mathrm{bulge}})$ as $(\sim 4\hskip
2.84544pt\mathrm{kpc},\sim 2\hskip 2.84544pt\mathrm{kpc},\sim 0.78\hskip
2.84544pt\mathrm{kpc})$ for the mock galaxy $I_{1}$, and $(\sim 3.85\hskip
2.84544pt\mathrm{kpc},\sim 2\hskip 2.84544pt\mathrm{kpc},\sim 0.65\hskip
2.84544pt\mathrm{kpc})$ for our model inferred density (See appendix D for
details ). Our model does not have spiral structure. And the disk in our model
is not thin enough in the outer regions, which is a common issue for disk
deprojections if not seen perfectly edge-on. Our model also does not match the
peanut shape of the barred bulge seen edge-on.

The resulting density distributions for other mock galaxies ($I_{2}$ to
$I_{6}$) are shown in appendix A. The deprojection of disks are hard for
galaxies with low inclination angles due to lack of information about the
disk’s intrinsic shapes, thus we do not try the method for galaxies with
$\theta_{\mathrm{disk}}\lesssim 45^{\circ}$, while it is easier for edge-on
galaxies. In contrast, the deprojection of barred bulges prefer galaxies with
lower inclination angles, in which the shapes of the barred bulges are better
revealed.

Overall, our method works reasonably well to match the basic shapes of disk
and barred bulge for galaxies with moderate inclination angles and with
different bar orientations, although the intrinsic shape of the barred bulge
is harder to find when it is projected nearly end-on.

Figure 5:  Comparison of the model inferred density distribution of mock
galaxy $I_{1}$ with the simulation. (a) Surface density distributions of the
simulation projected on $x-y$, $x-z$, and $y-z$ planes. (b) Surface density
distributions of a model with ($\theta=58^{\circ}$,
$\varphi_{\mathrm{bar}}=-47^{\circ}$, $\psi_{\mathrm{bar}}=92^{\circ}$).
Dashed lines in both panels mark the full length of the bar in simulation
($\sim 8\hskip 2.84544pt\mathrm{kpc}$). (c) The residuals of (5) and (5).

##  4 Verification of the deprojected model

Before this model-inferred 3D density distribution could be used to build the
gravitational potential of a dynamical model, we have to figure out how much
uncertainty/bias the model might introduce. We use AGAMA
333https://github.com/GalacticDynamics-Oxford/Agama Vasiliev (2019) to
calculate the potential, forces, and orbits with our model inferred 3D density
and then compare to those calculated with the true density distribution of the
simulation. We illustrate the results with mock galaxy $I_{1}$.

###  4.1 Potential

To obtain the gravitational potential given a density distribution, a
numerical integration is required to solve the Poisson equation.

$\nabla^{2}\Phi(\boldsymbol{x})=4\pi G\rho(\boldsymbol{x}).$ (21)

We freeze the N-body system at the snapshot that is chosen and calculate
potential from the particle distribution using multipole expansion of
spherical harmonics Binney & Tremaine (2008). A similar method is adopted to
compute the potential of model-inferred 3D density distribution given by Fig.
5. The true and model potentials in 2D planes of $x-y$, $x-z$, $y-z$, and
along each principal axis are presented in Fig. 6 and 6 respectively. The
relative difference of the model and true potentials is less than $10\%$ in
all regions.

Figure 6: Comparison of potentials of (6) the simulation in 2D planes and
along the principle axes (solid black lines) and (6) our deprojected 3D model.
(6) The residuals of (6) and (6).

###  4.2 Force

The force is the derivative of potential at each position.

$\vec{F}=-\nabla\Phi,\hskip
14.22636pt|F_{T}|=\sqrt{F_{x}^{2}+F_{y}^{2}+F_{z}^{2}}.$ (22)

Fig. 7 and 7 show the total force in the planes of $x-y$, $x-z$, $y-z$, and
force components along each axis for the true and model derived densities
respectively. The relative difference is less than $15\%$ in all components
and the maximum error occurs around the central point with $20\%$.

Figure 7: Force comparison of (7) the simulation in 2D planes and along
principle axes (solid black lines) and (7) our deprojected 3D model. (7) The
residuals of (7) and (7).

###  4.3 Orbital analysis

We further check if the potential derived from our model can support the key
orbital families of the bar. A rigid pseudo-isothermal DM halo potential (as
described in §2) is added to the true and model potential, respectively. The
scale velocity and scale radius are adopted as $V_{\mathrm{c}}=250\hskip
1.42271pt\mathrm{km}\mathrm{s}^{-1}$ and $R_{c}=15\hskip
1.42271pt\mathrm{kpc}$ respectively. We randomly select $15,000$ initial
conditions corresponding to the positions and velocities of particles from the
snapshot at $\mathrm{t}=2.4\hskip 1.42271pt\mathrm{Gyr}$. We choose particles
inside the corotation radius ($R<4.7$ $\mathrm{kpc}$), then integrate the
orbits in our model inferred potential and in the true potential.

AGAMA uses a modified version of the $8$th order Runge-Kutta integrator
$\mathrm{DOP853}$ Hairer et al. (1993). We calculate the orbits for the chosen
initial conditions in a rotating frame. So the appropriate Coriolis and
centrifugal pseudo-forces are taken into account and determined by the bar
pattern speed $\Omega_{\mathrm{p}}$. We adopt $\Omega_{\mathrm{p}}=39\hskip
2.27626pt\mathrm{km\hskip 1.13791pts^{-1}\hskip 1.13791ptkpc^{-1}}$, which is
the true pattern speed of the simulated bar Shen (2014).

####  4.3.1 Typical orbit families

First, we visually check if orbits with the same initial conditions are
similar in the true and model potentials. Here, orbits are integrated for
$2.5$ $\mathrm{Gyr}$. In the true potential, the time period of a circular
orbit at the end of the bar region is around $0.1$ $\mathrm{Gyr}$. Among the
$15,000$ selected orbits, we randomly plot $200$ orbits with the same initial
conditions in both potentials. Then we check the similarity in appearance of
the true and model orbits in $x-y$, $x-z$, and $y-z$ planes. Orbits with
similar trajectories are considered as the matched ones. We repeat this
process a few more times and average the percentage of matched orbits in each
selection. Generally, we conclude that $85\%$ of all orbits in our model
sample are visually matched with those in the true potential and around $15\%$
are unmatched.

Figure 8: Typical matched orbits in the true potential (black), and in our
model (blue), plotted in the $x-y$, $x-z$, and $y-z$ . Each pair of orbits are
calculated with the same initial conditions.

The most important bar supporting orbit families are prograde $x_{1}$ family,
which are elongated along the major axis of the bar. $x_{1}$ orbits in bars
originate from the same parent as the box orbit family in triaxial potentials
Contopoulos (1980); Schwarzschild (1982); Martinet & de Zeeuw (1988); Valluri
et al. (2016). Other periodic orbit families in the bars are prograde $x_{2}$
and unstable $x_{3}$ orbits, which are elongated perpendicular to the bar and
primarily found at small radii. In addition, retrograde $x_{4}$ orbits are
also perpendicular to the bar and become nearly round at larger radii.

Some typical cases of matched orbits are shown in Fig. 8. The orbits
calculated in the true simulated and our model potential are shown in black
and blue respectively. The orbits in Fig. 8 and 8 are elongated along the bar
major axis. They are $x_{1}$ and box orbit parented by $x_{1}$ orbits,
respectively. While the orbits in Fig. 8 and 8 are retrograde $x_{4}$ orbits
that are elongated perpendicular to the bar. In true and model potentials we
did not find prograde $x_{2}$ and unstable $x_{3}$ families in the inner
region of the bar, similar to Valluri et al. (2016).

There are a few types of resonant orbits usually found in the N-body bar
models: the orbits with $(\Omega_{x},\Omega_{y},\Omega_{z})=(3:-2:0)$ called
’fish/pretzel’ Valluri et al. (2016), the orbits with
$(\Omega_{x},\Omega_{z})=(1:2)$ known as ’banana’ Pfenniger & Friedli (1991)
and the orbits with $(\Omega_{x},\Omega_{y},\Omega_{z})=(3:0:-5)$ called
’brezel’ Portail et al. (2015). The latter two types are proposed as the
backbone of X-shaped structure Patsis et al. (2002); Portail et al. (2015).

We find a few resonant orbits (< 2 %) in the bar regions of our simulation.
Fig. 8 shows a banana orbit in $x-z$ plane that is well-matched in the true
and model potentials. Some other types of resonant orbits can also be found in
our model, but they could be originated from a different starting point than
that similar orbit in the real potential (See appendix B). Some of the
apparently chaotic orbits can sweep similar regions in the true and model
potential, as shown in Fig. 8.

Figure 9: Cartesian frequency map colored by number density, for orbits that
are integrated for $13$ $\mathrm{Gyr}$ in true potential (top row) and our
model (bottom row). Columns from left to right: frequency map for all orbits,
boxes, retrograde $x_{4}$, short-axis tubes (SAT) excluding retrograde
$x_{4}$, long-axis tubes (LAT), and disk orbits. The total number of all
orbits is $15,000$.

####  4.3.2 Frequency analysis and orbit classification

Frequency analysis of orbits was first introduced by Binney & Spergel (1982,
1984) and later developed by Laskar (1990, 1993a). It is a powerful way to
understand the features of orbits in large samples. The fundamental orbital
frequencies are obtained via Fourier transform of their position and velocity
coordinates. A long time integration for orbits is required to get an accurate
frequency map (e.g., $20-50$ orbital periods) Valluri & Merritt (1998);
Valluri et al. (2016). We integrate the orbits in our samples for $13$
$\mathrm{Gyr}$. Then we use the NAFF software
444https://bitbucket.org/cjantonelli/naffrepo/src/master/ Valluri & Merritt
(1998); Valluri et al. (2016) to compute the fundamental frequencies and
perform automated bar orbit classification. Our computation is done in the
Cartesian coordinates, which enables better classification of the bar orbits
than in the cylindrical coordinates (see appendix B in Valluri et al. (2016)).

Figure 10: Projected surface density in $x-y$ and $x-z$ plane extracted from
the 15,000 selected orbits in true potential (top rows) and model potential
(bottom rows). Columns from left to right: surface densities for all orbits,
boxes, retrograde $x_{4}$, short-axis tubes (SAT) excluding retrograde
$x_{4}$, long-axis tubes (LAT), and disk orbits.

Frequency maps of our samples are shown in Fig. 9. The top and bottom rows are
true and model orbits, respectively. The first column from left to right
represent the frequency maps of all the 15,000 orbits in our sample. The
frequency map of orbits in our model has generally the same features as that
of orbits in the true potential, except there is a small frequency offset for
the peak number density. The peak number density in frequency map of orbits in
the true potential is around $\Omega_{x}/\Omega_{z}=0.65$,
$\Omega_{y}/\Omega_{z}=0.87$ and for orbits in our model potential is around
$\Omega_{x}/\Omega_{z}=0.60$, $\Omega_{y}/\Omega_{z}=0.80$.

We classify the orbits in our samples into boxes, retrograde $x_{4}$, short-
axis tubes (SAT), long-axis tubes (LAT), and disk orbits. We emphasize here
SAT refer to all z-tube orbits excluding the retrograde $x_{4}$. The frequency
maps of each class of orbits are shown in Fig. 9 from the second to the fifth
columns (left to right), respectively. Comparison of fractions of different
classes of orbits in our model and in the true potential are shown in Table 4.
Box orbits contribute more than $50$ percent of the orbits, while retrograde
$x_{4}$ orbits contribute $\sim 18$ percent. The SAT and LAT orbits have lower
fractions. The orbital fractions in the true potential and in our model are
generally consistent with each other. We use ($\mathrm{R_{apo}}>4\hskip
2.84544pt\mathrm{kpc}$) to determine the disk orbits, which means the orbits
with apocenter radii larger than half-length of the bar ($4\hskip
2.84544pt\mathrm{kpc}$). The contributions of disk orbits are $18.6\%$ in the
true potential and $21\%$ in our model potential.

We use the frequency drift parameter to determine the chaotic orbits Valluri
et al. (2010). In this method, the orbital time series are divided into two
equal parts and the orbital fundamental frequencies are computed for each time
segment. Since regular orbits have time-independent frequencies, the change in
the frequency measured in the two time segments can be used to determine the
diffusion rate in frequency space Laskar (1993b); Valluri & Merritt (1998).
The frequency drift for each frequency component $\Omega_{i}$ ($i$ = $x$, $y$,
$z$) can be computed as Valluri et al. (2010):

$\log_{10}\left(\Delta
f_{i}\right)=\log_{10}\left|\frac{\Omega_{i}\left(t_{1}\right)-\Omega_{i}\left(t_{2}\right)}{\Omega_{i}\left(t_{1}\right)}\right|.$
(23)

The frequency drift parameter $\log_{10}(\Delta f)$ is defined as the value
associated with the fundamental frequency $\Omega_{i}$ with the largest
amplitude in the Fourier spectrum. A larger frequency drift parameter
indicates the orbit to be more chaotic. To separate chaotic orbits from
regular orbits in frozen N-body potentials, $\log_{10}(\Delta f)>-1.2$ is a
good empirical choice as tested in Valluri et al. (2010). This is also a good
choice for our samples (see appendix C).

Overall $16\%$ of all orbits in the true potential and $17.6\%$ in our model
potential have chaotic features. The contributions of all types of orbits in
our model potential are similar to these in the true potential.

class/type |  true |  model
---|---|---
boxes |  55.2 % |  51.7%
retrograde $x_{4}$ |  17.4 % |  18.3%
SAT |  6.5 % |  7.4%
LAT |  2.3 % |  1.6%
disk orbits |  18.6 % |  21%
sum |  100 % |  100 %
chaotic |  16 % |  17.6 %
regular |  84 % |  82.4 %

Table 4: Classification of $15,000$ selected orbits in the true and our model
potentials. Orbits are classified into boxes, retrograde $x_{4}$, short-axis
tubes (SAT) excluding retrograde $x_{4}$, long-axis tubes (LAT) and disk
orbits. The bottom part of the table presents the fraction of orbits that have
chaotic features in general.

####  4.3.3 Structures constructed with orbits

To check if the structures could be built with the orbits in our model, we
store each orbit in equal-time steps and sum up the density of the 15,000
selected orbits.

The top rows in Fig. 10 show surface densities in $x-z$ and $x-z$ planes
constructed with orbits in the true potential. The bottom rows are with the
orbits in the model potential. Columns from left to right are surface
densities built by all orbit, boxes, retrograde $x_{4}$, SAT, LAT, and disk
orbits, respectively. In general, the structures generated by our model orbits
are similar to the real ones. The retrograde $x_{4}$ orbits support a
structure perpendicular to the bar in face-on view, while the SAT orbits
support a structure elongated along the bar in the face-on map and a
boxy/peanut-shaped structure in the edge-on map. The LAT orbits support a
structure elongated perpendicular to the bar in the face-on map and an
X-shaped structure in the edge-on map. As we discussed in §3, the boxy/peanut-
shaped structure is missing in our model inferred density distribution,
however, as shown in the bottom panel of Fig. 10, the combination of SAT and
LAT orbits in our model potential could still support such a boxy/peanut-
shaped structure in the $x-z$ plane.

In summary, we find that our model potential can support the major orbits
families in a boxy/peanut-shaped bar. Although the fractions of various types
of orbits in our model potential slightly differ from those in the true
potential, we still have all types of orbits that could build the backbone of
the bar.

##  5 Conclusions

We present a method of deprojecting 2D image of a barred galaxy to construct
its 3D intrinsic density distribution. We adopt a two-step process, first, we
decompose the galaxy into an elliptical, Sersic bulge, and a flat, exponential
disk. By subtracting the disk from the original image we get a barred bulge,
then we fit the barred bulge and the disk separately with MGE and deproject
the 2D MGE to get 3D intrinsic density distribution.

We assume the disk is axisymmetric, the barred bulge is triaxial, and the
major axis of the barred bulge is aligned in the disk plane. Additional
constraints are imposed to narrow down the parameter space of allowed viewing
angles: (a) The position angle of the barred bulge is restricted to be close
to that of the reference disk, subtracting the apparent difference
$\Delta\psi^{\prime}$ measured from MGE fitting. (b) the inclination angle of
the disk is restricted to be close to the value obtained by assuming an
intrinsic flattening of $\mathrm{q}_{\mathrm{int}}=0.26$, (c) the intrinsic
scale height of the outer barred bulge is assumed to be similar to the disk in
the same region.

By combining the 3D density distribution of a barred bulge and a disk, we
construct the 3D density distribution of the whole galaxy. We validate the
method by applying it to mock images created from a simulated barred galaxy.
By comparing with the true simulation, we find that:
1) In general, the 3D density distribution from our model matches the true
simulation, with a triaxial barred bulge located at the center of a disk.
However, our model does not match the boxy/peanut shape of the bulge when seen
edge-on. Meanwhile, the disk in our model does not have spiral structure. And
the disk is not very thin at the outer regions in our model, which is a common
issue for disk deprojection if not seen edge-on.
2) We verify this method by comparing the potential and force inferred from
the model constructed 3D density to the true value of the simulation. The
residuals of subtracting the model potential from the true one can be up to
$10\%$ in an extended region within the barred bulge. While the difference in
forces can be as large as $15\%$ in the barred bulge region and $20\%$ around
the very center.
3) We find that 85$\%$ of our sample orbits, including the major families of
bar-supporting orbits, resonant and chaotic orbits, in the model potential
turn out to be very similar to those in the true potential. The unmatched 15%
are mostly chaotic or resonant orbits, their morphologies are easily altered
due to small differences in the gravitational field.
4) The orbits in the model inferred potential distribute similarly in the
frequency maps as those in the true potential. We classify orbits into boxes,
retrograde $x_{4}$, short-axis tubes excluding retrograde $x_{4}$, long-axis
tubes, and disk orbits. The contributions of different classes of orbits in
our model potential are close to these in the true potentials. The short-axis
tubes in our model build an elongated bar with box/peanut-shaped structure,
and the long-axis tubes build a X-shaped structure. These match perfectly the
structures in the true potential.

We have shown that this method can construct 3D density of barred galaxies
from their 2D images, Although our model inferred density does not match the
boxy/peanut shape exactly, the potential still supports the major orbit
families reproducing the boxy/peanut-shaped and X-shaped structures. The
fractions of orbits that building the backbone of bar in our model slightly
differ from those in the true potential. In a dynamical model, we have the
freedom of giving different weights to different orbits by fitting the data.
Thus the fraction of bar-supporting orbits in our model should not in
principle matter once we have sampled all typical orbits. In the future, we
will test to construct Schwarzschild/M2M models of nearby barred galaxies.

## Acknowledgements

We thank Eugene Vasiliev and Monica Valluri for useful discussions. The
research presented here is partially supported by the National Key R&D Program
of China under grant No. 2018YFA0404501; by the National Natural Science
Foundation of China under grant Nos. 12025302, 11773052, 11761131016; by the
“111” Project of the Ministry of Education of China under grant No. B20019;
and by the Chinese Space Station Telescope project, and by the Deutsche
Forschungsgemeinschaft under grant GZ GE 567/5-1(OG). This work made use of
the Gravity Supercomputer at the Department of Astronomy, Shanghai Jiao Tong
University, and the facilities of the Center for High Performance Computing at
Shanghai Astronomical Observatory. LZ acknowledges the support from National
Natural Science Foundation of China under grant No. Y945271001. BT
acknowledges support from CAS-TWAS President’s Fellowship for international
PhD students, awarded jointly by the Chinese Academy of science and The World
Academy of Sciences.

Software: Agama Vasiliev (2019), Astropy Astropy Collaboration et al. (2013,
2018), Galfit Peng et al. (2010), Jupyter Notebook Kluyver et al. (2016),
MgeFit Cappellari (2002), matplotlib Hunter (2007), numpy Harris et al.
(2020), photutils Bradley et al. (2016), NAFF Valluri & Merritt (1998);
Valluri et al. (2016), scipy Virtanen et al. (2020).

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Aguerri et al. (2009) Aguerri J. A. L., Méndez-Abreu J., Corsini E. M., 2009, A&A, 495, 491
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Athanassoula (2003) Athanassoula E., 2003, MNRAS, 341, 1179
  * Barazza et al. (2008) Barazza F. D., Jogee S., Marinova I., 2008, ApJ, 675, 1194
  * Bendinelli (1991) Bendinelli O., 1991, ApJ, 366, 599
  * Binney (1985) Binney J., 1985, MNRAS, 212, 767
  * Binney & Spergel (1982) Binney J., Spergel D., 1982, ApJ, 252, 308
  * Binney & Spergel (1984) Binney J., Spergel D., 1984, MNRAS, 206, 159
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition
  * Bissantz & Gerhard (2002) Bissantz N., Gerhard O., 2002, MNRAS, 330, 591
  * Blaña Díaz et al. (2018) Blaña Díaz M., et al., 2018, MNRAS, 481, 3210
  * Bradley et al. (2016) Bradley L., et al., 2016, Photutils: Photometry tools (ascl:1609.011)
  * Breddels et al. (2013) Breddels M. A., Helmi A., van den Bosch R. C. E., van de Ven G., Battaglia G., 2013, MNRAS, 433, 3173
  * Brown et al. (2013) Brown J. S., Valluri M., Shen J., Debattista V. P., 2013, ApJ, 778, 151
  * Bureau & Athanassoula (2005) Bureau M., Athanassoula E., 2005, ApJ, 626, 159
  * Cappellari (2002) Cappellari M., 2002, MNRAS, 333, 400
  * Cappellari (2008) Cappellari M., 2008, MNRAS, 390, 71
  * Cappellari et al. (2006) Cappellari M., et al., 2006, MNRAS, 366, 1126
  * Contopoulos (1956) Contopoulos G., 1956, Z. Astrophys., 39, 126
  * Contopoulos (1980) Contopoulos G., 1980, A&A, 81, 198
  * Cretton & van den Bosch (1999) Cretton N., van den Bosch F. C., 1999, ApJ, 514, 704
  * Cretton et al. (1999) Cretton N., de Zeeuw P. T., van der Marel R. P., Rix H.-W., 1999, ApJS, 124, 383
  * Debattista & Sellwood (1998) Debattista V. P., Sellwood J. A., 1998, ApJ, 493, L5
  * Emsellem et al. (1994a) Emsellem E., Monnet G., Bacon R., 1994a, A&A, 285, 723
  * Emsellem et al. (1994b) Emsellem E., Monnet G., Bacon R., Nieto J.-L., 1994b, A&A, 285, 739
  * Erwin (2018) Erwin P., 2018, MNRAS, 474, 5372
  * Erwin & Debattista (2016) Erwin P., Debattista V. P., 2016, ApJ, 825, L30
  * Eskridge et al. (2000) Eskridge P. B., et al., 2000, AJ, 119, 536
  * Fragkoudi et al. (2015) Fragkoudi F., Athanassoula E., Bosma A., Iannuzzi F., 2015, MNRAS, 450, 229
  * Friedli & Benz (1993) Friedli D., Benz W., 1993, A&A, 268, 65
  * Gadotti (2009) Gadotti D. A., 2009, MNRAS, 393, 1531
  * Gadotti (2011) Gadotti D. A., 2011, MNRAS, 415, 3308
  * Gadotti et al. (2007) Gadotti D. A., Athanassoula E., Carrasco L., Bosma A., de Souza R. E., Recillas E., 2007, MNRAS, 381, 943
  * Gebhardt et al. (2000) Gebhardt K., et al., 2000, The Astronomical Journal, 119, 1157
  * Hairer et al. (1993) Hairer E., Wanner G., Norsett S., 1993, Solving Ordinary Differential Equations. Springer
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hunt et al. (2013) Hunt J. A. S., Kawata D., Martel H., 2013, MNRAS, 432, 3062
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Kluyver et al. (2016) Kluyver T., et al., 2016, in Loizides F., Schmidt B., eds, Positioning and Power in Academic Publishing: Players, Agents and Agendas. pp 87 – 90
  * Kormendy & Kennicutt (2004) Kormendy J., Kennicutt Robert C. J., 2004, ARA&A, 42, 603
  * Kowalczyk et al. (2017) Kowalczyk K., Łokas E. L., Valluri M., 2017, MNRAS, 470, 3959
  * Laskar (1990) Laskar J., 1990, Icarus, 88, 266
  * Laskar (1993a) Laskar J., 1993a, Celestial Mechanics and Dynamical Astronomy, 56, 191
  * Laskar (1993b) Laskar J., 1993b, Physica D Nonlinear Phenomena, 67, 257
  * Li & Shen (2012) Li Z.-Y., Shen J., 2012, ApJ, 757, L7
  * Li et al. (2011) Li Z.-Y., Ho L. C., Barth A. J., Peng C. Y., 2011, ApJS, 197, 22
  * Li et al. (2018) Li Z.-Y., Shen J., Bureau M., Zhou Y., Du M., Debattista V. P., 2018, ApJ, 854, 65
  * Long et al. (2013) Long R. J., Mao S., Shen J., Wang Y., 2013, MNRAS, 428, 3478
  * Magorrian (1999) Magorrian J., 1999, MNRAS, 302, 530
  * Martinet & de Zeeuw (1988) Martinet L., de Zeeuw T., 1988, A&A, 206, 269
  * Masters et al. (2011) Masters K. L., et al., 2011, MNRAS, 411, 2026
  * Méndez-Abreu et al. (2010) Méndez-Abreu J., Simonneau E., Aguerri J. A. L., Corsini E. M., 2010, A&A, 521, A71
  * Méndez-Abreu et al. (2018) Méndez-Abreu J., Costantin L., Aguerri J. A. L., de Lorenzo-Cáceres A., Corsini E. M., 2018, MNRAS, 479, 4172
  * Menéndez-Delmestre et al. (2007) Menéndez-Delmestre K., Sheth K., Schinnerer E., Jarrett T. H., Scoville N. Z., 2007, ApJ, 657, 790
  * Monnet et al. (1992) Monnet G., Bacon R., Emsellem E., 1992, A&A, 253, 366
  * Noordermeer & van der Hulst (2007) Noordermeer E., van der Hulst J. M., 2007, MNRAS, 376, 1480
  * Patsis et al. (2002) Patsis P. A., Skokos C., Athanassoula E., 2002, MNRAS, 337, 578
  * Peng et al. (2010) Peng C. Y., Ho L. C., Impey C. D., Rix H.-W., 2010, AJ, 139, 2097
  * Pfenniger & Friedli (1991) Pfenniger D., Friedli D., 1991, A&A, 252, 75
  * Portail et al. (2015) Portail M., Wegg C., Gerhard O., Martinez-Valpuesta I., 2015, MNRAS, 448, 713
  * Portail et al. (2017) Portail M., Gerhard O., Wegg C., Ness M., 2017, MNRAS, 465, 1621
  * Qin et al. (2015) Qin Y., Shen J., Li Z.-Y., Mao S., Smith M. C., Rich R. M., Kunder A., Liu C., 2015, ApJ, 808, 75
  * Richstone & Tremaine (1985) Richstone D. O., Tremaine S., 1985, ApJ, 296, 370
  * Rodríguez & Padilla (2013) Rodríguez S., Padilla N. D., 2013, MNRAS, 434, 2153
  * Rybicki (1987) Rybicki G. B., 1987, in de Zeeuw P. T., ed., IAU Symposium Vol. 127, Structure and Dynamics of Elliptical Galaxies. p. 397, doi:10.1007/978-94-009-3971-4_41
  * Schwarzschild (1982) Schwarzschild M., 1982, ApJ, 263, 599
  * Shen (2014) Shen J., 2014, in Feltzing S., Zhao G., Walton N. A., Whitelock P., eds, Vol. 298, Setting the scene for Gaia and LAMOST. pp 201–206, doi:10.1017/S1743921313006376
  * Shen et al. (2010) Shen J., Rich R. M., Kormendy J., Howard C. D., De Propris R., Kunder A., 2010, ApJ, 720, L72
  * Sheth et al. (2005) Sheth K., Vogel S. N., Regan M. W., Thornley M. D., Teuben P. J., 2005, ApJ, 632, 217
  * Stark (1977) Stark A. A., 1977, ApJ, 213, 368
  * Valluri & Merritt (1998) Valluri M., Merritt D., 1998, ApJ, 506, 686
  * Valluri et al. (2004) Valluri M., Merritt D., Emsellem E., 2004, ApJ, 602, 66
  * Valluri et al. (2010) Valluri M., Debattista V. P., Quinn T., Moore B., 2010, MNRAS, 403, 525
  * Valluri et al. (2016) Valluri M., Shen J., Abbott C., Debattista V. P., 2016, ApJ, 818, 141
  * Van den Bosch et al. (2008) Van den Bosch R. C. E., van de Ven G., Verolme E. K., Cappellari M., de Zeeuw P. T., 2008, MNRAS, 385, 647
  * Vasiliev (2013) Vasiliev E., 2013, MNRAS, 434, 3174
  * Vasiliev (2019) Vasiliev E., 2019, MNRAS, 482, 1525
  * Vasiliev & Valluri (2020) Vasiliev E., Valluri M., 2020, ApJ, 889, 39
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Wang et al. (2013) Wang Y., Mao S., Long R. J., Shen J., 2013, Monthly Notices of the Royal Astronomical Society, 435, 3437
  * Weinberg (1985) Weinberg M. D., 1985, MNRAS, 213, 451
  * Weiner et al. (2001) Weiner B. J., Sellwood J. A., Williams T. B., 2001, ApJ, 546, 931
  * Zhao (1996) Zhao H., 1996, MNRAS, 283, 149
  * Zhu et al. (2018a) Zhu L., et al., 2018a, Nature Astronomy, 2, 233
  * Zhu et al. (2018b) Zhu L., et al., 2018b, MNRAS, 473, 3000
  * Zou et al. (2014) Zou Y., Shen J., Li Z.-Y., 2014, ApJ, 791, 11
  * de Zeeuw & Franx (1989) de Zeeuw T., Franx M., 1989, ApJ, 343, 617

##  Appendix A Deprojection of more cases with different orientations

To show the effectiveness of our method for galaxies with various
observational orientations, we obtain the 3D density maps for mock galaxies
with different viewing angles as listed in Table 1. The results are shown in
Figs. 11 and 12. The first and second columns show the mock images and
deprojected density maps respectively. The 3D density inferred for these
galaxies generally match the true 3D density distribution of the simulation
well, in a similar manner to $I_{1}$ which we described in detail.

Figure 11:  Deprojection of more mock galaxies, with different angles of
$\varphi_{\mathrm{bar}}$: 11) Mock data $I_{2}$; 11) mock data $I_{3}$. Dashed
lines mark the full length of the bar in simulation ($\sim 8\hskip
2.84544pt\mathrm{kpc}$).

Figure 12:  Deprojection of edge-on mock galaxies with different angles of
$\varphi_{\mathrm{bar}}$: 12) Mock data $I_{4}$ ; 12) mock data $I_{5}$ ; 12)
mock data $I_{6}$. Dashed lines mark the full length of the bar in simulation
($\sim 8\hskip 2.84544pt\mathrm{kpc}$).

##  Appendix B Unmatched orbits

Fig. 13 shows some typical unmatched orbits. The first row shows the orbits
which are boxes in $x-y$ plane under the true and model potentials, while in
$x-z$ plane they have different shapes. In Fig. 13 our model potential
generates a brezel orbit while it happens inversely in Fig. 13. Also in Fig.
13 and 13 the fish orbit is generated once in our model and once in the true
potential for different starting points. The third row shows the orbits which
have different trajectories in all planes. Around $1/3$ of visually chaotic
orbits follow totally different trajectories in the true and model potentials,
Fig. 13 shows such an example. Overall most of the unmatched orbits are from
chaotic or resonant orbits.

Figure 13: Typical unmatched orbits in the true potential (black), and in our
model (blue), plotted in the $x-y$, $x-z$, and $y-z$. Each pair of orbits is
calculated with the same initial conditions. Resonant orbits and chaotic
orbits are easy to be unmatched, they are sensitive to the minor changes of
the potential.

##  Appendix C Chaotic orbits criterion

We examined if $\log_{10}(\Delta f)>-1.2$ used in Valluri et al. (2010, 2016)
is a fair criterion to define chaotic orbits in our potential. In Fig. 14, we
show the histogram distribution of $\log_{10}(\Delta f)$ for the orbits in our
model and in the true potential, which are similar to each other, especially
at the high end with $\log_{10}(\Delta f)>-1.2$. We thus conclude it is a fair
criterion.

Figure 14: Distributions of frequency drift parameter $\log_{10}(\Delta f)$
for $15,000$ orbits in the true (red) and model (blue) potentials. The dashed
line indicates the threshold value of $\log_{10}(\Delta f)=-1.2$. Orbits with
$\log_{10}(\Delta f)>-1.2$ are classified as chaotic orbits.

##  Appendix D measurement of the bulge-dominated region

Fig. 15 shows the surface density distribution along the three principal axes,
blue is the simulation, and orange is our model inferred density from mock
galaxy $I_{1}$.

It is not straightforward to determine the bar size in three directions for a
given density. To quantitatively compare the bulge-dominated region in
simulation and in our model inferred density, we use the surface density of
the disk in our model as a reference. We define the bulge-dominated region
along each direction as the positions where
$2\times\Sigma_{\mathrm{disk}}=\Sigma_{\mathrm{total}}$. We obtain
$(x_{\mathrm{bulge}},y_{\mathrm{bulge}},z_{\mathrm{bulge}})$ as $(\sim 4\hskip
2.84544pt\mathrm{kpc},\sim 2\hskip 2.84544pt\mathrm{kpc},\sim 0.78\hskip
2.84544pt\mathrm{kpc})$ for the simulation and $(\sim 3.85\hskip
2.84544pt\mathrm{kpc},\sim 2\hskip 2.84544pt\mathrm{kpc},\sim 0.65\hskip
2.84544pt\mathrm{kpc})$ for our model inferred density.

Figure 15:  Surface density distribution of the simulation (blue) and our
model for mock galaxy $I_{1}$ (orange); along the major axis (top),
intermediate axis (middle), and minor axis (bottom). The green solid line and
dashed line indicate $\Sigma_{\mathrm{disk}}$ and
$2\times\Sigma_{\mathrm{disk}}$ of our model. We define the bulge-dominated
region at the radius where $\Sigma_{\mathrm{total}}$ =
$2\times\Sigma_{\mathrm{disk}}$. The blue and orange vertical lines mark the
size of the bulge-dominated region in the simulation and in our model,
respectively.