The SAMI Galaxy Survey: The Internal Orbital Structure and Mass Distribution of Passive Galaxies from Triaxial Orbit-superposition Schwarzschild Models

Dynamical models are crucial for uncovering the internal dynamics of galaxies; however, most of the results to date assume axisymmetry, which is not representative of a significant fraction of massive galaxies. Here, we build triaxial Schwarzschild orbit-superposition models of galaxies taken from the SAMI Galaxy Survey, in order to reconstruct their inner orbital structure and mass distribution. The sample consists of 161 passive galaxies with total stellar masses in the range 109.5–1012 M ⊙. We find that the changes in internal structures within 1R e are correlated with the total stellar mass of the individual galaxies. The majority of the galaxies in the sample (73% ± 3%) are oblate, while 19% ± 3% are mildly triaxial and 8% ± 2% have triaxial/prolate shape. Galaxies with logM⋆/M⊙>10.50 are more likely to be non-oblate. We find a mean dark matter fraction of f DM = 0.28 ± 0.20, within 1R e. Galaxies with higher intrinsic ellipticity (flatter) are found to have more negative velocity anisotropy β r (tangential anisotropy). β r also shows an anticorrelation with the edge-on spin parameter λRe,EO , so that β r decreases with increasing λRe,EO , reflecting the contribution from disk-like orbits in flat, fast-rotating galaxies. We see evidence of an increasing fraction of hot orbits with increasing stellar mass, while warm and cold orbits show a decreasing trend. We also find that galaxies with different (V/σ – h 3) kinematic signatures have distinct combinations of orbits. These results are in agreement with a formation scenario in which slow- and fast-rotating galaxies form through two main channels.

The assembly history of a galaxy is thought to be one of the major factors that determines its internal kinematic structure (e.g., White 1979;Fall & Efstathiou 1980;Park et al. 2019) and so observations of the in-ternal kinematic structure should give an indication of a galaxy's past.
Our current understanding of galaxy formation suggests that massive galaxies form in a two-phase process (e.g., Naab et al. 2009;Oser et al. 2010).During the first phase, at high redshift, they grow by a rapid episode of in-situ star formation, resulting in compact massive systems.After z ≈ 2, these massive, log 10 (M /M ) > 10.5, compact galaxies are predicted to be quiescent and grow mostly by accreting mass through gas-poor galaxy mergers that add stars mainly to their outskirts.
Early-type galaxies (ETGs) have been separated into two classes, based on their stellar kinematics: fast rotators and slow rotators (e.g.Emsellem et al. 2004Emsellem et al. , 2007;;Cappellari et al. 2007;Emsellem et al. 2011).Cappellari (2016) suggested that these two classes also indicate two major channels of galaxy formation where fast-rotating ETGs start their life as star-forming disks and evolve through a set of processes dominated by gas accretion, bulge growth and quenching.In contrast, slow-rotating ETGs assemble near the centers of massive halos, via intense star formation at high redshift, and evolve from a set of processes dominated by gas-poor mergers.However, Naab et al. (2014) showed that the detailed formation history of a galaxy cannot be constrained from the slow-fast rotator classification alone, but when combined with the higher-order kinematic signatures, different merger scenarios can be distinguished.
In order to understand the evolutionary history of galaxies, we need a detailed analysis of its intrinsic structure.The Schwarzschild orbit-superposition method (Schwarzschild 1979) is a powerful dynamical modelling technique that allows dynamical substructures in galaxies to be revealed.Several different implementations of the Schwarzschild method, with varying degrees of symmetry, have been described (e.g.Cretton et al. 1999;Gebhardt et al. 2003;Valluri et al. 2004;van den Bosch et al. 2008;Vasiliev & Athanassoula 2015;Vasiliev & Valluri 2020;Neureiter et al. 2021).The Schwarschild method has been used to model supermassive black holes (van der Marel et al. 1998;Verolme et al. 2002;Gebhardt et al. 2003;Valluri et al. 2004;Krajnović et al. 2009;Rusli et al. 2013;Seth et al. 2014;Thater et al. 2017Thater et al. , 2019;;Liepold et al. 2020;Quenneville et al. 2021), the internal orbital structures of globular clusters (van de Ven et al. 2006;Feldmeier-Krause et al. 2017), earlytype galaxies (Cappellari et al. 2006;Thomas et al. 2007;van de Ven et al. 2008;Thomas et al. 2014;Fahrion et al. 2019;Poci et al. 2019;Jin et al. 2020;den Brok et al. 2021;Thater et al. 2022) and recently expanded to galaxies of all morphologies (Vasiliev & Athanassoula 2015;Zhu et al. 2018b,c;Vasiliev & Valluri 2020;Lipka & Thomas 2021).The orbit distributions obtained by these models have also been used to identify different dynamical components in these stellar systems (e.g.van de Ven et al. 2006;Cappellari et al. 2007;van den Bosch et al. 2008;Lyubenova et al. 2013;Breddels & Helmi 2014;Krajnović et al. 2015).Zhu et al. (2018b) separated orbits into four different components: a cold component with near circular orbits (with strong rotation), a hot component with near radial orbits (characterized by random motions), a warm component in-between (characterized by weak rotation) and a counter-rotating component (similar to the warm and cold components, but with reversed angular momentum).The inferred internal orbital distributions were then used to reconstruct the observed photometry and stellar kinematics of each component.However, the majority of these studies only had a few objects available (less than 30 galaxies).A large sample of galaxies, observed with good radial coverage and spatial resolution, is required in order to understand the average evolution history of the general galaxy population.
In the last two decades, Integral Field Spectroscopy (IFS) surveys such as SAURON (Spectroscopic Areal Unit for Research on Optical Nebulae; de Zeeuw et al. 2002), ATLAS 3D (Cappellari et al. 2011), CALIFA (Calar Alto Legacy Integral Field Array survey; Sánchez et al. 2012), SAMI (Sydney-Australian-Astronomical-Observatory Multi-object Integral-Field Spectrograph) Galaxy Survey (Croom et al. 2012;Bryant et al. 2015;Croom et al. 2021), MASSIVE (Ma et al. 2014), MaNGA (Mapping Nearby Galaxies at Apache Point Observatory; Bundy et al. 2015) and the Fornax 3D survey (Sarzi et al. 2018) have provided us with rich observational datasets of galaxies, allowing their structure and evolution to be investigated in detail through the mapping of stellar kinematics across individual galaxies.These IFS surveys have made possible the use of techniques such as Schwarschild orbit-superposition method to dynamically decompose IFS observations to estimate the internal mass distribution, intrinsic stellar shapes and orbit distributions of galaxies across the Hubble sequence (e.g., Zhu et al. 2018a,b,c;Zhuang et al. 2019;Jin et al. 2020;Aquino-Ortíz et al. 2020).Zhu et al. (2018c) studied a sample of 250 galaxies in the CALIFA survey, with total stellar masses between 10 8.5 and 10 12 M , spanning all morphological types.About 95% of the galaxies in their sample had stellar kinematic maps with R max > 1R e , and ∼ 8% with R max > 3R e .They found that, within 1 R e , galaxies have more stars in warm orbits than in either cold or hot orbits.Similar results were also found in a sample of 149 early-type galaxies in the MaNGA survey (Jin et al. 2020), with stellar masses ranging between 10 9.9 and 10 11.8 M and observations up to 1.5 -2.5 R e per galaxy.These studies also found that the changes of internal structures within 1R e are correlated with the stellar mass of the galaxies.
The number of galaxies considered for Schwarzschild model studies to date has been limited and they have often not incorporated higher-order kinematic moments to further constrain the orbital models.Higher-order kinematic signatures are defined as the deviations from a Gaussian line-of-sight velocity distribution (LOSVD).When the LOSVD is parametrized as a Gauss-Hermite series (van der Marel & Franx 1993;Gerhard 1993), its skewness and excess kurtosis are parametrized by the coefficients of the 3 rd -and 4 th -order Hermite polynomials (h 3 and h 4 , respectively).Given the connection between the higher-order stellar kinematic moments and a galaxy's assembly history (Naab et al. 2014), their inclusion in dynamical modelling can help distinguish between different formation scenarios.
In this paper we will apply Schwarzschild modelling to the SAMI Galaxy Survey (Croom et al. 2012;Bryant et al. 2015;Owers et al. 2017) to investigate the evolutionary histories of passive galaxies by studying their internal structures.The SAMI Galaxy Survey data allows us to study the internal orbits of a significant number of galaxies for the first time and allows us to further constrain the Schwarzschild models by adding information on the higher-order kinematic moments.

OBSERVATIONS
The Sydney-AAO Multi-object Integral field spectrograph (SAMI) Galaxy Survey is a large, optical Integral Field Spectroscopic (Croom et al. 2012;Bryant et al. 2015;Owers et al. 2017) survey of low-redshift (0.04 < z < 0.095) galaxies covering a broad range in stellar mass, 7 < log 10 (M /M ) < 12, morphology and environment.The sample, with ≈ 3000 galaxies, is selected from the Galaxy and Mass Assembly survey (GAMA; Driver et al. 2011) regions (field and group galaxies), as well as eight additional clusters to probe higher-density environments (Owers et al. 2017).
The SAMI instrument (Croom et al. 2012), on the 3.9m Anglo-Australian telescope, consists of 13 "hexabundles" (Bland-Hawthorn et al. 2011;Bryant et al. 2014), across a 1-degree field of view.Each hexabundle consists of 61 individual 1. 6 fibres, and covers a ∼ 15 diameter region on the sky.In the typical configuration, 12 hexabundles are used to observe 12 science targets, with the 13 th one allocated to a secondary standard star used for calibration.Moreover, SAMI also has 26 individual sky fibres, to enable accurate sky subtraction for all observations without the need to observe separate blank sky frames.The SAMI fibres are fed to the dualbeam AAOmega spectrograph (Sharp et al. 2006).
2.1.IFS Spectra and kinematic maps SAMI data consist of 3D data cubes: two spatial dimensions and a third spectral dimension.
The wavelength coverage is from 3750 to 5750 Å in the blue arm, and from 6300 to 7400 Å in the red arm, with a spectral resolution of R = 1812 (2.65 Å full-width half maximum; FWHM) and R = 4263 (1.61 Å FWHM), respectively (van de Sande et al. 2017a), so that two data cubes are produced for each galaxy target.
Each galaxy field was observed in a set of approximately seven 30 minute exposures, that are aligned together by fitting the galaxy position within each hexabundle with a two-dimensional Gaussian and by fitting a simple empirical model describing the telescope offset and atmospheric refraction to the centroids.The exposures are then combined to produce a spectral cube with regular 0.5 spaxels, with a median seeing of 2.1 .More details of the Data Release 3 reduction can be found in Croom et al. (2021) 1 .
Stellar kinematic measurements were derived using the penalized pixel fitting code (pPXF; Cappellari & Emsellem 2004;Cappellari 2017), after combining the blue-and red-arm spectra by matching their spectral resolution.A detailed description of the method used to derive the stellar kinematic measurement can be found in van de Sande et al. (2017a,b).In particular, for our analysis, we use the Voronoi-binned kinematic measurements.Bins are adaptively generated to contain a target S/N of 10 Å−1 , using the Voronoi binning code of Cappellari & Copin (2003).
The available stellar kinematic measurements consist of 2D maps of stellar rotational velocity V , velocity dispersion σ, and the high kinematic orders (h 3 and h 4 ).In addition, each kinematic map has kinematic position angle and FWHM of the Point Spread Function (PSF -taken from a star observed at the same time as the galaxies) provided.

