{'fraction_non_alphanumeric': 0.03678269989512061, 'fraction_numerical': 0.022124556759726313, 'mean_word_length': 4.304238410596026, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 5, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Deprojection of external barred galaxies from photometry

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.


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 ℎ 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 Hunt et al. 2013;Portail et al. 2015Portail et al. , 2017. There are several commonly used implementations of the Schwarzschild's orbitsuperposition 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 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 = √︂ 2 + 2 2 + 2 2 which requires the axial ratios ≤ ≤ 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 ( < 60 • ) 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 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.

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 × 10 10 with 10 6 equal-mass particles. A rigid pseudo-isothermal DM halo potential is adopted Φ = 1 2 2 0 ln 1 + 2 2 c , in which the scale velocity and scale radius are 0 = 250 kms −1 and = 15 kpc, respectively. The bar has a half-length of ∼ 4 kpc and rotates with a pattern speed of Ω p ∼ 39 km s −1 kpc −1 (corotation radius ∼ 4.7 kpc). The end-to-end separation between the outer two edges of the X-shaped structures is ∼ 4 kpc along the major axis and ∼ 2.4 kpc along the vertical minor axis (Li & Shen 2012).
Throughout the paper, we use the coordinate ( , , ) to describe the intrinsic 3D structure, where , , are aligned with the long, intermediate and short axes of the galaxy. While we use the coordinate ( , ) to describe the projected structure to the 2D observational plane, and is along the line-of-sight. The orientation of a projection is specified by the viewing angles ( , , ). and 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 ( = 90 • , = 0 • ), ( = 90 • , = 90 • ) and ( = 0 • , = 0 • , ..., 90 • ). The angle 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): where matrix P is responsible for the projection onto the sky plane defined as: and matrix R expresses the rotation around the line-of-sight by angle : To simplify the description of multiple components of a galaxy, we always align the major axis of the disk with axis and minor with . Note that this could be different from the natural coordinate, i.e. aligning with north, when dealing with real observational data. Here is defined from the axis, thus we always have disk = 90 • 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 Mpc and then create a surface brightness map with a spatial resolution of 1 arcsec pixel −1 . At this distance, 1 kpc equals 5 arcsec. We produce several mock images from the simulation adopting different viewing angles of and as listed in Table 1, which are labeled as true values. To keep the disk major axis aligned with -axis, all of our mock galaxies have true = 90 • . 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 1 .

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 Inclination angle true is the orientation of the disk, true = 90 • means edgeon and true = 0 • means face-on. true describes the orientation of the bar, true = 90 • means side-on and true = 0 • means end-on. To keep the disk major axis aligned with -axis, we have true = 90 • for all cases.
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.