Multi Gaussian Expansion profiles and effective radius
Multi-Gaussian Expansion (MGE; Emsellem et al. 1994;Cappellari 2002) profile fits for the SAMI Galaxy 1 Reduced data-cubes and stellar kinematic data products for all galaxies are available on: https://datacentral.org.au.
Survey are produced from the r−band photometry by D' Eugenio et al. (2021).The MGE method consists of a series expansion of galaxy images using 2D Gaussian functions.This method enables us to take the PSF into account; given a value of the inclination and assuming an intrinsic shape, the MGE model can be deprojected analytically, which is orders of magnitude faster than the general, integral-based method.
The fits are applied to re-analysed Sloan Digital Sky Survey (SDSS; York et al. 2000) images for GAMA galaxies, reprocessed as described in Hill et al. (2011), as well as VST/ATLAS (VLT Survey Telescope -ATLAS; Shanks et al. 2015) and SDSS DR9 (Ahn et al. 2012) observations for cluster galaxies, with VST/ATLAS data reprocessed as described in Owers et al. (2017).The images are square cutouts with 400 side, centerd on the center of the galaxy, and the MGE fits are calculated using MgeFit (Cappellari 2002) and the regularisation feature described in Scott et al. (2009).For a given galaxy, each Gaussian component has its PA fixed to that of the host's major axis.As such, the stellar mass distribution is assumed to be axisymmetric in projection, but can be intrinsically triaxial.A more detailed description of the fitting process can be found in D 'Eugenio et al. (2021).From the MGE best fit, we use the projected luminosity, size, and flattening of each Gaussian component to model the surface density of each galaxy and to deproject the stellar component into a 3D density.The effective radius, R e , used here is that of the major axis in the r-band.The semi-major axis values were taken from MGE fits.

Stellar Mass
Stellar masses are estimated assuming a Chabrier (2003) initial mass function (IMF), from the K-corrected g− and the i− magnitudes using an empirical proxy developed from GAMA photometry (Taylor et al. 2011;Bryant et al. 2015).For cluster galaxies, stellar masses are derived using the same approach (Owers et al. 2017).We use the photometric stellar masses for our analysis in order to be consistent with previous SAMI studies and to have consistent comparisons with previous results in the literature (e.g. from CALIFA and MaNGA).

Sample Selection
We use data from the final SAMI data release (described in the Data Release 3 publication Croom et al. 2021).This data release consists of 3068 unique galaxies.Of these, we have MGE profiles from D 'Eugenio et al. (2021) for 2957 galaxies (r−band images are not available for some galaxies or they have been affected by a bright star in the field of view).Following van de Sande et al. (2017a), we exclude all galaxies whose kinematics are influenced by mergers, that have strong bars or that have a bright secondary object within one effective radius in their stellar velocity field.This leaves us with 2834 galaxies with stellar kinematic and MGE measurements.
We exclude all galaxies with masses below log 10 (M /M ) = 9.5, because the incompleteness of the stellar kinematic sample is larger than 50% of the SAMI galaxy survey sample observed in this mass range.We further exclude 433 galaxies where R e < 2 (due to their spatial size being smaller than the instrumental spatial resolution).This leaves us with 1649 galaxies.
Following the recommendations of van de Sande et al. (2017a), for each galaxy we select spaxels that meet the following quality criteria: 2017a) is for measurements with S/N < 20 Å−1 and σ obs < 70 km/s.We cautiously include these in this analysis and increase the errors on the measurements that do not meet this criterion to down-weight their contributions.The 1589 galaxies that meet these criteria are shown in Fig. 1.
In this paper we focus on passive galaxies, because the long-term goal of this project is to study the effects of galaxy environment on passive galaxies (Santucci et al. in prep).We use the SAMI spectroscopic classification presented in Owers et al. (2019) to select a homogeneous sample.The SAMI spectroscopic classification labelled galaxies as star-forming, passive, or Hδ-strong, using the absorption-and emission-line properties of each SAMI spectrum.We select 738 passive galaxies.

Radial coverage and spatial sampling selection
We compare the spatial resolution and radial extent of our sample to the sample from (Zhu et al. 2018b) who used CALIFA data to derive orbital parameters using the Schwarzschild method.SAMI Voronoi bins are generated to contain a target S/N of 10 Å−1 .Since the target S/N is the only requirement for the bins, individual spaxels of 0.5 are left unbinned when they meet this requirement.For these single-spaxel bins, the covariance is larger (since they are smaller than the SAMI spatial resolution).In Fig. 1 we show the number of Voronoi bins within 1R e versus the radial coverage avalable (in units of R e ) for the 1589 SAMI galaxies (in grey) that meet our quality criteria.CALIFA galaxies (in purple; from Zhu et al. 2018b) have a similar distribution in number of bins to SAMI, however their bins were generated with different criteria (their minimum S/N = 20 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 We calculate the marginalized fractions of galaxies to the total number in each sample, by mass and size, and show them in the top and left panels of the figure.Grey lines are for SAMI galaxies, while the violet lines are for CALIFA galaxies.The CALIFA and the SAMI samples have similar distributions in Voronoi bins and radial coverage, although there are more CALIFA galaxies with measurements up to 2Re.For this analysis we select galaxies in the top right corner (Rmax > 1Re and Voronoi bins > 85).
and their spaxel size is consistent with their spatial resolution), therefore a direct comparison is not possible.CALIFA and SAMI also show a similar distribution in radial coverage, although there are more CALIFA galaxies with measurements up to 2R e .In this analysis, the first in a series, we select a highquality subsample of SAMI galaxies, identified by good spatial resolution and good radial coverage (top right corner of Fig. 1).This region is selected as the optimal compromise between best quality data and reasonable sample size, and corresponds to galaxies with 85 Voronoi bins within 1R e and R max > R e .More details about the radial coverage tests we performed can be found in Appendix A.
This quality cut gives us a sample of 179 passive galaxies.We visually inspect the galaxies in this sample using HSC images and exclude the face-on strongly barred galaxies that were not identified as barred from the square cutouts used for the MGE modelling.This cut gives us a final sample of 161 galaxies.These are shown in Fig. 2 and used hereafter in this analysis.The majority of the galaxies in our sample are earlytype galaxies (∼ 85%), ∼ 11% are S0/Early-spirals and ∼ 4% are late-type galaxies (visual morphological classification from Cortese et al. 2016).We note that our final sample is biased toward galaxies that are more massive and larger than the general SAMI passive population.This bias is caused by selecting galaxies with at least 85 Voronoi bins within 1 R e .Effective radius, Re, versus stellar mass.Blue circles are the passive galaxies in the SAMI sample with log 10 (M /M ) > 9.5 and Re > 2 (738), orange squares are the galaxies included in the final sample (161).We calculate the marginalized fractions of galaxies with respect to the total number in each sample, by mass and size, and show them in the top and left panels of the figure.Blue lines are for the passive galaxies in the SAMI Galaxy Survey, while the orange lines are for our final sample.The two samples are slightly different in the marginalized mass and size distributions, so that we have higher fractions of massive and large galaxies in the final sample compared to the initial sample.This is due to selecting galaxies with more than 85 Voronoi bins.

SCHWARZSCHILD ORBIT-SUPERPOSITION TECHNIQUE
3.1.Schwarzschild's models and free parameters We use the Schwarzschild orbit-superposition technique (Schwarzschild 1979) to model our individual galaxies, using the implementation from van den Bosch et al. (2008) with correct orbital mirroring from Quenneville et al. (2022).This code allows us to model triax-ial stellar systems2 , while many previous applications of this technique assume axisymmetry.There are three main steps required to create a Schwarzschild model: 1. Construct a model for the underlying gravitational potential; 2. Calculate a representative library of orbits using the gravitational potential previously modelled; 3. Find a combination of orbits that can reproduce the observed kinematic maps and luminosity distribution.
These steps are fully described in van den Bosch et al. (2008) and Zhu et al. (2018a) and are summarized in the following subsections.

Gravitational Potential
The model gravitational potential of each galaxy is generated using the combination of three components: a stellar and a dark matter distribution and a central super-massive black hole.The triaxial stellar component mass is calculated from the best-fit two-dimensional MGE luminosity density (from D'Eugenio et al. 2021) which is de-projected assuming the orientation in space of the galaxy, described by three viewing angles (θ, φ, ψ), to obtain a three-dimensional luminosity density.The space orientation (θ, φ, ψ) can be converted directly to the intrinsic shape (p i , q i , u i ), where p i = B i /A i , q i = C i /A i and u i = σ obs Gauss,i /σ Gauss,i .A i , B i , C i represent the major, medium and minor axes of the 3D triaxial Gaussian component and σ Gauss,i represents the size of each Gaussian component.Moreover, the flattest Gaussian component, having the minimum observed flattening q min , dictates the allowed space orientation for the de-projection, so that we can take (p min , q min , u min ) as our free parameters.The 3D density defined by these intrinsic shapes is then converted into a stellar mass distribution using a radially constant stellar mass-to-light ratio M /L (note that M /L is a free parameter in our modelling).The corresponding stellar gravitational potential Φ is calculated using the classical formula from Chandrasekhar (1969).
The dark matter halo distribution is assumed to follow a spherical Navarro-Frenk-White profile (NFW; Navarro et al. 1996).The mass, M 200 (mass enclosed within a radius, R 200 , where the average density is 200 times the critical density), in a NFW dark matter halo is determined by two parameters.These are the concentration parameter, c, and the fraction of dark matter within R 200 , f = M 200 /M (where M 200 is as defined above and M is the total stellar mass).
The spatial resolution of SAMI data is poorer than the influence radius of the black hole, so its mass leaves no imprint on the stellar kinematic maps and therefore does not affect our results.We therefore fix the black hole mass to the value derived from the stellar velocity dispersion, measured within an aperture of 1R e , assuming the relation between black hole mass and the stellar velocity dispersion of a galaxy from McConnell et al. (2011).
Combining the components used to describe the gravitational potential, we have six free parameters (stellar mass-to-light ratio, M /L, the intrinsic shape of the flattest Gaussian component (p min , q min , u min ), the dark matter halo concentration, c, and dark matter fraction, f ) that must be determined.To determine these bestfit parameters for each galaxy, we run an optimized grid-based parameter search as described in Zhu et al. (2018a) and summarized in Sec.3.4.

Orbit library
To fit a model to our observed data we need an orbit library.To create the orbit library we use a separable triaxial potential, where all orbits are regular and conserve three integrals of motion (energy E, second integral I 2 and third integral I 3 ) which can be calculated analytically.Our Schwarzschild implementation considers four different families of orbits: three types of tube orbits (short axis tubes, outer and inner long axis tubes) and box orbits.We create initial conditions for our orbits by sampling from the three integrals of motion.We refer to van den Bosch et al. (2008) for the details of the orbit sampling.
The number of points we sample across the three integrals is n E × n θ × n R = 21 × 10 × 7, where n E , n θ , n R are the number of intervals taken across the energy E, the azimuthal angle θ and radius R on the (x, z) plane.However, this orbit library includes mostly short axis tubes, long axis tubes and a relatively low fraction of box orbits in the inner region.Since box orbits are essential for creating triaxial shapes, we construct an additional set of box orbits.Box orbits always touch equipotentials (Schwarzschild 1979), so they can be described by combining the energy E with two spherical angles (θ and φ).The number of points included in the box orbit set are n E × n θ × n φ = 21 × 10 × 7.
We add an additional set of orbits to account for retrograde stars commonly found in early-type galaxies (Bender 1988;Kuijken et al. 1996).This set contains 21 × 10 × 7 orbits to describe the initial conditions for counter-rotating orbits.To summarize, we use three sets of 21 × 10 × 7 orbits: a typical set of (E, I 2 , I 3 ), a box orbits set of (E, θ, φ) and a counter-rotating set of also (E, −I 2 , I 3 ).
As in van den Bosch et al. (2008) and Zhu et al. (2018b), we dither every orbit to give 5 3 orbits by perturbing the initial conditions slightly, in order to smooth the model.The orbit trajectories created by the dithering will be co-added to form a single orbit bundle in our orbit library.
We then use Schwarzschild's method to weight the various orbit contributions to the LOSVD in each bin to construct a model with observational parameters that can be fit to the data (the description of how kinematic maps are fitted can be found in Zhu et al. 2018a).The quantities that will be compared to observations are spatially convolved with the same PSF as the observations.The model and the observed values are then divided by the observational error so that a χ 2 comparison is achieved.The weights are determined by the van den Bosch et al. (2008) implementation, using the Lawson & Hanson (1974) non-negative least squares (NNLS) implementation.

Best-fit model
In order to find the best-fit model, which contains six free parameters, we run a grid based parameter search.We use a parameter grid with intervals of 0.5, 0.1, 0.2, 0.05, 0.05 and 0.01 in M /L, log(c), log(f ), q min , p min and u min , respectively, and perform an iterative search for the best-fitting models.After each iteration, the best-fit model is selected by using a χ 2 comparison.The best-fit model is defined as the model with minimum kinematic χ 2 : where V n mod , σ n mod , h n 3, mod and h n 4, mod are the model values for each bin n, V n obs , σ n obs , h n 3, obs and h n 4, obs are the observed values in each bin and V n obserr , σ n obserr , h n 3, obserr and h n 4, obserr represent the observational errors.N kin is the number of bins in the kinematic maps.We define a confidence level around that minimum value and select all the models whose χ 2 is within that confidence level: kin , as we use V, σ, h 3 and h 4 as model constraints, and N par is the number of free parameters (6 here).We then create new models around the existing models with lower kinematic χ 2 values by walking two steps in every direction of the parameter grid from each of the selected models.In this way, the searching process goes in the direction of smaller χ 2 on the parameter grid, and it stops when the minimum χ 2 model is found.Next, we continue the iteration by using a larger value of χ 2 s , to ensure all the models within 1σ confidence are calculated before the iteration finishes.The values of χ 2 s are chosen empirically so that it is neither too small (finding only local minimums) nor too large.For the final step, we reduce the parameter intervals by half to get a better estimate of the best-fit parameters.The models whose χ 2 are within the confidence level are included for calculating the statistical uncertainties of the model parameters for single data analysis.The maximum and minimum values of the parameters or properties in these models are treated as upper and lower limits in 1σ error regions.
The kinematic maps for the best-fit models of example galaxies 9403800123, 9011900793, 220465 and 9008500323 are presented in Fig. 3, Fig. 4, Fig. 5 and Fig. 6.We selected these four galaxies as representative of the sample, with 9403800123 being a non edgeon oblate galaxy (with 255 Voronoi bins within 1R e ), 9011900793 an edge-on oblate galaxy (with 87 Voronoi bins within 1R e ), 220465 a triaxial galaxy (142 Voronoi bins within 1R e ) and 9008500323 a prolate galaxy (with 104 Voronoi bins within 1R e ).Even when the spatial sampling is low, as in the case of 9011900793, the model is able to reproduce the observed kinematic maps well (χ 2 red = 2.22 for galaxy 9403800123, χ 2 red = 1.72 for galaxy 9011900793, χ 2 red = 1.79 for galaxy 220465 and χ 2 red = 1.99 for galaxy 9008500323 3 ).We also show the explored parameter grids and the obtained internal mass distribution, orbit circularity, triaxiality and tangential anisotropy for the best fits of these four galaxies in Appendix B. These parameters are fully described in the following section.

RESULTS
3 The reduced χ 2 is defined as 4N kin −Npar , with χ 2 calculated following Eq. 1.The values of χ 2 red are not always equal to 1 for the best-fit models of the galaxies in our sample.This is because the input kinematic maps of the galaxies in our sample were not symmetrized.Therefore, comparing the observed maps to the model maps, which are symmetric, can result in values of χ 2 red higher than 1.In this section and the next we present the results we obtain modelling a sample of 161 passive galaxies in the SAMI Galaxy Survey with the Schwarzschild orbitsuperposition technique.For each galaxy we explore a range in parameter space by building on average 1250 different models.This is consistent with previous analyses that used an iterative grid search in ∼6 dimensions.For example Jin et al. (2019Jin et al. ( , 2020) ) required 1000 to 2000 separate Schwarzschild models per galaxy to be run.By comparing the 2D maps of the flux and kinematic parameters derived from each model and observations we determine the best-fit parameters.From the best-fit model we derive the intrinsic properties of the inner mass distribution (for both stellar and dark matter components), intrinsic stellar shape (axis ratios and ellipticity), velocity anisotropy and the orbit circularity distribution.We take as our best-fit values the parameters calculated at or averaged within an aperture of 1R e , depending on the parameter.Uncertainties on the measured values are calculated using Monte Carlo realizations, as described in Appendix C, combined with the 1σ confidence levels for the parameters fluctuations from the best-fit model that we describe in Sec.3.4.

Inner mass distribution
The total mass (M tot ) radial distribution is one of the fundamental parameters of the Schwarzschild model, which includes a stellar component and a dark matter component (M dark ).A black hole mass component is included as well, but not discussed here as its contribution to the total mass distribution is negligible.The distribution of the fraction of dark matter (f DM = M dark /M tot ) within 1R e for the galaxies in our sample is shown in Fig. 7.The average value of the dark matter fraction is 0.28, with a standard deviation of 0.20.Similar to Cappellari et al. (2013), we fit a quadratic function to the f DM versus stellar mass distribution.The best-fit rela- tion follows f DM ∼ 0.10 + 0.17 × (log M /M − 10.59) 2 , although the 1-σ scatter along this relation is as high as δf DM = 0.24.
Above a stellar mass of log(M /M ) ∼ 10.75 we see a hint of an increasing f DM as a function of stellar mass.To test whether this trend is statistically significant, we use the Kendall's correlation coefficient τ , using the Python package scipy.stats.kendalltau(Virtanen et al. 2019).This correlation coefficient is robust to small sample sizes.A τ value close to 1 indicates strong correlation, whereas a value close to −1 indicates strong anti-correlation.For galaxies with log(M /M ) 10.75 we find a value of τ = 0.17, with a probability of correlation of 99.73%.While the trend of increasing fraction of dark matter with increasing stellar mass is mild, it is significant at the 3-σ level.

Intrinsic stellar shape
Next, we investigate the intrinsic shapes of the galaxies in our sample.As shown in Sec.3.2, three parameters are used to model the dynamically-based intrinsic stellar shape of each galaxy: p, q and u.The intrin-sic shape has been shown to be connected to various other galaxy properties such as: stellar mass (Sánchez-Janssen et al. 2010), luminosity (Sánchez-Janssen et al. 2016), spin parameter (e.g.Foster et al. 2017), mean stellar population age (van de Sande et al. 2018) and its environment (Fasano et al. 2010;Rodríguez et al. 2016).Furthermore, theoretical simulations suggest that intrinsic shape depends on a galaxy's formation history (Jesseit et al. 2009;Li et al. 2018b,a).
Here, in particular, we analyse the triaxial parameter T Re , calculated at 1R e and defined as: (2) We show an example of the best-fit intrinsic shape parameters p, q and T as a function of radius in Appendix B, Fig. 30.Based on the triaxiallity parameter T Re , we separate galaxies into three groups according to their dynamically-based intrisic shape: oblate (T Re = 0), prolate (T Re = 1) and triaxial (T Re = 0, 1).In Fig. 8 we show the triaxial parameter T Re as a function of stellar mass log(M /M ).The majority of the galax- The best-fit model maps (χ 2 red = 1.79) accurately reconstruct the structures seen in the observations, not only for the velocity and velocity dispersion maps, but also for h3 and h4.
ies in our sample are close to oblate (118 out of 161 galaxies; 73% ± 3%), 30 galaxies (19% ± 3%) show evidence of being mildly triaxial (0.1 < T Re ≤ 0.3) and 13 galaxies (8% ± 2%) have triaxial/prolate shapes (with T Re > 0.3).There is evidence of a slight increase of triaxiality with increasing stellar mass (τ = 0.1), however, this trend is only significant at a 1-σ level (with a probability of 82.96%).However, if we consider the fraction of galaxies that have T Re > 0.1 (non-oblate galaxies), we find a clear increase of the fraction with stellar mass, with a sharp change at ∼ 10 10.50 M /M , with the fraction of non-oblate galaxies increasing from 12% ± 4% to 29% ± 2% at this mass.
Non-oblate galaxies are often dispersion-dominated, with their shape reflecting the anisotropic velocity dispersion.In contrast, oblate galaxies may have varying degrees of rotation support and anisotropy (e.g.Kireeva & Kondratyev 2019).We analyse the distribution of the velocity dispersion anisotropy in the next section.

Velocity anisotropy
Velocity dispersion anisotropy parameters (e.g.β r , β z ) are widely used as indicators of the underlying orbit distribution of a galaxy.However, various definitions and approaches exist in the literature.The velocity dispersion anisotropy parameter used in more recent literature, β z , is in cylindrical coordinates and has been used in particular to describe the global anisotropy in fastrotating axisymmetric galaxies (Cappellari et al. 2007).This parameter measures the velocity anisotropy along the radius on the disk plane, in cylindrical coordinates, following the idea of cylindrically aligned stellar velocity ellipsoids ellipsoids in oblate galaxies.However, for triaxial galaxies β z (< R e ) will have a contribution from both circular orbits (which have cylindrically-aligned velocity dispersion ellipsoids) as well as radial and box orbits (which have spherically-aligned velocity dispersion ellipsoids).Recent results (Thater et al. 2022) show that the velocity dispersion ellipsoids for the elliptical galaxy NGC 6958 are more closely aligned with spherical coordinates.The misalignment between the measured ellipsoids and the cylindrical coordinates can reach angles as high as 80 • .This misalignment can even occur in disk galaxies, most notably, our own Milky Way (Büdenbender et al. 2015;Hagen et al. 2019).Following Thater et al. (2022), we measure the misalignment of the velocity ellipsoids for the galaxies in our sample and find that they are more closely aligned with spherical coordinates.For this reason, we focus on the radial velocity anisotropy parameter, β r , in the results presented here.For completeness, we also include the results for β z in Appendix D.
We define the velocity anisotropy parameter β r , in spherical coordinates, following Binney & Tremaine (2008): with (r, θ, φ) the standard spherical coordinates, and with σ k the velocity dispersion along the direction k at a given location inside the galaxy.The summation defines how we computed this quantity from our Schwarzschild models.M n is the mass contained in each of the N polar grid cells in the meridional plane of the model, and σ k,n is the corresponding mean velocity dispersion along the direction k.
We calculate the value of β r within 1R e , excluding the inner regions (r < 2 ) since this region is affected by atmospheric seeing.β r > 0 indicates radial anisotropy, β r < 0 indicates tangential anisotropy and β r = 0 indicates isotropy.Figure 9 shows the derived values of β r , for each galaxy, as a function of intrinsic ellipticity.Here, we derive ε using the intrinsic flattening, q Re , from the best-fit model of the galaxy, measured at 1R e ; ε intr = 1 − q Re .In general, galaxies with high ellipticity (flat galaxies, ε intr > 0.7) are close to isotropic or tangentially anisotropic (supported by rotation).We also find that radially anisotropic galaxies are typically more massive than tangentially anisotropic galaxies.The blue solid line is the parabolic best-fit to the data -fDM ∼ 0.10 + 0.17 × (log M M − 10.59) 2 .The fit is calculated by taking into account the errors on the data points, so that it is critically determined by the data points with small errors.The shaded region represents the error on the best-fit.

Spin Parameter
The proxy for the spin parameter, λ r , has previously been used to separate slow-rotating galaxies from fastrotating galaxies (Emsellem et al. 2007(Emsellem et al. , 2011;;Cappellari 2016).We use the Cortese et al. ( 2016) definition of the spin parameter to calculate λ r for each galaxy: where i refers to each spaxel within the ellipse with semimajor axis R e and ellipticity ε, F i is the corresponding flux of the i th spaxel, V i is its stellar velocity, σ i is the velocity dispersion and R i is the semi-major axis of the ellipse in which the spaxel lies.Since λ r is calculated within 1R e , it will be referred to as λ Re hereafter.
For completeness, we also measure the ratio of ordered to random motion V /σ, also measured within 1R e , using the definition from Cappellari et al. (2007): Re ) as a function of stellar mass.Galaxies with TRe = 0 are classified as oblate, galaxies with TRe = 1 as prolate and those in-between as triaxial.Grey dashed lines represent TRe = 0.1, TRe = 0.3 and TRe = 0.8.The majority of the galaxies in our sample are oblate, with a few galaxies with non-oblate shape.The average values of the triaxiality parameter for each of the 4 mass bins are shown as dark blue squares, with the error bars marking the 1-σ scatter.There is a weak increase in the triaxiality parameter with increasing stellar mass.The percentage of galaxies that are non-oblate (TRe > 0.1) increases with increasing stellar mass, going from 12%±4% below 10 10.50 M /M to 29%±2% above this mass.
Results obtained using V /σ are similar to those obtained for λ Re and are shown in Appendix E.
Inclination has a strong impact on the observed λ Re and V /σ quantities, in particular when the viewing angle is close to face-on (e.g.Binney et al. 1990).While inclination corrections are now commonly applied to λ Re measurements (e.g.Querejeta et al. 2015;van de Sande et al. 2018;Falcón-Barroso et al. 2019;del Moral-Castro et al. 2020;Fraser-McKelvie et al. 2021) these methods cannot be applied for slow rotating or triaxial galaxies (for a detailed discussion see van de Sande et al. 2021a).Our triaxial Schwarzschild models now allow us, irrespective of galaxy type, to deproject each galaxy to a consistent edge-on view and reconstruct a best-fit internal orbital distribution for that viewing angle.
In order to reconstruct the edge-on maps, we recalculate and store the orbit library for each galaxy, with a specific projection.Schwarzschild models take into account the PSF of the observations when reproducing the kinematics.To construct 2D maps without the impact of seeing within the Schwarzschild routine, we set the PSF FWHM to 0.01 for the model to use when projecting the galaxy.
Once we have constructed the edge-on projected maps, we measure the spin parameter within 1R e by applying Equation 6.In order to produce results comparable to observations, we remeasure the MGE model on our edge-on projected maps to derive R e , using the MgeFit python package (Cappellari 2002).We then derive the ellipticity by finding the model isophote with area A = πR 2 e , and use its ellipticity as the galaxy ellipticity (D'Eugenio et al. 2021).We show the derived edge-on λ Re,EO values as a function of the edge-on intrinsic ellipticity from our MGE fit, ε intr,EO , in Fig. 10, color-coded by their velocity anisotropy β r .The magenta line corresponds to the relation β z = 0.65 × ε for edge-on galaxies as in Cappellari et al. (2007).
We find that λ Re,EO increases with increasing intrinsic ellipticity.In particular, galaxies that have low values of λ Re,EO are rounder than galaxies with higher values of λ Re,EO .Moreover, we find that galaxies that are radially anisotropic (positive values of β r ) show lowto mid-values of ellipticity and λ Re,EO , while galaxies with high ellipticity and λ Re,EO are more isotropic or tangentially anisotropic.This is seen more clearly when a locally weighted regression algorithm (LOESS -Cappellari et al. 2013) is applied to the data to recover any mean underlying trend in β r (Fig. 10, panel b).In general, the variation in β r seems to mostly be driven by the spin parameter, λ Re,EO .
The anti-correlation between λ Re,EO and β r can be seen in Fig. 11.Testing the correlation using the Kendall's correlation coefficient τ , we find a value of τ = −0.27,with a probability of correlation of 99.99% that β r decreases with increasing λ Re,EO .This means that fast-rotating galaxies are, as expected, more tangentially anisotropic than slow-rotating systems, which are more radially anisotropic.

Orbital structure
Stellar orbits can be characterized by two main properties: the time-averaged radius r, representing the size of each orbit, and the circularity λ z = L z /(r×V c ), where L z is the time averaged z-component of the orbit's angular momentum (xv y − yv x ), r = x 2 + y 2 + z 2 , and The denominator represents the angular momentum of a typical circular orbit associated with the original orbit.Using the ratio of these two angular momentum terms, we can quantify the orbit circularity.|λ z | = 1 represents highly-rotating short-axis tube orbits (circular orbits), while λ z = 0 represents mostly box or radial orbits.Taking the radius, r, and the circularity, λ z , of each orbit, and considering their weights given by the solution from the best-fit model, we can use the orbit circularity distribution in the phase space to obtain the probability density of orbits within 1R e , for each galaxy.
Figure 12 shows the overall orbit circularity distribution for all the galaxies in our sample, sorted by increasing stellar mass (shown in the top x-axis).The orbit circularity distribution is calculated by integrating the probability distribution of λ z over all radii within 1R e and normalizing it to unity.The color of each square represents the normalized density, ω, of the orbits on the phase space.We divide the orbits into four broad categories (similar to Zhu et al. 2018a,c): cold orbits, λ z ≥ 0.80 (close to circular orbits); warm orbits, 0.25 < λ z < 0.80 (short-axis tube orbits with a component of rotation but also contribution of random motions); hot orbits, −0.25 ≥ λ z ≤ 0.25 (mostly box orbits and long-axis tube orbits); counter-rotating orbits, λ z < −0.25, (similar to the warm and cold components, but with opposite rotation).Overall, the amount of hot The two parameters are anti-correlated, so that βr decreases with increasing λRe,EO.This means that fast-rotating galaxies are more likely to be tangentially anisotropic.
orbits increases with increasing stellar mass, while the number of warm and cold orbits becomes smaller with increasing mass.
To better visualize these trends with stellar mass, we calculate the luminosity-weighted fractions of each com-ponent within 1R e as a function of stellar mass in Fig. 13, panel a.We also divide the sample into 4 mass bins with 29 galaxies each and we show the median values for each mass bin as bold points.We find a clear increase in the fraction of hot orbits with increasing stellar mass (τ = 0.16, with a probability of correlation of 99.71%), while the fraction of warm orbits decreases with increasing stellar mass (τ = −0.19,with a probability of correlation of 99.95%), both of them showing a large scatter.In particular, the fraction of hot and warm orbits seem to have a sharp change above log M /M = 11.The fraction of cold orbits only have a weak correlation with mass (τ = −0.10,with a probability of correlation of 94.21%), declining towards more massive galaxies.The fraction of counter-rotating orbits does not seem to depend on stellar mass (τ = 0.05, with a probability of correlation of 61.44%).
We also explore the correlation of the fractions of the orbital components with the bulge to total flux ratio, B/T in panel b, with the intrinsic spin parameter λ Re,EO in panel c and with the intrinsic ellipticity ε intr in panel d.B/T ratios are calculated from the r-band photometry, performing a 2D photometric bulge-disk decomposition (Barsanti et al. 2021 for the decomposition of cluster galaxies and Casura et al., in prep, for the galaxies in the GAMA region).Only 97 galaxies in our sample have reliable B/T values for the 2 component decomposition.The orbital fractions show a correlation with the B/T ratios similar to that with stellar mass.Figure 12.Overall orbit circularity distribution (calculated by integrating the probability distribution of λz, over all radii within 1Re and normalizing it to unity), for all the galaxies in our sample, sorted by increasing stellar mass (shown in the top x-axis).The color indicates the normalized density, ω, of the orbits on the phase space.The orbits are divided into four categories: cold orbits (λz ≥ 0.80), warm orbits (0.25 < λz < 0.80), hot orbits (−0.25 ≥ λz ≤ 0.25) and counter-rotating orbits (λz < −0.25).Darker colors indicate higher probabilities as illustrated by the color bar.The right-hand panel shows the average orbit-circularity distribution within the mass range.Overall, the fraction of hot orbits seems to increase with increasing stellar mass, while the fraction of warm and cold orbits becomes smaller with increasing mass.
Looking at λ Re,EO and ε intr the orbital fractions have similar trends: hot orbits decrease with increasing λ Re,EO and ε intr , warm orbits increase with increasing λ Re,EO and ε intr and cold orbits show an increase in the fractions, while there is a significant change (τ = −0.21,with a probability of correlation of 99.99%) in the fraction of counter-rotating orbits only with λ Re,EO , so that the fraction decreases with increasing λ Re,EO .In particular, we note that the trends with λ Re,EO are tighter than those with stellar mass (average 1-σ scatter ∼ 0.09 compared to the average 1-σ scatter ∼ 0.12 with stellar mass).4.6.Higher-order stellar kinematics and orbital components van de Sande et al. (2017a) used the higher-order stellar kinematic moments (h 3 and h 4 ) to classify galaxies in the SAMI Galaxy Survey into 5 distinct classes based on each galaxy's individual h 3 versus V /σ signature.Galaxies belonging to Class 1 are typically the most massive, large and red.Most of Class 1 galaxies are also classified as slow rotators, indicating that they have more complex dynamical structures as compared to fast rotators.Galaxies in Class 2-5 are all consistent with being oblate rotating axisymmetric spheroids as based on λ Re and ε, but have a range of higher-order kinematic signatures.Galaxies in Class 2 are less massive, but still red, and reside in between slow and fast rotators.True fast rotators are in Class 3 and 4, with galaxies showing a strong anti-correlation of V /σ and h 3 .Galaxies in Class 5 have very high V /σ and ellipticity, but they do not show any anti-correlation with h 3 .Here, we examine the connection between the distributions of these classes and the orbital components of the galaxies in our sample.
In Fig. 14 we show the overall orbit circularity distribution for all the galaxies in our sample, grouped by their kinematic classes.The orbit circularity distribution is calculated by integrating the probability distribution of λ z over all radii within R max,h3h4 and normalizing it to unity, similarly to Fig. 12. R max,h3h4 is the radius within which the h 3 versus V /σ signatures were derived for each galaxy, due to S/N restrictions (van de Sande et al. 2017a).Within each subpanel in Fig. 14, we have ordered the galaxies by their intrinsic λ Re,EO values.The color indicates the normalized density, ω, of the orbits on the phase space.There is a clear distinction between the orbital distributions, depending on the galaxy kinematic class.In general, hot orbits are more dominant in galaxies belonging to Class 1, and they decrease going towards Class 5, with Class 4 showing the lowest values.The contribution of cold orbits becomes There is a clear increase of hot orbits (red diamonds) with increasing stellar mass (and B/T ratio), while the fraction of warm orbits (orange circles) decreases with increasing stellar mass (and B/T ratio), both of them showing a large scatter.Hot orbits decrease with increasing λRe,EO (and εintr), while the fraction of warm orbit increases with increasing λRe,EO (and εintr).The fraction of cold orbits (blue triangles) is also declining towards more massive galaxies and increases with galaxies becoming flatter.The fraction of counter-rotating orbits (green squares) does not show any significant trend with B/T ratio or εintr, but it does decrease with increasing λRe,EO.The correlation between the orbital fractions and λRe,EO shows very little scatter.
more important in Classes 3, 4 and 5, while warm orbits can also be a significant fraction for galaxies in Class 2. Counter-rotating orbits do not have any significant contribution for Class 3 and 5.The distribution of orbits in each class is clearer if we look at their integrated distributions, shown in Fig. 15.Within each class, there are also clear trends of the orbital components with λ Re , so that, as expected, cold orbits are increasing with increasing λ Re (rotationally supported galaxies).Similarly, warm orbits also increase with increasing λ Re,EO .In contrast, the hot component becomes less important with increasing λ Re,EO , while the counter-rotating orbits do not show any par-ticular trend.In particular, in slow-rotating galaxies, the main contribution is given by hot orbits.This is not unexpected, since these galaxies are expected to be pressure-supported.The warm component starts to become important for galaxies in Class 2, with its contribution increasing with increasing λ Re,EO .Galaxies in Class 3, 4 and 5 show higher contributions from warm and cold orbits for all the galaxies (compared to Class 1 and 2).We do not find strong evidence for a difference in the orbital distribution between the higher-order kinematic Classes 3-5 as derived from the circularity diagram.Nonetheless, the existence of the different signatures in the higher-order moment maps points to kine-matic features that are not captured in the λ z -r space.This will be explored further in future work, but is beyond the scope of this paper.

DISCUSSION
We have constructed Schwarzschild orbit-superposition models of 161 passive galaxies from the SAMI Galaxy Survey in order to derive intrinsic properties such as the internal mass distribution, intrinsic stellar shape, velocity anisotropy and orbit circularity distribution.We find that changes in the internal structures are mostly correlated with the stellar mass of the galaxies.

Fractions of dark matter
We find an average value of the dark matter fraction of f DM = 0.28, with a standard deviation of 0.20, within 1R e .In general, our results for f DM are broadly consistent with previous stellar dynamic determinations within 1R e found in the literature which also all assume a NFW dark matter halo distribution (Fig. 16).For example, Gerhard et al. (2001) found f DM = 0.1−0.4 from spherical dynamical modelling of 21 ETGs, Cappellari et al. (2006) inferred a median f DM ≈ 0.3 by comparing dynamics and population masses of 25 ETGs, and assuming a universal IMF, Thomas et al. (2007Thomas et al. ( , 2011) ) measured f DM = 0.23 ± 0.17 via axisymmetric dynamical models of 17 ETGs, Cappellari et al. (2013) measured a f DM of 0.15 for early-type galaxies in ATLAS 3D using Jeans Anisotropic Modelling (JAM), with galaxies showing an increasing fraction of dark matter with increasing mass for masses log(M /M ) > 10.6, consistent with our findings here.Similar results were also found by Posacki et al. (2015), f DM = 0.14 for 55 earlytype galaxies from stellar dynamics and lensing, and by Poci et al. (2017) -f DM = 0.19 using JAM to model a sample of 258 early-type galaxies in ATLAS 3D .For the Milky Way, Bland-Hawthorn & Gerhard (2016) found a f DM = 0.3.Overall, these studies show that baryons dominate the centers of galaxies, especially in our mass range, where the efficiency of galaxy building is peaking.Jin et al. (2020) found a similar trend for early-type galaxies in the MaNGA sample, with the f DM for the most massive galaxies (11.0 < log(M /M ) < 11.5) generally above 0.4, similar to what we see for galaxies in the same mass bin.However, we note that, as suggested by model tests with mock data from the Illustris simulations (Jin et al. 2019), estimations of f DM can have a systematic offset as a result of modelling the dark matter halos assuming that galaxies follow a NFW profile, which may not be correct.This is an interesting aspect that will need to be explored further and tested with a range of simulations and datasets.The trend we see in the f DM with stellar mass is also consistent with predictions from simulations, where galaxies with log(M /M ) ∼ 10.6 are the most efficient at forming stars (e.g.Behroozi et al. 2010Behroozi et al. , 2013;;Henriques et al. 2019).The physical interpretation of this behavior is the interplay between the feedback processes that impact star formation efficiency at different galaxy masses.Supernova feedback is more effective at reheating and expelling gas in low-mass galaxies, while AGN feedback is more effective in high-mass galaxies.

Intrinsic shape distribution
As seen in Fig. 8, the majority of our galaxies are very close to oblate axisymmetric (73% ± 3%), with T Re ≤ 0.1 , with varying degrees of intrinsic flattening, with 19% ± 3% being mildly triaxial (0.1 < T Re ≥ 0.3) and a small percentage (8% ± 2%) being triaxial/prolate (T Re > 0.3).There is a weak increase in the triaxiality parameter with increasing stellar mass.The percentage of galaxies that are non-oblate (T Re > 0.1) increases with increasing stellar mass, going from 12% ± 4% below 10 10.50 M /M to 29% ± 2% above this mass.
Triaxial Schwarzschild orbit-superposition dynamical models allow to measure intrinsic shapes directly.Previous studies used statistical methods to derive intrinsic shape properties; for example Kimm & Yi (2007) studied a sample of 3922 galaxies from SDSS (Adelman-McCarthy et al. 2006) and found that more massive galaxies are more likely to be triaxial than lower-mass galaxies.Foster et al. (2017) derived the intrinsic shape of 845 galaxies in the SAMI Galaxy Survey using an algorithm to simultaneously invert the distributions of apparent ellipticities and kinematic misalignments using the methodology of Weijmans et al. (2014).They find the majority (∼ 85%) of the galaxies in their sample to be oblate axisymmetric, in good agreement with Weijmans et al. (2014) and our results.Our result is also in agreement with previous results from the Illustris simulations, where only a very small fraction of galaxies are found to have prolate shapes, with the fraction decreasing to zero prolate galaxies below log(M /M ) = 11.48 (Li et al. 2018b).Jin et al. (2020) found higher fractions of triaxial and prolate galaxies in a sample of 149 early-type galaxies from the MaNGA survey.This discrepancy is partly explained by their higher stellar mass range analysed (their stellar masses ranged between 10 9.9 and 10 11.8 M ), and their different sample selection.Jin et al. (2020) also find an increase of the fraction of non-oblate galaxies with increasing stellar mass, in agreement with our results.The color indicates the normalized density, ω, of the orbits on the phase space.Galaxies in Class 1 are dominated by hot orbits.Warm orbits become important for galaxies in Class 2, in particular at higher values of λRe,EO, with the warm orbits contribution increasing for Classes 3, 4 and 5. Hot orbits become less important with increasing λRe,EO.

Velocity Anisotropy
We find that galaxies with higher ellipticities have, in general, more negative values of β r .This means that flatter galaxies are more tangentially anisotropic than rounder galaxies, while the latter are more likely to be supported by radial anisotropy.Moreover, we find a tight relationship of β r with λ Re,EO .This is not unexpected, since both parameters are a measure of rotation.The idea that the most giant early-type galaxies are not flattened by rotation but by anisotropy was proposed in the late 1970s (Bertola & Capaccioli 1977;Illingworth 1977;Binney 1978), however most of the dynamical modelling methods available to date do not allow for triaxiality, which is needed for a significant fraction of massive galaxies in order to construct accurate models.
Our results are also in agreement with more recent studies.For example, Gerhard et al. (2001) found that most of the galaxies in their sample of 21 ETGs were moderately radially anisotropic (β r ≈ 0.2), in agreement with the values we find in this study.

Orbital structures
We find that the hot orbital component generally dominates within R e , becoming the most prevalent component among galaxies with total stellar mass log(M /M ) > 11.As expected, bulge-dominated galaxies have high fractions of hot orbits (consistent with a pressure-supported bulge).In most galaxies a substantial number of stars within R e are on warm orbits, with the contribution becoming more important at lower stellar masses.The cold component rarely dominates within R e and its importance decreases with increasing stellar mass.The counter-rotating component is roughly constant for galaxies at all masses.
Stellar orbit distributions have only been derived explicitly before for two large (N>100) samples of galaxies, in the CALIFA (Zhu et al. 2018c) and MaNGA (Jin et al. 2020) surveys.We show the orbital fractions derived for early-type CALIFA and MaNGA galaxies, as well as the results from this work, in Fig. 17.The variations of the fraction of orbits is in good agreement with the general trends with stellar mass seen by Zhu et al. (2018c) and Jin et al. (2020).Jin et al. (2020) also found an increase in the fraction of hot orbits for massive (log(M /M ) > 11) galaxies, similar to what we find.
Previous studies that did not have access to stellar orbit modelling, commonly used the proxy for the spin parameter λ Re , and the flattening of galaxies, to shed light on galaxy intrinsic properties.Schwarzschild dynamical models allow us to explain the trends in λ Re by showing the contributions from different orbital components, providing a new insight into how λ Re is built-up.We measured the edge-on λ Re,EO from our model fits and compared it to the orbital fractions, shown in Fig. 14.We find a clear trend of the fractions of orbits with λ Re,EO : hot orbits show a rapid decrease in fraction with increasing λ Re,EO , while warm orbits have the opposite behaviour (increasing rapidly with increasing λ Re,EO ).
Counter-rotating orbits have slightly lower fractions for galaxies with higher spin parameter, while cold orbits show low fractions up to λ Re,EO ≈ 0.3, after which their importance starts to increase.This confirms that λ Re,EO is a good indicator of the underlying orbit distribution of a galaxy.The observed spin parameter used in the literature (e.g.Cappellari et al. 2007;Emsellem et al. 2007Emsellem et al. , 2011;;van de Sande et al. 2017a) is a projected quantity along an often-unknown line-of-sight viewing angle.Slow rotators are found to be more massive, dominating above 2 × 10 11 M (e.g.Emsellem et al. 2011;Cappellari 2016;Brough et al. 2017;Veale et al. 2017;Greene et al. 2017;van de Sande et al. 2017bvan de Sande et al. , 2021b)).This is in agreement with our more direct orbit-based finding of an increase of the hot component with increasing galaxy mass and the hot component starting to be dominating for galaxies with log(M /M ) > 10.75.

Implications for galaxy formation
While the degeneracy due to deprojection impacts the reliability of the recovered shape (Rybicki 1987;Krajnović et al. 2005;de Nicola et al. 2020), the Schwarzschild orbit-superposition method is still the best method that exists to derive the true threedimensional structure of individual galaxies.In this paper we find that the changes of internal structures within 1R e are correlated with the total stellar mass of individual galaxies.In particular, we find a rapid change in structure for galaxies above a stellar mass log(M /M ) ∼ 11.Below this stellar mass, galaxies tend to be oblate and with a substantial number of stars within R e on warm orbits, while higher-mass galaxies with log(M /M ) > 11 tend to be more triaxial and dominated by hot orbits.A similar change is also seen in the fraction of dark matter (Fig 7).The change in the hot and warm orbital fractions that we observe in Fig. 13 at stellar masses higher than ∼ 10 11 M and the change in intrinsic shape at similar mass that we see in Fig. 8 could be interpreted as an indication of different formation channels.In particular, major and minor mergers are found to be the main driver of triaxial and prolate shapes, while exclusively very minor mergers are largely associated with triaxial systems and oblate slow rotators are formed in the absence of mergers (Lagos et al. 2020).The increasing fractions of hot orbits with increasing stellar mass supports a scenario where the most massive slow rotators form via gas-poor major mergers (Li et al. 2018b).
The trends we observe in the inner parts of passive galaxies (within 1R e ) are generally consistent with the two formation paths of early-type galaxies (e.g.see recent review by Cappellari 2016).In this picture slowrotating ETGs assemble near the center of massive dark matter halos via intense star formation at high redshift, and their evolution is dominated by gas-poor mergers.
These galaxies are more likely to be triaxial and more massive, in agreement with what we find.By comparison, low-mass fast-rotating ETGs grow via gas accretion and their structures show similarities with that of spiral galaxies.Moreover, since the warm component can be interpreted as being similar to a thick disk, the increasing contribution that we see from warm orbits in fastrotating galaxies provides further evidence for disk-like components in these systems as indicated by Krajnović et al. (2008).
Simulations suggest that stars on different orbits have different formation paths.The cold components are mostly young stars formed in-situ, the warm component likely traces old stars formed in-situ, or stars being heated from cold disks via secular evolution, and a small fraction of the warm component stars could be accreted (Gómez et al. 2017;Park et al. 2021).The stars on hot orbits in the outer regions should mostly be accreted (Gómez et al. 2017;Tissera et al. 2017) via minor or major mergers, while stars on hot orbits in the inner regions are predicted to have formed at high-redshift.Further comparison with simulations will help us to understand the physical processes that lead to the orbit distribution observed at present times.

Evidence of early accretion from stellar populations
Resolved stellar dynamics trace the change in angular momentum and orbital distribution of stars due to mergers, but major mergers are likely to have obscured the effects of earlier interactions.However, evidence of these earlier interactions can be found in the stellar populations.In particular, a galaxy's mean stellar age provides information on when the stars were formed (e.g.Tinsley 1980;Bender et al. 1993;Park et al. 2021).So combining stellar population and stellar kinematic studies can provide unique but complementary insights into how galaxies build-up their stellar mass and angular momentum.
van de Sande et al. ( 2018) studied a sample of galaxies in the SAMI Galaxy Survey and found that there is a strong relation between V /σ Re and mean stellar age, such that galaxies with young stellar populations are predominantly rotationally supported, whereas galaxies with old stellar populations are more pressure supported by random orbital motion of stars.For the large majority of galaxies that are oblate-rotating spheroids, they found that characteristic stellar age is related to the intrinsic ellipticity of galaxies.They studied a full range of morphologies, but showed that this trend is still observed when galaxies are in early-type or late-type subsamples.
To check whether this relation holds for our parameters derived using Schwarzschild models, we color-code our data in the λ Re,EO − ε intr,EO plot by luminosityweighted, mean stellar population age (see Scott et al. 2017) in Fig. 18 and use LOESS smoothing to recover any mean underlying trend.We find a good match to the trends as found by van de Sande et al. ( 2018), with slow-rotating galaxies being generally older and rounder than fast-rotating galaxies.This relationship is consistent with predictions from hydrodynamical cosmological simulations and observations, where slow-rotating galaxies form via intense star formation at high redshift, and evolve from a set of processes dominated by gas-poor mergers (Cappellari 2016).
All the results presented here are in agreement with a formation scenario in which passive galaxies form through two main channels and where the changes of internal structures within 1R e are generally correlated with the total stellar mass of the individual galaxies.Re, EO 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 log Age [Gyr] Figure 18.λRe,EO as a function of the intrinsic ellipticity εintr,EO, calculated at 1Re.The magenta line represents the relation between the anisotropy parameter βz and the intrinsic ellipticity εintr,EO, for galaxies viewed edge-on, that bounds all regular rotating galaxies.Galaxies are colorcoded by log Age, and LOESS smoothed to recover any mean underlying trend.Older galaxies are generally slow-rotating and rounder than younger systems.

CONCLUSION
We constructed Schwarzschild orbit-superposition models of 161 passive galaxies, from the SAMI Galaxy Survey, with stellar masses raging from 9.5 < log(M /M ) < 11.4.We derived intrinsic properties such as the internal mass distribution (for both stellar and dark matter), intrinsic stellar shape (axis ratios and ellipticity), velocity anisotropy and orbit circularity distribution which gives us the most detailed insight into their assembly history.We draw the following conclusions: • Passive galaxies have an average dark matter fraction f DM = 0.28 ± 0.20, consistent with previous results (Fig. 16).
The fraction of non-oblate galaxies increases with increasing stellar mass, with a sudden change at ∼ 10 10.50 M /M (Fig. 8).
• Galaxies with high intrinsic ellipticity (flat galaxies, ε > 0.7) are found to be more isotropic (β r ∼ 0) or more tangentially anisotropic (β r < 0; Fig. 9).β r is anti-correlated with the spin parameter λ Re,EO , so that β r decreases with increasing λ Re,EO , consistent with slow-rotating galaxies being more radially anisotropic and fast-rotating galaxies being more tangentially anisotropic (Fig. 11).
• By dividing the stellar orbital distribution into cold, warm, hot, and counter-rotating components, we find that the hot component generally dominates within R e , becoming the most prevalent component among galaxies with total stellar mass log(M /M ) > 11.In most galaxies a substantial number (∼ 40% of stars within R e ) are on warm orbits, with the warm contribution becoming more important at lower stellar masses.The contribution from the cold orbital components is small across stellar mass, with its contribution decreasing further with increasing mass.The counterrotating component is roughly constant for galaxies at all masses (Fig. 13).
• The changes of internal structures (fraction of dark matter, f DM , intrinsic shape and orbital distribution) within 1R e are correlated with the total stellar mass of the individual galaxies.
• The fractions of orbits show tight correlations with the intrinsic λ Re,EO , with hot orbits being dominant for slow-rotating galaxies and contributions from warm and cold orbits becoming more important with increasing λ Re,EO .We also find a clear distinction between the orbital distributions of galaxies, depending on their kinematic class (from van de Sande et al. (2017a) based on the higher-order V /σ -h 3 signatures).Class 1 is dominated by hot orbits, with little contribution from other components.The contribution of warm orbits increases from Class 2 to 5, while the contribution from hot orbits become less important.Class 4 and 5 also show contributions from cold and counter rotating components (Fig. 14 and 15).
• These results are in agreement with a formation scenario in which galaxies form through two main different channels.Slow-rotating ETGs assemble near the center of massive dark matter halos via intense star formation at high redshift, and their evolution is dominated by gas-poor mergers.These galaxies are more likely to be triaxial and more massive, dominated by radial anisotropy, in agreement with what we find.By comparison, lowmass fast-rotating ETGs grow via gas accretion and their structures show similarities with that of spiral galaxies.Moreover, the intrinsic shapes of slow rotators could point to different type of mergers in their evolutionary history.
The work presented here expands on previous analyses by including the higher-order stellar kinematic moments.We found that including the higher-order kinematic moments h 3 and h 4 can improve the model fits, even if the h 3 and h 4 measurements have high uncertainties.We therefore recommend the inclusion of h 3 and h 4 in future works.Moreover, since h 3 and h 4 are quantities that are predicted to be connected with a galaxy's assembly history (Naab et al. 2014), studying their relation to the internal orbital structure of galaxies provides an extra tool to help disentangle the different possible formation scenarios.We did not find a significant difference between the orbital components of fast-rotating galaxies in the different high-order classes from van de Sande et al. (2017a) (determined using the V/sigma-h3 correlation) in the λ z − r space (Fig. 14), but this is an interesting aspect that should be further explored in future works with larger samples (e.g. the forthcoming Hector survey; Bryant et al. 2016).

ACKNOWLEDGEMENTS
We thank the anonymous referee for their comments that helped to improve this manuscript.We thank the DYNAMITE team for their invaluable support with the code and optimisation and Michele Cappellari for useful discussions.The SAMI Galaxy Survey is based on observations made at the Anglo-Australian Telescope.Software: pPXF (Cappellari & Emsellem 2004;Cappellari 2017), MgeFit (Cappellari 2002), Voronoi binning where f orb T OT is the orbital fraction for cold, warm, hot or counter-rotating (CR) -derived from the total map and f orbmap is the orbital fraction derived for one of the 4 kinematic maps -R max = 0.5R e , 1R e , 1.5R e and 2R e .Each point in Fig. 22 shows the average of the 4 residuals (one for each orbital component), color-coded by the R max of the maps.We also show the residuals for the fraction of dark matter within 1R (f DM ), the mass-to-light ratio in the r-band (M/L r ) and the intrinsic axis ratios at 1R e -p Re and q Re in Fig. 23.The average residuals for each of the maps is shown in Table 1.
Comparing the derived values within 1R e of the different maps for each galaxy, we find a general good agreement for all input R max maps, with the exception of those retrieved from the R max = 0.5R e maps, which show a large scatter.We are therefore confident in the values estimated within 1R e calculated using maps that extend to at least 1R e for the analysis presented here.1. Average residuals between the derived orbital fractions, fraction of dark matter within 1Re, mass-to-light ratio in the r-band (M /Lr) and intrinsic axis ratios pRe and qRe of each of the 4 maps from the values derived from the total maps for the galaxies in our CALIFA test sample.Comparing the derived values within 1Re of the different maps for each galaxy, we find a general good agreement for all input Rmax maps, with the exception of those retrieved from the Rmax = 0.5Re maps, which have larger average residual.
B. EXAMPLE GALAXIES 9403800123, 9011900793, 220465 AND 9008500323 The parameter space for the complete model runs for example galaxies 9403800123, 9011900793, 220465 and 9008500323 (Fig. 3,4,5 and 6) are shown in Fig. 24, Fig. 25, Fig. 26 and in Fig. 27, respectively.The dots represent the parameters we have explored.Models within the best-fit region are color-coded according to their χ 2 values.The largest red dot highlighted with a black cross indicates the best-fit model.Fig. 28 to Fig. 31 show the obtained internal mass distribution, orbit circularity, triaxiality and tangential anisotropy for the four galaxies.
To test whether the parameter grid is well sampled in our iterative grid search, we run a super-sampled grid search for the 4 example galaxies.The best-fit parameters for both the default (∼1250 models) and the "super-sampled" (> 6000) grid search are consistent with one another within the 1-σ confidence level, converging on global minima.
We also tested whether including h 5 and h 6 make a significant difference to our best-fit model for the example galaxies.Fixing h 5 and h 6 to 0 and allowing the model to fit these higher moments does not significantly improve the fit.The variations in h 5 and h 6 are quite small (∼ 0.06) and there are no significant changes in the kinematic fit, nor in the χ 2 (derived from the fit to the measured moments) level (for example χ 2 red changed from 2.22 to 2.18 for example galaxy 9403800123) or morphology.

C. UNCERTAINTIES ON THE MODEL BEST-FIT PARAMETERS
In addition to the 1-σ fluctuations from the best-fit model, we use Monte Carlo realizations to estimate the uncertainties on our best-fit values.This is particularly important to derive the uncertainties in the underlying model properties which are not accessible from the 1-sigma confidence level directly.This approach factors in convergence issues and grid sampling, with no assumptions about how the model parameter uncertainties are distributed -only that the kinematic errors are Gaussian (a common assumption).To this end, we select 16 SAMI galaxies (∼ 10% of the total sample), spanning different regions in the size -stellar mass plane.We apply Monte Carlo realizations, as described below, to each one of them, and we use the resulting variations from the best-fit parameters as the uncertainties for galaxies located in similar locations of the galaxy mass-size plane.For each galaxy, we take the kinematic  values from the best-fit model and perturb them by adding noise, taken from a Gaussian distribution with standard deviation equal to the mean error of each observed kinematic moment (V, σ, h 3 , h 4 ).We keep the standard deviation as the uncertainty for each perturbed value.We tested repeating this process to have 30, 50 and 100 different realizations.
We then derive the best fit for each of the perturbed kinematic maps, using the same iterative grid search described in Sec.3.4.We compare the orbital weights retrieved from each realization and we find that there is in general good We find that the fraction of the orbits in passive galaxies follow a unimodal distribution.This becomes more evident when considering 50 or 100 realizations (right-hand panels of Fig. 32).We  for β r and ∼ 15% for λ Re,EO .The uncertainties for the fraction of dark matter are around 8%.To the uncertainty of each parameter derived with this method we also add, in quadrature, the 1σ confidence level from the parameter grid, which represents the model fluctuations.This method is applied to derive the uncertainties of all the quantities presented in this work.(2007).Galaxies with higher ellipticities have higher values of βz.This means that flatter galaxies are more anisotropic than rounder galaxies.

E. RATIO OF ORDERED TO RANDOM MOTION
For completeness, we also measure the ratio of ordered to random motion V /σ, also measured within 1R e , using the definition from Cappellari et al. (2007): V /σ as a function of the ellipticity εintr,EO derived from MGE fits to the edge-on projected maps, calculated at 1Re.The magenta line corresponds to the relation βz = 0.7ε for edge-on galaxies (Emsellem et al. 2007).Data points are color-coded by the velocity anisotropy βr.As expected, V /σ increases with increasing intrinsic ellipticity.

Figure 1 .
Figure1.Number of Voronoi bins within 1Re that meet our quality criteria versus the maximum radius available for stellar kinematics (in units of Re) for the galaxies in the SAMI Galaxy Survey (1589 galaxies; grey circles) and in the CALIFA survey (259 galaxies; violet diamonds).Black dashed lines indicate Rmax/Re = 1 and Voronoi bins = 85.We calculate the marginalized fractions of galaxies to the total number in each sample, by mass and size, and show them in the top and left panels of the figure.Grey lines are for SAMI galaxies, while the violet lines are for CALIFA galaxies.The CALIFA and the SAMI samples have similar distributions in Voronoi bins and radial coverage, although there are more CALIFA galaxies with measurements up to 2Re.For this analysis we select galaxies in the top right corner (Rmax > 1Re and Voronoi bins > 85).
Figure2.Effective radius, Re, versus stellar mass.Blue circles are the passive galaxies in the SAMI sample with log 10 (M /M ) > 9.5 and Re > 2 (738), orange squares are the galaxies included in the final sample (161).We calculate the marginalized fractions of galaxies with respect to the total number in each sample, by mass and size, and show them in the top and left panels of the figure.Blue lines are for the passive galaxies in the SAMI Galaxy Survey, while the orange lines are for our final sample.The two samples are slightly different in the marginalized mass and size distributions, so that we have higher fractions of massive and large galaxies in the final sample compared to the initial sample.This is due to selecting galaxies with more than 85 Voronoi bins.

Figure 3 .
Figure3.Example of a galaxy with excellent spatial sampling: SAMI CATID 9403800123 in the cluster Abell 4038.This galaxy (log M /M = 11.05 and Re = 5.52 ) is a non edge-on oblate galaxy and has stellar kinematic measurements up to 1.36 Re and counts 255 spatial bins within 1Re (black ellipse).Columns show 2D maps for, from left to right, flux, velocity, velocity dispersion, h3 and h4.First row shows the observed maps, second row shows the best-fit maps derived from the Schwarzschild modelling and the third row shows the residuals, calculated as the difference between the observation and the model, divided by the observational uncertainties.The best-fit model maps (χ 2 red = 2.22) accurately reconstruct the structures seen in the observations, not only for the velocity and velocity dispersion maps, but also for h3 and h4.

Figure 4 .
Figure 4. Example of a galaxy near the minimum requirement of 85 spatial bins: SAMI CATID 9011900793 in the cluster Abell 119.This galaxy (log M /M = 10.34 and Re = 5.19 ) is a edge-on oblate galaxy and has stellar kinematic measurements up to 1.45 Re and 87 spatial bins within 1Re.Panels are as in Fig. 3.The best-fit model maps (χ 2 red = 1.72) accurately reconstruct the structures seen in the observations.

Figure 5 .
Figure 5. Example galaxy SAMI CATID 220465, in the GAMA region.This galaxy (log M /M = 11.31 and Re = 5.00 ) is a triaxial galaxy and has stellar kinematic measurements up to 1.5 Re and 142 spatial bins within 1Re.Panels are as in Fig. 3.The best-fit model maps (χ 2 red = 1.79) accurately reconstruct the structures seen in the observations, not only for the velocity and velocity dispersion maps, but also for h3 and h4.

Figure 6 .
Figure 6.Example galaxy SAMI CATID 9008500323, in the cluster Abell 85.This galaxy (log M /M = 10.78 and Re = 4.15 ) is a prolate galaxy and has stellar kinematic measurements up to 1.81 Re and 104 spatial bins within 1Re.Panels are as in Fig. 3.The best-fit model maps (χ 2 red = 1.99) accurately reconstruct the structures seen in the observations, not only for the velocity and velocity dispersion maps, but also for h3 and h4.

Figure 7 .
Figure7.Fraction of dark matter to total mass (fDM = M dark /Mtot) within 1Re as a function of stellar mass.The median values of the fraction of dark matter for each of the 4 mass bins are shown as dark blue squares, with the error bars marking the 25 th and 75 th percentiles.The blue solid line is the parabolic best-fit to the data -fDM ∼ 0.10 + 0.17 × (log M M − 10.59) 2 .The fit is calculated by taking into account the errors on the data points, so that it is critically determined by the data points with small errors.The shaded region represents the error on the best-fit.

Figure 9 .
Figure9.Velocity dispersion anisotropy in spherical coordinates βr within 1Re as a function of intrinsic ellipticity (εintr = 1 − q) at 1 Re, color-coded by stellar mass.The grey dashed line represents isotropy, βr = 0. Negative βr indicates tangentially-anisotropic systems (supported by rotation), while positive βr are radially-anisotropic systems (supported by random motions).Galaxies with very high ellipticity are close to isotropic or tangentially anisotropic.Radially-anisotropic galaxies are generally more massive.

rFigure 11 .
Figure10.λRe,EO as a function of the ellipticity εintr,EO derived from MGE fits to the edge-on projected maps, calculated at 1 Re.The magenta line corresponds to the relation βz = 0.65 × ε for edge-on galaxies as inCappellari et al. (2007).Galaxies are colored by their velocity anisotropy βr in panel a and LOESS smoothed in panel b.As expected, λRe,EO increases with increasing intrinsic ellipticity.Galaxies that are radially anisotropic show low-to mid-values of ellipticity and λRe, while galaxies with high ellipticity and λRe,EO are more isotropic or tangentially anisotropic.

Figure 13 .
Figure13.Fractions of orbital components as a function of: a) stellar mass, b) bulge to total flux ratio, B/T, c) λRe,EO, d) εintr.Bold points show the median values for each mass bin, with error bars representing the 1σ scatter around the median value.There is a clear increase of hot orbits (red diamonds) with increasing stellar mass (and B/T ratio), while the fraction of warm orbits (orange circles) decreases with increasing stellar mass (and B/T ratio), both of them showing a large scatter.Hot orbits decrease with increasing λRe,EO (and εintr), while the fraction of warm orbit increases with increasing λRe,EO (and εintr).The fraction of cold orbits (blue triangles) is also declining towards more massive galaxies and increases with galaxies becoming flatter.The fraction of counter-rotating orbits (green squares) does not show any significant trend with B/T ratio or εintr, but it does decrease with increasing λRe,EO.The correlation between the orbital fractions and λRe,EO shows very little scatter.