2D bulge-to-disk decomposition
We use GALFIT 3.0.5 1 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: where is the scale length of disk and = 2 + 2 2 1 2 , ( , ) are aligned with the major and minor axis of the disk, and is the axial ratio.
We use a Sérsic function to describe the surface brightness of an elliptical bulge: Where is the effective radius (or half-light radius), Σ is the surface brightness at the effective radius , and is called the "Sérsic index", which controls the concentration of the profile. When is small, it has a shallow inner profile and a steep truncation at a large radius. Inversely, When is large, it has a steep inner profile and a highly extended outer wing. is not a free parameter, it coupled to by (2 ; ) = 1 2 Γ(2 ), Γ and are the Gamma function and lower incomplete Gamma function, respectively. Here is defined as: = / , where and are the minor and major axis of the elliptical bulge in the observational plane. and are aligned with the major 1 https://users.obs.carnegiescience.edu/peng/work/galfit/ galfit.html where proj bar is the position angle of the elliptical bulge measured counter-clockwise from the -axis. Fig. 1 presents the 2D surface brightness profile of (a) the mock galaxy 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 70 arcsec, and brighter at 70 arcsec. Our best-fit parameters are = 20.90 arcsec, Σ 0 = 11.63 mag arcsec −2 , = 0.51 for the disk component and = 1.093, = 5.88 arcsec, Σ = 11.69 mag arcsec −2 , = 0.46, and proj bar = 71 • for the elliptical bulge. Note that we have disk = 90 • , thus the position angle of the bar is −19 • different from the disk.

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 MgeFit software from Cappellari (2002) 2 . The MGE describes the surface brightness (in the unit of L sun pc −2 converted from surface brightness in the unit of mag arcsec −2 ) written as (Cappellari 2002): with ( refers to each Gaussian).
where ( , ) are polar coordinates in the sky plane. indicates the observed total luminosity, is the projected flattening and we assume 0 1, is the scale length along the projected major axis, and is the position angle measured counter-clockwise from the -axis to the major axis of each Gaussian component. We denote: where Δ is the isophotal twist of each Gaussian that can be measured directly. For the disk component, we always have = disk = 90 • by aligning its major axis with , and Δ = 0 • 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 (Δ = − 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. 3b (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).

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 ( , , ) is aligned with the galaxy's principal axes.   Fig. 3a. is the number of each individual Gaussian for which, is the central flux in the unit of (L pc −2 ), presents the size in the unit of (arcsec), indicates the flattening. Position angles of all disk Gaussians components are fixed to be 90 • . The triaxial MGE luminosity density is defined as (Cappellari 2002): where is the number of Gaussian functions. = / and = / are the axial ratios. , , and are the intrinsic major, intermediate, and minor axes of the Gaussians.
Where = / , 2 = 1 − 2 , and = + Δ . Note that, in principle, the absolute value of the position angle of a triaxial structure could be a free parameter, however, the relative difference of position angle of the Gaussians Δ are measured directly from the MGE fitting.
All Gaussians are supposed to have the same viewing angles ( , , ), 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 and the maximum difference of twist Δ 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 and too large Δ 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 Δ of the barred bulge Gaussians and disk Gaussians. Therefore we need to deproject the barred bulge and disk separately.

Deprojection of disk: axisymmetric case
We consider the axisymmetric oblate shape for the disk. We always align the disk major axis with -axis, thus we have disk = 90 • and Δ = 0 for all Gaussians of the disk. For an axisymmetric oblate system is irrelative, therefore Eqs. 13 and 14 are simplified to: and where ( > 0 • ) is the inclination angle of the disk, is the intrinsic flattening and is the observed flattening of each Gaussian component in the disk. Axisymmetric MGE deprojection above is only valid up until (cos 2 < 2 ) 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 of disks in late-type galaxies have a narrow distribution centered at ∼ 0.26 (Rodríguez & Padilla 2013). We can roughly derive the inclination of the disk with its observed flattening obs , by assuming an intrinsical flattening of 0.26 from Eq.16. We denote the disk inclination angle derived in this way as derive disk . Fig. 4 shows the allowed viewing angles ( vs. ) 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 as shown in Fig. 4b.

Allowed viewing angles
The intrinsic position angle 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 disk = 90 • . For the barred bulge fixed in the disk plane, we restricted it to be: Note that it is = bar + Δ that goes into Eqs. (13-15), the large isophotal twist between the barred bulge and disk is considered in Δ (see Table 3).
We can further constrain the inclination angle with the disk as shown in Fig. 4d, based on what we roughly derived in Section 3.3.2: 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: Where ( ) bar , and ( ) 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 ( ) bar of outer most Gaussian of the barred bulge (5th Gaussian in Table. 3) to be close to ( ) 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.

3D density distribution
We compare the model inferred density distribution with the true density distribution of the simulation. In Fig. 5a, we show the surface density distribution of the simulation projected on the − , − , and − planes. In Fig. 5b, we show one model chosen from the allowed regions of viewing angles with ( , bar , bar ) = (58 • , −47 • , 92 • ). We emphasize that is the same for the disk and the barred bulge as we align the barred bulge in the disk plane. Any is allowed for disk because it is assumed to be an oblate axisymmetric structure. We always have disk = 90 • because the disk major axis is aligned with the -axis in the observational plane. The residual between the true and model inferred density distribution is shown in Fig. 5c. 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 ( bulge , bulge , bulge ) as (∼ 4 kpc, ∼ 2 kpc, ∼ 0.78 kpc) for the mock galaxy 1 , and (∼ 3.85 kpc, ∼ 2 kpc, ∼ 0.65 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 ( 2 to 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 disk 45 • , 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.

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 3 (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 1 .

Potential
To obtain the gravitational potential given a density distribution, a numerical integration is required to solve the Poisson equation.
3 https://github.com/GalacticDynamics-Oxford/Agama 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. 5b. The true and model potentials in 2D planes of − , − , − , and along each principal axis are presented in Fig. 6a and 6b respectively. The relative difference of the model and true potentials is less than 10% in all regions.

Force
The force is the derivative of potential at each position. Fig. 7a and 7b show the total force in the planes of − , − , − , 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%.

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 c = 250 kms −1 and = 15 kpc respectively. We randomly select 15, 000 initial conditions corresponding to the positions and velocities of particles from the snapshot at t = 2.4 Gyr. We choose particles inside the corotation radius ( < 4.7 kpc), then integrate the orbits in our model inferred potential and in the true potential.
AGAMA uses a modified version of the 8th order Runge-Kutta integrator 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 Ω p . We adopt Ω p = 39 km s −1 kpc −1 , which is the true pattern speed of the simulated bar (Shen 2014). 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 − , − , and − 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. The most important bar supporting orbit families are prograde 1 family, which are elongated along the major axis of the bar. 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 2 and unstable 3 orbits, which are elongated perpendicular to the bar and primarily found at small radii. In addition, retrograde 4 orbits are also perpendicular to the bar and become nearly round at larger radii.

Typical orbit families
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. 8a and 8b are elongated along the bar major axis. They are 1 and box orbit parented by 1 orbits, respectively. While the orbits in Fig. 8c and 8d are retrograde 4 orbits that are elongated perpendicular to the bar. In true and model potentials we did not find prograde 2 and unstable 3 families in the inner region of the bar, similar to Valluri et al. (2016).
We find a few resonant orbits (< 2 %) in the bar regions of our simulation. Fig. 8e shows a banana orbit in − 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. 8f.

Frequency analysis and orbit classification
Frequency analysis of orbits was first introduced by Binney & Spergel (1982, 1984 and later developed by Laskar (1990Laskar ( , 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 Gyr. Then we use the NAFF software 4 (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)).
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 4 https://bitbucket.org/cjantonelli/naffrepo/src/master/ 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 Ω /Ω = 0.65, Ω /Ω = 0.87 and for orbits in our model potential is around Ω /Ω = 0.60, Ω /Ω = 0.80. We classify the orbits in our samples into boxes, retrograde 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 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 4 orbits contribute ∼ 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 (R apo > 4 kpc) to determine the disk orbits, which means the orbits with apocenter radii larger than half-length of the bar (4 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 timeindependent 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 Ω ( = , , ) can be computed as (Valluri et al. 2010): The frequency drift parameter log 10 (Δ ) is defined as the value associated with the fundamental frequency Ω 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 (Δ ) > −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.

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 − and − 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 4 , SAT, LAT, and disk orbits, respectively. In general, the structures generated by our model orbits are similar to the real ones. The retrograde 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 − 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.

CONCLUSIONS
We present a method of deprojecting 2D image of a barred galaxy to construct its 3D intrinsic density distribution. We adopt a twostep 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, Figure 9. Cartesian frequency map colored by number density, for orbits that are integrated for 13 Gyr in true potential (top row) and our model (bottom row).
Columns from left to right: frequency map for all orbits, boxes, retrograde 4 , short-axis tubes (SAT) excluding retrograde 4 , long-axis tubes (LAT), and disk orbits. The total number of all orbits is 15, 000. Figure 10. Projected surface density in − and − 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 4 , short-axis tubes (SAT) excluding retrograde 4 , long-axis tubes (LAT), and disk orbits. 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 Δ 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 q 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 4 , short-axis tubes excluding retrograde 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/peanutshaped 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. a b Figure A1. Deprojection of more mock galaxies, with different angles of bar : a) Mock data 2 ; b) mock data 3 . Dashed lines mark the full length of the bar in simulation (∼ 8 kpc).

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. A1 and A2. 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 1 which we described in detail. Fig. B1 shows some typical unmatched orbits. The first row shows the orbits which are boxes in − plane under the true and model potentials, while in − plane they have different shapes. In Fig. 1a our model potential generates a brezel orbit while it happens inversely in Fig. 1b. Also in Fig. 1c and 1d 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. 1e shows such an example. Overall most of the unmatched orbits are from chaotic or resonant orbits.

APPENDIX C: CHAOTIC ORBITS CRITERION
We examined if log 10 (Δ ) > −1.2 used in Valluri et al. (2010Valluri et al. ( , 2016 is a fair criterion to define chaotic orbits in our potential. In Fig. C1, we show the histogram distribution of log 10 (Δ ) for the orbits in our model and in the true potential, which are similar to each other, a b c Figure A2. Deprojection of edge-on mock galaxies with different angles of bar : a) Mock data 4 ; b) mock data 5 ; c) mock data 6 . Dashed lines mark the full length of the bar in simulation (∼ 8 kpc). especially at the high end with log 10 (Δ ) > −1.2. We thus conclude it is a fair criterion. Fig. D1 shows the surface density distribution along the three principal axes, blue is the simulation, and orange is our model inferred density from mock galaxy 1 .

APPENDIX D: MEASUREMENT OF THE BULGE-DOMINATED REGION
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 × Σ disk = Σ total . We obtain ( bulge , bulge , bulge ) as (∼ 4 kpc, ∼ 2 kpc, ∼ 0.78 kpc) for the simulation and (∼ 3.85 kpc, ∼ 2 kpc, ∼ 0.65 kpc) for our model inferred density. a b c d e f Figure B1. Typical unmatched orbits in the true potential (black), and in our model (blue), plotted in the − , − , and − . 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.