Figure 14 .
Figure14.Orbit circularity distribution calculated by integrating the probability distribution of λz over all radii within R max,h3h4 and normalizing it to unity, for all the galaxies in our sample, grouped by their kinematic classes from van de Sande et al. (2017a) based on the higher-order (V /σ -h3) signatures.Each class has been ordered by the intrinsic λRe,EO values.The color indicates the normalized density, ω, of the orbits on the phase space.Galaxies in Class 1 are dominated by hot orbits.Warm orbits become important for galaxies in Class 2, in particular at higher values of λRe,EO, with the warm orbits contribution increasing for Classes 3, 4 and 5. Hot orbits become less important with increasing λRe,EO.

Figure 15 .
Figure 15.Normalized density of orbits as a function of the orbit circularity λz for each kinematic class from van de Sande et al. (2017a) based on the higher-order (V /σ -h3) signatures.There is a decrease in the contribution from hot orbits going from Class 1 to 5, with Class 4 having the lowest value.Warm orbits become more important from Class 2 to Class 5. Counter-rotating orbits do not have any significant contribution for Class 3 and 5.

Figure 16 .
Figure 16.Median values of the fractions of dark matter (fDM = M dark /Mtot) within 1 Re as a function of stellar mass for: SAURON (green triangle, Cappellari et al. 2006), ATLAS 3D (green dotted line, Cappellari et al. 2013 -derived with a cosmologically-motivated NFW halo) galaxies, the Milky Way (red cross, Bland-Hawthorn & Gerhard 2016), CALIFA (purple stars, Zhu et al. 2018b), (MaNGA (black diamonds, Jin et al. 2020) and SAMI (dark blue squares and blue solid line).Horizontal error bars delimit the mass range covered by each study.Vertical error bars mark the 25 th and 75 th percentiles, when available.The shaded region represents the error on the best-fit line for SAMI galaxies.Our results are in good agreement with the results presented in the literature.

Figure 17 .
Figure17.Median values of the fractions of orbital components as a function of stellar mass.The cold component is shown in blue, the warm component in orange, hot component in red and counter rotating in green.SAMI passive galaxies are shown as filled points, MaNGA early-type galaxies as open points(Jin et al. 2020) and the shaded areas represent the median values of early-type galaxies in the CALIFA sample(Zhu et al. 2018c).Horizontal error bars delimit the mass range covered.Vertical error bars mark the 25 th and 75 th percentiles, when available.The distribution of the fractions of orbits in SAMI and MaNGA are similar.All three samples show similar trends of orbital fractions with stellar mass.

Figure 20 .Figure 21 .
Figure 20.Best-fit model for CALIFA test galaxy NGC5888 using 2-moments maps.Top: Observed luminosity, velocity and velocity dispersion Bottom: best-fit model luminosity, velocity and velocity dispersion.The model does not reproduce the velocity dispersion well.

Figure 23 .
Figure23.Average residuals between the derived fraction of dark matter within 1Re (fDM; top left), mass-to-light ratio in the r-band (M /Lr; top right) and intrinsic axis ratios pRe (bottom left) and qRe (bottom right) of each of the 4 maps from the values derived from the total maps for the galaxies in our test sample, as a function of Re.For each map, the residuals of the four orbital components are calculated following A1.Each point is color-coded by the value of Rmax of the map as shown in the bottom right corner.Comparing the derived values within 1Re of the different maps for each galaxy, we find a general good agreement for all input Rmax maps, with the exception of those retrieved from the Rmax = 0.5Re maps, which show a large scatter.

Figure 24 .
Figure24.Example galaxy 9403800123: model parameter grid.There are six free parameters: stellar mass-to-light ratio, M /Lr in solar units, the intrinsic shape of the flattest Gaussian component (pmin, qmin, umin), the dark matter halo concentration, log c, and dark matter fraction, logM200/M .The diamonds represent the parameters explored, with the best-fit model highlighted with a black cross.Models within the best-fit region are color-coded according to their χ 2 values shown in the color bar.The best-fit values are well constrained.

Figure 25 .
Figure 25.Example galaxy 9011900793: model parameter grid.The diamonds represent the parameters explored, with the best-fit model highlighted with a black cross.Models within the best-fit region are color-coded according to their χ 2 values shown in the color bar.

Figure 26 .
Figure 26.Example galaxy 220465: model parameter grid.The diamonds represent the parameters explored, with the best-fit model highlighted with a black cross.Models within the best-fit region are color-coded according to their χ 2 values shown in the color bar.

Figure 28 .
Figure 28.Example galaxy 9403800123 (top left panel), 9011900793 (top right panel), 220465 (bottom left panel) and 9008500323 (bottom right panel): enclosed mass.Cumulative total mass (in black), stellar mass (in red) and dark matter mass (in blue) as a function of the radius of the galaxy.Solid lines are the cumulative profiles calculated from the best-fit, while the filled regions indicate the errors.Grey dotted and dash-dotted lines are located at 1Re and at Rmax, respectively.At larger radii the dark matter contribution becomes more important.

Figure 32 .Figure 33 .
Figure 32.Distribution of the orbital weights for the Monte Carlo realizations around the best-fit model values found for example galaxy 91963.Left-hand plot: 30 realizations; central plot: 50 realizations; right-hand plot: 100 realizations.The dashed lines represent the best-fit values.The unimodal distributions of the orbital components become more evident when increasing the number of realizations.We use 50 Monte Carlo realizations to derive the uncertainties for our galaxies to optimize the model run-time required.

)
Figure34.V /σ as a function of the ellipticity εintr,EO derived from MGE fits to the edge-on projected maps, calculated at 1Re.The magenta line corresponds to the relation βz = 0.7ε for edge-on galaxies(Emsellem et al. 2007).Data points are color-coded by the velocity anisotropy βr.As expected, V /σ increases with increasing intrinsic ellipticity. ) The Sydney-AAO Multi-object Integral field spectrograph (SAMI) was developed jointly by the University of Sydney and the Australian Astronomical Observatory.The SAMI input catalogue is based on data taken from the Sloan Digital Sky Survey, the GAMA Survey and the VST ATLAS Survey.The SAMI Galaxy Survey is supported by the Australian Research Council center of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013, the Australian Research Council center of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020, and other participating institutions.The SAMI Galaxy Survey website is http://samisurvey.org/.
Figure22.Average residuals between the derived orbital fractions of each of the 4 maps from the values derived from the total maps for the galaxies in our test sample, as a function of Re.For each map, the residuals of the four orbital components are calculated following A1 and then averaged over the orbital components.Each point corresponds to the average value, color-coded by the value of Rmax of the map as shown in the bottom right corner.Comparing the derived values within 1Re of the different maps for each galaxy, we find a general good agreement for all input Rmax maps, with the exception of those retrieved from the Rmax = 0.5Re maps, which show a large scatter.