Jeans modeling of axisymmetric galaxies with multiple stellar populations
MMNRAS , 1–14 (2019) Preprint 19 February 2021 Compiled using MNRAS L A TEX style file v3.0
Jeans modeling of axisymmetric galaxies with multiple stellar populations
Caterina Caravita, , ★ Luca Ciotti and Silvia Pellegrini , Department of Physics and Astronomy, University of Bologna, via P. Gobetti 93/2, 40129 Bologna, Italy INAF-OAS of Bologna, via P. Gobetti 93/3, 40129 Bologna, Italy
19 February 2021
ABSTRACT
We present the theoretical framework and the numerical setting of JASMINE2, a code designed to efficiently solve theJeans equations for multi-component axisymmetric stellar systems. The models may include an arbitrary number of stellardistributions, a dark matter halo, and a central supermassive black hole; each stellar distribution is implicitly described by atwo-integral distribution function, and the stellar components can have different structural (density profile, flattening, mass,scale-length), dynamical (rotation, velocity dispersion anisotropy), and population (age, metallicity, initial mass function, mass-to-light ratio) properties. In order to determine the ordered rotational velocity and the azimuthal velocity dispersion fields ofeach component, we introduce a decomposition that can be used when the commonly adopted Satoh decomposition cannotbe applied. The numerical implementation of JASMINE2, and of the post-processing procedures (including projection), areoptimised to fully exploit the scalings allowed by the Poisson and the Jeans equations. For illustrative purposes, we presentthree multi-component galaxy models with a central black hole and a dark matter halo; one of the models is also used to testJASMINE2 against available analytical solutions.
Key words: galaxies: structure - galaxies: kinematics and dynamics - methods: analytical - methods: numerical
Axisymmetric galaxy models often represent an acceptable descrip-tion of real galaxies, beyond the zeroth-order approximation of spher-ical symmetry. Analytical models of one and multi-component ax-isymmetric galaxies are available (e.g., see Binney & Tremaine 2008,hereafter BT08, and references therein; see also Ciotti et al. 2021,hereafter CMPZ21), but these models, while important to highlightfundamental properties of the dynamics of axisymmetric systems,and to guide the construction of realistic galaxy models to be carriedout numerically, suffer from the restrictions imposed by the request ofanalytical tractability. From this point of view, analytical and numer-ical modeling should be seen as complementary approaches, each ofthem with their own merits and limitations.On the numerical side, the most common models are based on thesolution of the Jeans equations (e.g., Ciotti 2000, BT08). This ap-proach allows to model axisymmetric stellar systems, in the simplestassumption of a 2-integral phase-space distribution function (DF, e.g.Posacki et al. 2013), or a 3-integral DF (e.g. Cappellari 2008), start-ing from the assignment of the density components. As well known,the proper description of a stellar system should start from the as-signment of the phase-space DF of each separate component, andthe solution of the associated Poisson equation (Bertin 2014; Ciotti2000; BT08); however, in several applications, the Jeans approachis highly preferred, for its direct control on the density distributions(even if it leaves open fundamental issues such as the phase-spaceconsistency). ★ E-mail: [email protected]
Recently, we developed a numerical code for axisymmetric mod-eling of galaxies, based on the solution of the Jeans equations, namedJASMINE (
Jeans AxiSymmetric Models of galaxies IN Equilibrium ,Posacki et al. 2013). In this paper, we present a substantially improvedcode version, JASMINE2, complemented by a post-processing pro-cedure, especially designed to build (and project) multi-componentsystems, in a fast and flexible way. The new code allows for the choice,in input, of an arbitrary number of mass components: different stel-lar components, Dark Matter (DM) components, and a central BlackHole (BH). Then it computes numerically the gravitational potentialof each density component by using the well-known formula basedon complete elliptic integrals of the first kind (e.g., equation C.29in BT08). This numerical method is highly accurate but it can bequite time expensive, depending on the grid resolution and on thenumber of density components of the model. In order to efficientlyexplore the parameter space (that can be very large, especially formulti-component models), we fully exploit all the scalings allowedby the Poisson and the Jeans equations. The model construction isorganised in two logically distinct parts: in the first one, the
Potential& Jeans Solver , JASMINE2 computes the potential and then solvesthe Jeans equations for each scaled stellar density component; in thesecond, the
Post-Processing (PP) , the mass and luminosity weightsare assigned, and the kinematical decompositions imposed, and fi-nally the scaled solutions of the Jeans equations are combined, andthe resulting kinematical fields projected.This procedure allows to drastically reduce the computational timeneeded for the construction of a multi-component model: with a sin-gle run of the Potential & Jeans Solver, one can build a family ofgalaxy models, all characterised by the same set of scaled density © a r X i v : . [ a s t r o - ph . GA ] F e b C. Caravita, L. Ciotti and S. Pellegrini components; each specific model in the family is defined by choosingsuitable weights and kinematical decompositions in PP. In this waythe exploration of the parameter space is extremely fast and complete.Summarising, each stellar density component in a multi-componentmodel is characterised by different structural (density profile, flatten-ing, total mass, scale-length), dynamical (rotational support, velocitydispersion anisotropy), and stellar population (age, metallicity, Ini-tial Mass Function, mass-to-light ratio) properties. The addition of acentral BH and a DM halo is immediate.The paper is organised as follows. Section 2 presents the analyt-ical framework behind the numerical method, together with a newvelocity decomposition for the azimuthal velocity field, to be usedwhen the commonly adopted Satoh (1980) 𝑘 -decomposition cannotbe applied (a not uncommon case in multi-component systems). InSection 3 the scaling procedure is described in detail. In Section 4three illustrative galaxy models are built and few tests are mentioned.In Section 5 the main conclusions are summarised. The new JASMINE2 code has been considerably extended startingfrom the original JASMINE version (Posacki et al. 2013), to effi-ciently solve the Jeans equations for multi-component axisymmetricstellar systems, with a central BH, and embedded in a DM halo. Webegin by setting the general notation and the main theory behind thenumerical code: in this Section we focus on the multi-componentJeans equations, on a generalisation of the Satoh 𝑘 -decompositionfor azimuthal motions, and on the projections on the plane of the sky.In the following Section 3, we will describe the scaling procedure,and in particular how the scaled solutions of the Jeans equations areobtained (Potential & Jeans Solver), and then combined by adoptingsuitable weights (PP). JASMINE2 works in cylindrical coordinates ( 𝑅, 𝜑, 𝑧 ) , with the sym-metry axis of the models aligned with the 𝑧 -axis. In full generality,we consider models composed by 𝑁 different stellar density distri-butions 𝜌 ∗ 𝑖 ( 𝑅, 𝑧 ) , of total mass 𝑀 ∗ 𝑖 , so that the total stellar density 𝜌 ∗ and the total stellar mass 𝑀 ∗ of the system are given respectivelyby 𝜌 ∗ ( 𝑅, 𝑧 ) = ∑︁ 𝑖 𝜌 ∗ 𝑖 , 𝑀 ∗ = ∑︁ 𝑖 𝑀 ∗ 𝑖 , 𝑖 = , . . . , 𝑁. (1)From now on, sums over 𝑖 indicate sums over the 𝑁 stellar compo-nents. We assume that each 𝜌 ∗ 𝑖 is made of a Simple Stellar Population(see e.g. Renzini & Buzzoni 1986; Maraston 2005), i.e. by stars ofthe same age, chemical composition, Initial Mass Function, and inparticular the same mass-to-light ratio Υ ∗ 𝑖 . Therefore, the total stellardistribution 𝜌 ∗ can be considered a multiple stellar population; theluminosity density and the total luminosity of each stellar componentcan be written respectively as 𝜈 ∗ 𝑖 ( 𝑅, 𝑧 ) = 𝜌 ∗ 𝑖 Υ ∗ 𝑖 , 𝐿 𝑖 = 𝑀 ∗ 𝑖 Υ ∗ 𝑖 , (2)so that 𝜈 ∗ ( 𝑅, 𝑧 ) = ∑︁ 𝑖 𝜈 ∗ 𝑖 , 𝐿 = ∑︁ 𝑖 𝐿 𝑖 . (3) The local and average stellar mass-to-light ratios of the galaxy aregiven by Υ ∗ ( 𝑅, 𝑧 ) ≡ 𝜌 ∗ 𝜈 ∗ = (cid:205) 𝑖 𝜌 ∗ 𝑖 (cid:205) 𝑖 𝜌 ∗ 𝑖 / Υ ∗ 𝑖 , < Υ ∗ > ≡ 𝑀 ∗ 𝐿 = (cid:205) 𝑖 𝑀 ∗ 𝑖 (cid:205) 𝑖 𝑀 ∗ 𝑖 / Υ ∗ 𝑖 , (4)where it is apparent how in general the local stellar mass-to-lightratio in a multi-component model depends on position.From equation (1) the gravitational potential associated with thetotal stellar density is 𝜙 ∗ ( 𝑅, 𝑧 ) = ∑︁ 𝑖 𝜙 ∗ 𝑖 , (5)where 𝜙 ∗ 𝑖 ( 𝑅, 𝑧 ) is the potential originated by the density component 𝜌 ∗ 𝑖 . The presence of a central BH, of mass 𝑀 BH , produces thepotential 𝜙 BH ( 𝑟 ) = − 𝐺 𝑀 BH 𝑟 , 𝑟 = √︁ 𝑅 + 𝑧 , (6)and an axisymmetric DM halo, of density 𝜌 h ( 𝑅, 𝑧 ) and total mass 𝑀 h (when finite), produces the potential 𝜙 h ( 𝑅, 𝑧 ) . Therefore, in general,the total gravitational potential of the model is Φ ( 𝑅, 𝑧 ) = 𝜙 ∗ + 𝜙 h + 𝜙 BH = ∑︁ 𝑗 𝜙 𝑗 , 𝑗 = , . . . , 𝑁 + . (7)From now on, sums over 𝑗 indicate sums over all the 𝑁 + 𝑁 stellar components, the central BH, and theDM halo. In principle, also the DM distribution can be made ofdifferent components, with a trivial generalisation of the currentdiscussion, that is not necessary for the goal of this paper. JASMINE2fully exploits the linearity of the Jeans equations with respect to thestellar density (Section 2.2) and to the gravitational potential (Section3). In JASMINE2 we maintain the assumption that each stellar compo-nent is implicitly described by a 2-integral phase-space DF 𝑓 𝑖 ( 𝐸, 𝐽 𝑧 ) (in general different for each component), where 𝐸 and 𝐽 𝑧 are respec-tively the energy and the axial component of the angular momentumof each star (per unit mass) in the total potential Φ . Therefore, theDF of the total stellar distribution is the 2-integral function 𝑓 = ∑︁ 𝑖 𝑓 𝑖 . (8)As usual, we indicate with ( v 𝑅 , v 𝜑 , v 𝑧 ) the velocity componentsin the phase-space, and with a bar over a quantity the operationof average over the velocity-space. By construction, for each stellarcomponent v 𝑅𝑖 = v 𝑧𝑖 =
0, the only non-zero ordered velocity canoccur in the azimuthal direction 𝜐 𝜑𝑖 ≡ v 𝜑𝑖 , and finally for thevelocity dispersion tensor 𝜎 𝑅𝑖 = 𝜎 𝑧𝑖 ≡ 𝜎 𝑖 . Of course, from equation(8) similar relations hold for the kinematical fields of the total 𝜌 ∗ .The Jeans equations for each stellar component are obtained asvelocity averages of the Collisionless Boltzmann Equation over thecorresponding 𝑓 𝑖 (e.g. BT08), so that 𝜕𝜌 ∗ 𝑖 𝜎 𝑖 𝜕𝑧 = − 𝜌 ∗ 𝑖 𝜕 Φ 𝜕𝑧 ,𝜕𝜌 ∗ 𝑖 𝜎 𝑖 𝜕𝑅 = 𝜌 ∗ 𝑖 Δ 𝑖 𝑅 − 𝜌 ∗ 𝑖 𝜕 Φ 𝜕𝑅 , (9) MNRAS000
0, the only non-zero ordered velocity canoccur in the azimuthal direction 𝜐 𝜑𝑖 ≡ v 𝜑𝑖 , and finally for thevelocity dispersion tensor 𝜎 𝑅𝑖 = 𝜎 𝑧𝑖 ≡ 𝜎 𝑖 . Of course, from equation(8) similar relations hold for the kinematical fields of the total 𝜌 ∗ .The Jeans equations for each stellar component are obtained asvelocity averages of the Collisionless Boltzmann Equation over thecorresponding 𝑓 𝑖 (e.g. BT08), so that 𝜕𝜌 ∗ 𝑖 𝜎 𝑖 𝜕𝑧 = − 𝜌 ∗ 𝑖 𝜕 Φ 𝜕𝑧 ,𝜕𝜌 ∗ 𝑖 𝜎 𝑖 𝜕𝑅 = 𝜌 ∗ 𝑖 Δ 𝑖 𝑅 − 𝜌 ∗ 𝑖 𝜕 Φ 𝜕𝑅 , (9) MNRAS000 , 1–14 (2019) eans modeling with multiple stellar populations where Φ is the total potential in equation (7), and Δ 𝑖 ≡ v 𝜑𝑖 − 𝜎 𝑖 , 𝜎 𝜑𝑖 = v 𝜑𝑖 − 𝜐 𝜑𝑖 = Δ 𝑖 + 𝜎 𝑖 − 𝜐 𝜑𝑖 ; (10)therefore, in the isotropic case, Δ 𝑖 = 𝜐 𝜑𝑖 . Imposing the naturalboundary condition 𝜌 ∗ 𝑖 𝜎 𝑖 → 𝑧 → ∞ , the solution of equations(9) is 𝜌 ∗ 𝑖 𝜎 𝑖 = ∫ ∞ 𝑧 𝜌 ∗ 𝑖 𝜕 Φ 𝜕𝑧 (cid:48) 𝑑𝑧 (cid:48) , 𝜌 ∗ 𝑖 Δ 𝑖 = 𝑅 (cid:18) 𝜕𝜌 ∗ 𝑖 𝜎 𝑖 𝜕𝑅 + 𝜌 ∗ 𝑖 𝜕 Φ 𝜕𝑅 (cid:19) . (11)We notice that Δ 𝑖 can be also recast as a commutator-like integral(e.g., see equation 35 in CMPZ21), with some advantage for ana-lytical and numerical investigations; however we found by severalnumerical tests that Δ 𝑖 can also be accurately computed by (centred)numerical differentiation as in equation (11), so in JASMINE2 wemaintained this more direct way of evaluation.A central point for the working of JASMINE2 is the sum rule inthe phase-space imposed by the identity (8). In fact, with equation(8), we are assuming that the 𝑁 stellar components 𝜌 ∗ 𝑖 are physicallydistinct, each of them described by its own 𝑓 𝑖 , and so necessarily the 𝑁 pairs of equations (9) are the moment equations of each 𝑓 𝑖 in thetotal potential Φ . As usual, if 𝐹 ( x , v ) is a generic dynamical propertydefined over the phase-space, then 𝐹 𝑖 = ∫ 𝐹 𝑓 𝑖 𝑑 v 𝜌 ∗ 𝑖 , 𝐹 = (cid:205) 𝑖 𝜌 ∗ 𝑖 𝐹 𝑖 𝜌 ∗ , 𝐹 L = (cid:205) 𝑖 𝜈 ∗ 𝑖 𝐹 𝑖 𝜈 ∗ , (12)where the properties 𝐹 and 𝐹 L of the galaxy can be interpretedas the mass-weighted and the luminosity-weighted averages of the 𝐹 𝑖 , respectively. These considerations show how to combine thesolution for the single 𝜌 ∗ 𝑖 to obtain the dynamical fields associatedwith the total 𝜌 ∗ . Clearly, the Jeans equations for 𝜌 ∗ in equation (1)are obtained as the sum of equations (9) over the 𝑁 components, andtheir solution can be written as 𝜎 = (cid:205) 𝑖 𝜌 ∗ 𝑖 𝜎 𝑖 𝜌 ∗ , Δ = (cid:205) 𝑖 𝜌 ∗ 𝑖 Δ 𝑖 𝜌 ∗ , v 𝜑 = (cid:205) 𝑖 𝜌 ∗ 𝑖 v 𝜑𝑖 𝜌 ∗ , 𝜐 𝜑 = (cid:205) 𝑖 𝜌 ∗ 𝑖 𝜐 𝜑𝑖 𝜌 ∗ , (13)where the previous identities are of straightforward proof from equa-tions (12) and (10). Note that 𝜎 𝜑 is not given by the simple sumof the 𝜎 𝜑𝑖 of the single components, as 𝜎 in equation (13) above,because from equations (13) and (10) one has: 𝜎 𝜑 = v 𝜑 − 𝜐 𝜑 = Δ + 𝜎 − 𝜐 𝜑 = (cid:205) 𝑖 𝜌 ∗ 𝑖 ( 𝜎 𝜑𝑖 + 𝜐 𝜑𝑖 ) 𝜌 ∗ − 𝜐 𝜑 . (14)Similarly, from equation (12), we derive all the correspondingluminosity-weighted quantities; we do not give here their expres-sions, since they are just obtained by using as weights the luminositydensities 𝜈 ∗ 𝑖 of the components instead of the mass densities 𝜌 ∗ 𝑖 , inequations (13) and (14).Finally, the rotation curve in the equatorial plane is given in termsof the circular velocity 𝜐 c 𝑗 of each mass component as 𝜐 = ∑︁ 𝑗 𝜐 𝑗 . (15) As well known, equations (9) are degenerate in the azimuthal di-rection, i.e. they only provide v 𝜑𝑖 = 𝜎 𝜑𝑖 + 𝜐 𝜑𝑖 . The most commonphenomenological approach to break this degeneracy is the Satoh (1980) 𝑘 -decomposition. If Δ 𝑖 ≥ 𝜐 𝜑𝑖 = 𝑘 𝑖 √︁ Δ 𝑖 , 𝜎 𝜑𝑖 = 𝜎 𝑖 + (cid:0) − 𝑘 𝑖 (cid:1) Δ 𝑖 , (16)with negative values of 𝑘 𝑖 describing clockwise rotation. The spe-cial case 𝑘 𝑖 = 𝜎 𝜑𝑖 = 𝜎 𝑖 ),with flattening totally supported by rotation; while, if 𝑘 𝑖 =
0, thereis no net rotation ( 𝜐 𝜑𝑖 = 𝑘 𝑖 ( 𝑅, 𝑧 ) , also allowing for values greater than unity, up to a position-dependent maximum determined by the request 𝜎 𝜑𝑖 = 𝑘 𝑖 , so that from equations (16) and (13) the total 𝜌 ∗ willhave an effective Satoh parameter 𝑘 e given by: 𝑘 e ≡ 𝜐 𝜑 √ Δ = (cid:205) 𝑖 𝑘 𝑖 𝜌 ∗ 𝑖 √ Δ 𝑖 𝜌 ∗ √ Δ , 𝜎 𝜑 = 𝜎 + (cid:0) − 𝑘 (cid:1) Δ , (17)where 𝑘 e in general depends on position, even if the 𝑘 𝑖 do not.Clearly, in case of Δ 𝑖 < 𝜌 ∗ 𝑖 , the Satoh decomposition inequation (16) cannot be applied. The case of a negative Δ 𝑖 over someregions of space (or everywhere) is not frequently encountered in ap-plications, but it is not impossible; for example, it necessarily occursfor density components in multi-component systems with sphericallysymmetric total density, or in density distributions elongated alongthe symmetry axis (see e.g. Chapter 13 in Ciotti 2021). Indeed, Δ = Δ 𝑖 must be negative (excludingthe trivial case of all the subcomponents spherically symmetric, sothat Δ 𝑖 = Δ 𝑖 < 𝑓 𝑖 < 𝜑𝑖 = Δ 𝑖 + 𝜎 𝑖 < certainly the whole model must be discarded as unphysical, even if the solu-tion for the total stellar distribution are well-behaved. Therefore, iffor some component Δ 𝑖 <
0, then JASMINE2 checks the condition Δ 𝑖 + 𝜎 𝑖 ≥
0, and in case of positivity the Satoh decomposition isgeneralised to 𝜐 𝜑𝑖 = 𝑘 𝑖 √︃ Δ 𝑖 + 𝜎 𝑖 , 𝜎 𝜑𝑖 = (cid:0) − 𝑘 𝑖 (cid:1) (cid:0) Δ 𝑖 + 𝜎 𝑖 (cid:1) , 𝑘 𝑖 ≤ , (18)where again 𝑘 𝑖 can depend on position. We refer to this alternative de-composition as to the generalised 𝑘 -decomposition. The case 𝑘 𝑖 = 𝜐 𝜑𝑖 = 𝑘 𝑖 = 𝜎 𝜑𝑖 =
0; notice that no isotropic rotators can be realised fromequation (18) when Δ 𝑖 <
0, because isotropy ( 𝜎 𝜑𝑖 = 𝜎 𝑖 ) would cor-respond to 𝑘 𝑖 <
0. Moreover, while with the Satoh decomposition aspherical system cannot rotate and is isotropic independently of thevalue of 𝑘 𝑖 , with the generalised 𝑘 -decomposition one can model ro-tating (and anisotropic) spherical systems. An application of this lastcase can be found in exploratory numerical simulations of rotatinggas flows in galaxies of Yoon et al. (2019). Notice that the generaliseddecomposition applied to systems with Δ 𝑖 (cid:29) 𝜎 𝑖 (as for instance thecase of highly flattened discs) reduces to the standard Satoh formula.A more interesting (and delicate) case, requiring particular care in thechoice of the parameter 𝑘 𝑖 , is represented by systems with | Δ 𝑖 | (cid:28) 𝜎 𝑖 ,when we have 𝜐 𝜑𝑖 ∼ 𝑘 𝑖 𝜎 𝑖 , and 𝜎 𝜑𝑖 ∼ ( − 𝑘 𝑖 ) 𝜎 𝑖 . This means that,in order to avoid substantial rotation, for example in almost sphericalsystems (oblate or prolate), 𝑘 𝑖 must be kept small.We finally remark that, for a given multi-component system, it isalso possible to assume a Satoh decomposition for some components,and the generalised decomposition for the others; in analogy withequation (17) it is possible to define a total effective decomposition MNRAS , 1–14 (2019)
C. Caravita, L. Ciotti and S. Pellegrini parameter 𝑘 e as 𝑘 e ≡ 𝜐 𝜑 √ Δ + 𝜎 , 𝜎 𝜑 = (cid:0) − 𝑘 (cid:1) (cid:0) Δ + 𝜎 (cid:1) . (19) We recast here the projection formulae presented in Posacki et al.(2013) for the case of a multi-component system, focusing in partic-ular on how the solutions for the components must be summed toobtain the projected fields of the total stellar distribution. We indi-cate with <, > the scalar product, with n the line-of-sight direction(hereafter los) directed from the observer to the galaxy , and with 𝑙 the integration path along the los. For the ease of notation in thisSection we drop the subscript 𝑖 , so that all the following formulaemust be intended to hold separately for each stellar component 𝜌 ∗ 𝑖 (and of course also for the total 𝜌 ∗ ). We will resume the use of thesubscript 𝑖 at the end of the Section, when we give the expressions forthe projected fields of 𝜌 ∗ as functions of the projected fields of thecomponents. The projection of a stellar density, and of the orderedvelocity 𝝊 = 𝜐 𝜑 e 𝜑 , are Σ ∗ = ∫ ∞−∞ 𝜌 ∗ d 𝑙, Σ ∗ 𝜐 los = ∫ ∞−∞ 𝜌 ∗ 𝜐 𝜑 < e 𝜑 , n > d 𝑙, (20)where e 𝜑 = (− sin 𝜑, cos 𝜑, ) is the unitary vector in the tangentialdirection. From the adopted orientation of n , a positive/negative 𝜐 los indicates a motion receding from/approaching to the observer,respectively. The los velocity dispersion can be written as 𝜎 = 𝜎 + 𝑉 − 𝜐 = 𝑉 − 𝜐 (21)(e.g. Ciotti & Pellegrini 1996; Posacki et al. 2013, CMPZ21), wherefollowing Cappellari (2008) we also defined 𝑉 ≡ 𝜎 + 𝑉 , and Σ ∗ 𝜎 = ∫ ∞−∞ 𝜌 ∗ < 𝜎 n , n > d 𝑙, (22) Σ ∗ 𝑉 = ∫ ∞−∞ 𝜌 ∗ 𝜐 𝜑 < e 𝜑 , n > d 𝑙, (23)where in equation (22) 𝜎 is the 3 × 𝑉 rms and 𝜐 los in general depend on the specific direction n ,and 𝜐 los , 𝜎 P and 𝑉 P on the specific velocity decomposition adopted,but 𝑉 rms is independent of the velocity decomposition. The previousidentities are fully general and hold for a generic inclination of the loswith respect to the galaxy. As JASMINE2 deals with axisymmetricmodels, it is assumed without loss of generality that the los is parallelto the 𝑥 − 𝑧 plane, and the projection plane rotates around the 𝑦 axis.In particular, in the face-on projection, the los is parallel to the 𝑧 axis with n = − e 𝑧 , the projection plane is the 𝑥 − 𝑦 plane, and Σ ∗ = ∫ ∞ 𝜌 ∗ d 𝑧, Σ ∗ 𝜎 = ∫ ∞ 𝜌 ∗ 𝜎 d 𝑧, (24)because 𝜐 los = 𝑉 P =
0, and so 𝜎 los = 𝜎 P . In the edge-on projection(hereafter EO), the los is aligned with the 𝑥 axis with n = − e 𝑥 , theprojection plane coincides with the 𝑦 − 𝑧 plane, and ( cos 𝜑, sin 𝜑 ) = ( 𝑥 / 𝑅, 𝑦 / 𝑅 ) where 𝑅 = √︁ 𝑥 + 𝑦 . Then, from equation (20), Σ ∗ = ∫ ∞ 𝑦 𝜌 ∗ 𝑅 √︁ 𝑅 − 𝑦 d 𝑅, Σ ∗ 𝜐 los = 𝑦 ∫ ∞ 𝑦 𝜌 ∗ 𝜐 𝜑 √︁ 𝑅 − 𝑦 d 𝑅. (25) At variance with the convention adopted in CMPZ21, where n points fromthe galaxy to the observer. Moreover, with some algebra, from equations (22) and (23) we have: Σ ∗ 𝜎 = ∫ ∞ 𝑦 ( 𝑅 − 𝑦 ) 𝜎 + 𝑦 𝜎 𝜑 𝑅 √︁ 𝑅 − 𝑦 𝜌 ∗ d 𝑅, (26) Σ ∗ 𝑉 = 𝑦 ∫ ∞ 𝑦 𝜌 ∗ 𝜐 𝜑 𝑅 √︁ 𝑅 − 𝑦 d 𝑅, (27)so that, from equation (10), equation (21) can be recast in compactform as Σ ∗ 𝜎 = ∫ ∞ 𝑦 𝑅 𝜎 + 𝑦 Δ 𝑅 √︁ 𝑅 − 𝑦 𝜌 ∗ d 𝑅 − Σ ∗ 𝜐 , (28)where the independence of 𝑉 rms from the specific azimuthal velocitydecomposition is apparent.The projection formulae for a multi-component stellar system cannow be easily obtained, for a generic los, just by considering howthe intrinsic quantities add. From equations (1) and (13), and fromequations (20)–(23), it is immediate to see that Σ ∗ = ∑︁ 𝑖 Σ ∗ 𝑖 , 𝜐 los = (cid:205) 𝑖 Σ ∗ 𝑖 𝜐 los 𝑖 Σ ∗ , 𝑉 = (cid:205) 𝑖 Σ ∗ 𝑖 𝑉 𝑖 Σ ∗ , (29)and 𝜎 is given again by equation (21).No difficulty is encountered in the construction of the luminosity-weighted fields analogous to equations (29), by using the surfacebrightness distributions 𝐼 ∗ 𝑖 = Σ ∗ 𝑖 / Υ ∗ 𝑖 , and 𝐼 ∗ = (cid:205) 𝑖 𝐼 ∗ 𝑖 , so the pro-jected stellar mass-to-light ratio Υ ∗ los ≡ Σ ∗ / 𝐼 ∗ , defined in analogywith the local Υ ∗ in equation (4).Summarising, we now have the framework needed to determinethe solution of the Jeans equations once the solutions for the singlecomponents in the total potential are known. In the following Sectionwe describe how the solution of each stellar component is obtained,thanks to the scaling procedure adopted by JASMINE2. One of the most important updates in JASMINE2 is the full useof the scalings allowed by the Poisson and the Jeans equations, sothat with a single run of the Potential & Jeans Solver one can build afamily of solutions that can be combined in PP with different weights,producing several galaxy models with almost no effort.The basic idea is elementary. We recognise that equations (9), at fixed total potential Φ , are invariant for a mass scaling of the density 𝜌 ∗ 𝑖 , i.e. at fixed Φ the derived velocity fields would be independentof the value of 𝑀 ∗ 𝑖 . However, as Φ contains also 𝜙 ∗ 𝑖 , equations (9)obviously are not invariant to such scaling; nonetheless, the 𝑁 + 𝜌 ∗ 𝑖 in the potentials 𝜙 𝑗 𝜕𝜌 ∗ 𝑖 𝜎 𝑖 𝑗 𝜕𝑧 = − 𝜌 ∗ 𝑖 𝜕𝜙 𝑗 𝜕𝑧 ,𝜕𝜌 ∗ 𝑖 𝜎 𝑖 𝑗 𝜕𝑅 = 𝜌 ∗ 𝑖 Δ 𝑖 𝑗 𝑅 − 𝜌 ∗ 𝑖 𝜕𝜙 𝑗 𝜕𝑅 , (30)and their solutions 𝜌 ∗ 𝑖 𝜎 𝑖 𝑗 = ∫ ∞ 𝑧 𝜌 ∗ 𝑖 𝜕𝜙 𝑗 𝜕𝑧 (cid:48) 𝑑𝑧 (cid:48) , 𝜌 ∗ 𝑖 Δ 𝑖 𝑗 = 𝑅 (cid:32) 𝜕𝜌 ∗ 𝑖 𝜎 𝑖 𝑗 𝜕𝑅 + 𝜌 ∗ 𝑖 𝜕𝜙 𝑗 𝜕𝑅 (cid:33) , (31) MNRAS000
0, and so 𝜎 los = 𝜎 P . In the edge-on projection(hereafter EO), the los is aligned with the 𝑥 axis with n = − e 𝑥 , theprojection plane coincides with the 𝑦 − 𝑧 plane, and ( cos 𝜑, sin 𝜑 ) = ( 𝑥 / 𝑅, 𝑦 / 𝑅 ) where 𝑅 = √︁ 𝑥 + 𝑦 . Then, from equation (20), Σ ∗ = ∫ ∞ 𝑦 𝜌 ∗ 𝑅 √︁ 𝑅 − 𝑦 d 𝑅, Σ ∗ 𝜐 los = 𝑦 ∫ ∞ 𝑦 𝜌 ∗ 𝜐 𝜑 √︁ 𝑅 − 𝑦 d 𝑅. (25) At variance with the convention adopted in CMPZ21, where n points fromthe galaxy to the observer. Moreover, with some algebra, from equations (22) and (23) we have: Σ ∗ 𝜎 = ∫ ∞ 𝑦 ( 𝑅 − 𝑦 ) 𝜎 + 𝑦 𝜎 𝜑 𝑅 √︁ 𝑅 − 𝑦 𝜌 ∗ d 𝑅, (26) Σ ∗ 𝑉 = 𝑦 ∫ ∞ 𝑦 𝜌 ∗ 𝜐 𝜑 𝑅 √︁ 𝑅 − 𝑦 d 𝑅, (27)so that, from equation (10), equation (21) can be recast in compactform as Σ ∗ 𝜎 = ∫ ∞ 𝑦 𝑅 𝜎 + 𝑦 Δ 𝑅 √︁ 𝑅 − 𝑦 𝜌 ∗ d 𝑅 − Σ ∗ 𝜐 , (28)where the independence of 𝑉 rms from the specific azimuthal velocitydecomposition is apparent.The projection formulae for a multi-component stellar system cannow be easily obtained, for a generic los, just by considering howthe intrinsic quantities add. From equations (1) and (13), and fromequations (20)–(23), it is immediate to see that Σ ∗ = ∑︁ 𝑖 Σ ∗ 𝑖 , 𝜐 los = (cid:205) 𝑖 Σ ∗ 𝑖 𝜐 los 𝑖 Σ ∗ , 𝑉 = (cid:205) 𝑖 Σ ∗ 𝑖 𝑉 𝑖 Σ ∗ , (29)and 𝜎 is given again by equation (21).No difficulty is encountered in the construction of the luminosity-weighted fields analogous to equations (29), by using the surfacebrightness distributions 𝐼 ∗ 𝑖 = Σ ∗ 𝑖 / Υ ∗ 𝑖 , and 𝐼 ∗ = (cid:205) 𝑖 𝐼 ∗ 𝑖 , so the pro-jected stellar mass-to-light ratio Υ ∗ los ≡ Σ ∗ / 𝐼 ∗ , defined in analogywith the local Υ ∗ in equation (4).Summarising, we now have the framework needed to determinethe solution of the Jeans equations once the solutions for the singlecomponents in the total potential are known. In the following Sectionwe describe how the solution of each stellar component is obtained,thanks to the scaling procedure adopted by JASMINE2. One of the most important updates in JASMINE2 is the full useof the scalings allowed by the Poisson and the Jeans equations, sothat with a single run of the Potential & Jeans Solver one can build afamily of solutions that can be combined in PP with different weights,producing several galaxy models with almost no effort.The basic idea is elementary. We recognise that equations (9), at fixed total potential Φ , are invariant for a mass scaling of the density 𝜌 ∗ 𝑖 , i.e. at fixed Φ the derived velocity fields would be independentof the value of 𝑀 ∗ 𝑖 . However, as Φ contains also 𝜙 ∗ 𝑖 , equations (9)obviously are not invariant to such scaling; nonetheless, the 𝑁 + 𝜌 ∗ 𝑖 in the potentials 𝜙 𝑗 𝜕𝜌 ∗ 𝑖 𝜎 𝑖 𝑗 𝜕𝑧 = − 𝜌 ∗ 𝑖 𝜕𝜙 𝑗 𝜕𝑧 ,𝜕𝜌 ∗ 𝑖 𝜎 𝑖 𝑗 𝜕𝑅 = 𝜌 ∗ 𝑖 Δ 𝑖 𝑗 𝑅 − 𝜌 ∗ 𝑖 𝜕𝜙 𝑗 𝜕𝑅 , (30)and their solutions 𝜌 ∗ 𝑖 𝜎 𝑖 𝑗 = ∫ ∞ 𝑧 𝜌 ∗ 𝑖 𝜕𝜙 𝑗 𝜕𝑧 (cid:48) 𝑑𝑧 (cid:48) , 𝜌 ∗ 𝑖 Δ 𝑖 𝑗 = 𝑅 (cid:32) 𝜕𝜌 ∗ 𝑖 𝜎 𝑖 𝑗 𝜕𝑅 + 𝜌 ∗ 𝑖 𝜕𝜙 𝑗 𝜕𝑅 (cid:33) , (31) MNRAS000 , 1–14 (2019) eans modeling with multiple stellar populations with 𝜌 ∗ 𝑖 𝜎 𝑖 𝑗 → 𝑧 → ∞ , do have important scaling propertiesthat will be exploited in Section 3.1. We note that here and in thefollowing the double subscript in 𝜎 𝑖 𝑗 does not refer to the tensorialnature of the velocity dispersion, but just identifies the solution ofthe 𝑖 -th stellar component in the 𝑗 -th potential component.Leaving aside for the moment the scaling properties of equations(30) and (31), it is obvious that the sums 𝜎 𝑖 = ∑︁ 𝑗 𝜎 𝑖 𝑗 , Δ 𝑖 = ∑︁ 𝑗 Δ 𝑖 𝑗 , (32)are the solution of equations (9), as can be demonstrated, first bysumming over 𝑗 the 𝑁 + 𝜎 𝑖 and Δ 𝑖 performed in equations(30) over the 𝑁 + 𝜙 𝑗 , with the decompositionof 𝜎 and Δ performed in equations (9) over the 𝑁 stellar compo-nents 𝜌 ∗ 𝑖 , there is a fundamental conceptual difference between thetwo decompositions. In fact, equations (9) are true moments of theCollisionless Boltzmann Equation obeyed by the distribution func-tions 𝑓 𝑖 in the total potential, and so they have a sort of autonomousphysical meaning; equations (30), instead, are just a mathematicaldecomposition over the different 𝜙 𝑗 of the Jeans equations for 𝜌 ∗ 𝑖 .As a consequence, phase-space consistency arguments apply to thesolution of equations (9), but not to 𝜎 𝑖 𝑗 and Δ 𝑖 𝑗 separately: as far asthe fields 𝜎 𝑖 and Δ 𝑖 are physically acceptable, the model is also ac-ceptable, independently of the specific properties of its components 𝜎 𝑖 𝑗 and Δ 𝑖 𝑗 . We distinguish three groups of model parameters, which are theinputs of JASMINE2 for the construction a multi-component model,and are summarised in Table 1.In the first group we have the physical scales 𝑀 ∗ and 𝑟 ∗ , i.e. thetotal stellar mass and its scale-length. All the density and potentialcomponents are made dimensionless by scaling them to the quantities 𝜌 n ≡ 𝑀 ∗ 𝜋𝑟 ∗ , 𝜙 n ≡ 𝐺 𝑀 ∗ 𝑟 ∗ , (33)and the numerical grids used by JASMINE2 are normalised to 𝑟 ∗ ,with ˜ 𝑅 ≡ 𝑅 / 𝑟 ∗ , and ˜ 𝑧 ≡ 𝑧 / 𝑟 ∗ . The dimensionless grids (usuallybilogarithmic and with a few hundreds of points in ˜ 𝑅 and ˜ 𝑧 , rangingfrom ≈ − or less at the origin, up to ≈ or more at theouter edge) guarantee the same resolution independently of the actualphysical size of the model, measured by 𝑟 ∗ . Even though the physicalscales are logically introduced first, the values of 𝑀 ∗ and 𝑟 ∗ (andso of 𝜌 n and 𝜙 n ) are fixed in the last step of the model construction,at the end of the PP (see Table 1) . In this way, different physicalrealisations (in size and total mass) can be obtained for the samemulti-component galaxy model.In the second group of parameters we have the relative massweights R 𝑖 ≡ 𝑀 ∗ 𝑖 / 𝑀 ∗ , R h ≡ 𝑀 h / 𝑀 ∗ , R BH ≡ 𝑀 BH / 𝑀 ∗ of the The situation is somewhat similar to that faced when decomposing a posi-tive density distribution over some prescribed set of functions (e.g, sphericalharmonics), when the basis functions can present regions of negative densi-ties. different components, the mass-to-light ratios Υ ∗ 𝑖 , and the param-eters 𝑘 𝑖 appearing in equations (16) and (18) for the kinematicaldecomposition of the azimuthal motions. By definition ∑︁ 𝑖 R 𝑖 = , (34)and in full generality we write 𝜌 ∗ 𝑖 = 𝜌 n R 𝑖 ˜ 𝜌 ∗ 𝑖 , 𝜌 h = 𝜌 n R h ˜ 𝜌 h , ˜ 𝜌 ∗ = ∑︁ 𝑖 R 𝑖 ˜ 𝜌 ∗ 𝑖 . (35)where ˜ 𝜌 ∗ 𝑖 and ˜ 𝜌 h are the scaled density distributions, and ˜ 𝜌 ∗ = 𝜌 ∗ / 𝜌 n is the dimensionless total stellar density. Notice that from equation(33) the volume integrals of ˜ 𝜌 ∗ 𝑖 over the whole dimensionless nu-merical grid evaluate to 4 𝜋 by construction. Similarly, 𝜙 ∗ 𝑖 = 𝜙 n R 𝑖 ˜ 𝜙 ∗ 𝑖 , 𝜙 h = 𝜙 n R h ˜ 𝜙 h , 𝜙 BH = 𝜙 n R BH ˜ 𝜙 BH , ˜ Φ = ∑︁ 𝑗 R 𝑗 ˜ 𝜙 𝑗 , (36)where ˜ Φ = Φ / 𝜙 n , and so 𝜎 𝑖 = 𝜙 n ˜ 𝜎 𝑖 , Δ 𝑖 = 𝜙 n ˜ Δ 𝑖 . (37)Finally, from the assumption of a constant mass-to-light ratio Υ ∗ 𝑖 foreach stellar component, 𝜈 ∗ 𝑖 = 𝜌 n R 𝑖 Υ ∗ 𝑖 ˜ 𝜌 ∗ 𝑖 , 𝐿 𝑖 = R 𝑖 Υ ∗ 𝑖 𝑀 ∗ . (38) The values of the weights are chosen in PP (see Table 1), because achange in their values, and in the kinematical decompositions, doesnot require to recompute the potentials and solve again the Jeansequations.
This possibility allows for a fast construction of differentmodels belonging to the same family.One family indeed is characterised by the choice of the third groupof parameters, to be performed at the beginning of the model con-struction: the structural parameters of the scaled density compo-nents ˜ 𝜌 ∗ 𝑖 and ˜ 𝜌 h , that in full generality we indicate with the symbols 𝜉 𝑖 ≡ 𝑟 ∗ 𝑖 / 𝑟 ∗ and 𝜉 h ≡ 𝑟 h / 𝑟 ∗ for the different scale-lengths, and with 𝑞 𝑖 and 𝑞 h for other parameters that determine the shape of the scaleddensities (for example the flattenings in case of ellipsoidal densitydistributions). The values of the structural parameters must be as-signed in order to run the Potential & Jeans Solver (see Table 1),and in general a change in some of their values requires a newcomputation of the potentials and of the Jeans solutions.3.1.1 The Potential & Jeans Solver
For a proposed set of values for the structural parameters, the scaledJeans equations solved by JASMINE2 are obtained from equations(30) and equations (35) and (36). In practice, for 𝑁 assigned scaledstellar components ˜ 𝜌 ∗ 𝑖 and a scaled dark matter halo ˜ 𝜌 h , the Potential& Jeans Solver first computes the scaled potentials ˜ 𝜙 ∗ 𝑖 and ˜ 𝜙 h , andthen solves the 𝑁 × ( 𝑁 + ) pairs of scaled equations (30), one for each˜ 𝜌 ∗ 𝑖 in the potential ˜ 𝜙 𝑗 (including ˜ 𝜙 BH ), over the dimensionless grid ( ˜ 𝑅, ˜ 𝑧 ) ; thus the scaled fields ˜ 𝜎 𝑖 𝑗 and ˜ Δ 𝑖 𝑗 are obtained. The possibilityto solve equations (30) without choosing R 𝑖 and R 𝑗 is due to the factthat, on one hand, the weights R 𝑖 appear linearly in both sides ofequations (30); on the other hand, ˜ 𝜎 𝑖 𝑗 and ˜ Δ 𝑖 𝑗 scale linearly with R 𝑗 , once the boundary condition is fixed to zero at infinity.The details of the numerical procedure for the computation of thepotentials (based on the use of staggered grids and of complete ellipticintegrals) and for the solution of the Jeans equations are described inPosacki et al. (2013). As already remarked, the numerical evaluationof the potential is the most time consuming part of the construction of MNRAS , 1–14 (2019)
C. Caravita, L. Ciotti and S. Pellegrini
JASMINE2 INPUTSPotential & Jeans SolverStructural parametersScaled stellar and DM densities ˜ 𝜌 ∗ 𝑖 , ˜ 𝜌 h Scale-length ratios 𝜉 𝑖 = 𝑟 ∗ 𝑖 𝑟 ∗ , 𝜉 h = 𝑟 h 𝑟 ∗ , . . . Shape parameters 𝑞 𝑖 , 𝑞 h , . . . Post-ProcessingWeightsMass ratios R 𝑖 = 𝑀 ∗ 𝑖 𝑀 ∗ , R h = 𝑀 h 𝑀 ∗ , R BH = 𝑀 BH 𝑀 ∗ Mass-to-light ratios Υ ∗ 𝑖 = 𝑀 ∗ 𝑖 𝐿 𝑖 Kinematical decompositions 𝑘 𝑖 , 𝜆 𝑖 , 𝛿 𝑖 Post-ProcessingPhysical ScalesTotal stellar mass 𝑀 ∗ Total stellar density scale-length 𝑟 ∗ Table 1.
The three main steps involved in the construction of a multi-component model by JASMINE2, listed from top to bottom in the orderin which they are considered by the code. a model. For this reason JASMINE2 contains a continuously updatedlibrary of analytical density-potential pairs available in the literature(and in some cases also based on homoeoidal expansion, Ciotti &Bertin 2005), so that one can choose between the numerical compu-tation of the potential and (when available) the use of the analyticalpotential.
As described in the previous Section, for a given multi-componentmodel of assigned ˜ 𝜌 ∗ 𝑖 , ˜ 𝜌 h , and with a central BH, the Potential &Jeans Solver gives the solution ˜ 𝜎 𝑖 𝑗 and ˜ Δ 𝑖 𝑗 of the scaled form ofequations (30). These solutions are then combined in PP, with theassignment of the mass ratios R 𝑗 , so that the solution 𝜎 𝑖 and Δ 𝑖 of equations (9) for 𝜌 ∗ 𝑖 is obtained, according to equations (32) and(37), as˜ 𝜎 𝑖 = ∑︁ 𝑗 R 𝑗 ˜ 𝜎 𝑖 𝑗 , ˜ Δ 𝑖 = ∑︁ 𝑗 R 𝑗 ˜ Δ 𝑖 𝑗 . (39)At this stage, as discussed in Section 2.3, the PP performs a positivitycheck of ˜ Δ 𝑖 and ˜ Δ 𝑖 + ˜ 𝜎 𝑖 : in case of negativity of the last quantity, anew choice of the weights R 𝑗 is made, until positivity is reached. Ifpositivity cannot be obtained for acceptable choices of R 𝑗 , then themulti-component model is discarded as unphysical.Once the mass weights are assigned and the positivity check ispassed, the PP requires the parameter 𝑘 𝑖 for the kinematical decom-position, that gives the scaled azimuthal velocity fields 𝜐 𝜑𝑖 and 𝜎 𝜑𝑖 .In full generality, we define each decomposition parameter as 𝑘 𝑖 ( 𝑅, 𝑧 ) = 𝜆 𝑖 𝛿 𝑖 ( 𝑅, 𝑧 ) , (40)where 𝜆 𝑖 is a constant weight, and 𝛿 𝑖 ( 𝑅, 𝑧 ) is a position-dependentfunction; the standard Satoh parameter is obtained with 𝛿 𝑖 = 𝜆 𝑖 = 𝑘 𝑖 . The benefit of this factorisation is due to the fact that theprojection formula of 𝜐 𝜑𝑖 in equation (20) for a given 𝛿 𝑖 ( 𝑅, 𝑧 ) scaleswith 𝜆 𝑖 , so that we can set the value of 𝜆 𝑖 after having computedthe projection integral. As projections represent the second mosttime-consuming step of JASMINE2, the possibility to choose (andchange) 𝜆 𝑖 after projections is a significant advantage. Note that, at variance with what happens for the fields ˜ 𝜎 𝑖 , ˜ Δ 𝑖 , ˜ 𝜐 𝜑𝑖 and ˜ 𝜎 𝜑𝑖 , themass weights R 𝑗 enter the expression of ˜ 𝜐 𝜑𝑖 under a square root(see equations 16 and 18). This implies that the R 𝑗 must be chosenbefore calculating the projections that use 𝜐 𝜑𝑖 . In other words, thepossibility to modify the values of R 𝑗 in PP, allowed by the ” 𝑖 𝑗 -decomposition”, ends with the computation of the scaled fields inequations (39).Once we have obtained the intrinsic and projected fields of each 𝜌 ∗ 𝑖 , the last steps are to combine them to calculate the total (mass-and luminosity-weighted) intrinsic and projected fields of 𝜌 ∗ (respec-tively from equations 13, 14, and equations 20, 21, 29), and finallyto choose the physical scales 𝑀 ∗ and 𝑟 ∗ . Summarising, a family of multi-component galaxy models is definedby the choice of 𝑁 scaled stellar density components ˜ 𝜌 ∗ 𝑖 , a scaled DMhalo ˜ 𝜌 h , and a central BH. The Potential & Jeans Solver computes theassociated scaled potentials ˜ 𝜙 𝑗 , and then solves the 𝑁 × ( 𝑁 + ) Jeansequations (30) in their scaled form. In the subsequent PP, specificvalues of the mass ratios R 𝑖 , R h , R BH , of the mass-to-light ratios Υ ∗ 𝑖 ,and of the kinematical decompositions with the parameters 𝑘 𝑖 , arefixed, thus defining a specific model in the same family. The solutionof the Jeans equations for the total density distribution is recoveredas (mass- or luminosity-) weighted sums of the scaled solutions, andthe projections along a given line-of-sight are performed. The valuesof the total stellar mass 𝑀 ∗ , and of its scale-length 𝑟 ∗ , complete theconstruction of the model.There are at least two significant advantages in this procedure,when compared with a straightforward integration of the Jeans equa-tions for a multi-component galaxy model. First, the gravitationalpotentials of each stellar component and of the DM halo need not tobe recalculated every time the weights are changed in PP; thus therun of the most time expensive part of JASMINE2 is required justonce for all the models in the same family. Second, the possibility tochoose the weight parameters in PP allows for a fast exploration ofthe parameter space (that, for multi-component models, can be verylarge). Qualitatively, the 𝑁 × ( 𝑁 + ) set of the 𝑖 𝑗 -th scaled solutionsof the Jeans equations for each 𝑖 -th density component in each 𝑗 -th potential component, can be interpreted as basis vectors that aresuccessively linearly combined with different weights, to obtain aspecific solution belonging to a family of multi-component models. In order to illustrate the new features and potentialities of JASMINE2,we describe in some detail the building of three multi-componentgalaxy models. All three models are made of two stellar distributions,to which a DM halo with a spherical Navarro-Frenk-White profile(Navarro et al. 1996, hereafter NFW) and a central supermassiveBH are added. In the first model (hereafter JJE) the total sphericalstellar profile and an ellipsoidal stellar component, both with a Jaffe(1983) profile, are assigned; if the dark mass is set to zero, this modelreduces to the JJe models of CMPZ21. The second model (hereafter The effective radius 𝑅 e of the total stellar distribution is obviously anotherimportant quantity that cannot be obtained as a linear combination of theeffective radii of the stellar components, and it can only be computed afterthe choice of the weights R 𝑖 and Υ ∗ 𝑖 .MNRAS , 1–14 (2019) eans modeling with multiple stellar populations JHD) consists of an ellipsoidal Jaffe stellar density distribution, thatrepresents a light stellar halo, coupled with a heavy Miyamoto-Nagaistellar disc (Miyamoto & Nagai 1975, hereafter MN). In the thirdmodel (hereafter JLD), the ellipsoidal Jaffe component dominates,while a small MN inner disc is counter-rotating. These three models,discussed in the next Sections, are intended to represent featuresobserved in real galaxies, but they are not designed to reproducespecific objects.
Before the description of a representative JJE model (Section 4.2.1),we illustrate some general properties of JJE models, in particular incomparison with JJe models, and three experiments to test the code(Section 4.1.1). Indeed, JJE models are a natural generalisation ofJJe models presented in CMPZ21: as these latter describe quite wellreal elliptical galaxies, and several of their dynamical properties canbe expressed in analytical form, they also represent an obvious testfor the proper working of JASMINE2.To better appreciate the properties of JJE models, we recall themain properties (and limitations) of JJe models. These are con-structed by assigning a total density 𝜌 ∗ following the axisymmetricellipsoidal generalisation of the Jaffe model, and another axisymmet-ric ellipsoidal Jaffe distribution 𝜌 ∗ , with different flattening, scale-length and total mass; in CMPZ21 the density distribution 𝜌 ∗ − 𝜌 ∗ ( = 𝜌 ∗ in the current notation), is interpreted as a DM halo; finally,a central BH is added to the system. The analytical conditions on 𝜌 ∗ to guarantee the positivity of 𝜌 ∗ are given, and then the Jeansequations for 𝜌 ∗ are solved in analytical closed form, by using ho-moeoidal expansion, truncated at the linear order in the flattenings of 𝜌 ∗ and 𝜌 ∗ . Albeit several properties of JJe models can be expressedin analytical form (making these models quite useful in numericalsimulations of gas flows in galaxies, see e.g. Gan et al. 2019a,b), afew important shortcomings still affect them: i) the Jeans equationsfor 𝜌 ∗ are integrated in the homoeoidal expansion limit, retainingonly linear terms in the flattenings, and they have not been studiedfor the difference component 𝜌 ∗ ; ii) projected kinematical fields of 𝜌 ∗ can be obtained in analytical form only as asymptotic formulaeat the center and at large radii. JASMINE2 is then the obvious toolto address the two points above.Here we generalise the JJe models to JJE models, by consideringfor the total 𝜌 ∗ an ellipsoidal Jaffe profile, of total mass 𝑀 ∗ , scale-length 𝜉 , and flattening 𝑞 : 𝜌 ∗ ( 𝑅, 𝑧 ) = 𝜌 n 𝜉𝑞 𝑚 ( 𝜉 + 𝑚 ) , 𝑚 = ˜ 𝑅 + ˜ 𝑧 𝑞 . (41)Therefore, at variance with JJe models, in JJE models the total Jaffemass distribution is purely stellar. We then consider another ellip-soidal Jaffe density profile, of total mass 𝑀 ∗ = R 𝑀 ∗ , scale-length 𝑟 ∗ = 𝜉 𝑟 ∗ , and flattening 𝑞 : 𝜌 ∗ ( 𝑅, 𝑧 ) = 𝜌 n R 𝜉 𝑞 𝑚 ( 𝜉 + 𝑚 ) , 𝑚 = ˜ 𝑅 + ˜ 𝑧 𝑞 . (42)The second stellar component is then defined as 𝜌 ∗ ( 𝑅, 𝑧 ) = 𝜌 ∗ ( 𝑅, 𝑧 ) − 𝜌 ∗ ( 𝑅, 𝑧 ) , (43)with 𝑀 ∗ = 𝑀 ∗ − 𝑀 ∗ = ( − R ) 𝑀 ∗ = R 𝑀 ∗ , in agreement withequation (34). Notice that 𝜌 ∗ is not an ellipsoid, unless 𝑞 = 𝑞 , and In the spherical limit, and in the assumption of constant mass-to-light ratio,the Jaffe scale radius 𝑟 ∗ is related to the effective radius by 𝑅 e (cid:39) . 𝑟 ∗ . Model 𝜌 ∗ 𝜌 ∗ JJE Jaffe 𝜌 ∗ ( 𝜉 = , 𝑞 = ) − 𝜌 ∗ 𝜉 = . 𝑞 = . R = . R = . Υ ∗ = Υ ∗ = 𝑘 = . 𝑘 = . 𝜉 = 𝑏 = . 𝑞 = . 𝑞 = R = . R = . Υ ∗ = Υ ∗ = 𝑘 = . 𝑘 = . 𝜉 = 𝑏 = . 𝑞 = . 𝑞 = R = . R = . Υ ∗ = Υ ∗ = 𝑘 = . 𝑘 ( 𝑅, 𝑧 ) Table 2.
The parameters for the stellar components of the illustrative JJE, JHDand JLD models (Sections 4.1 and 4.2). In the JJE model, the component 𝜌 ∗ is obtained as difference between a total spherical ( 𝑞 =
1) Jaffe profile 𝜌 ∗ ,with scale-length 𝜉 =
1, and a small and light Jaffe ellipsoidal component 𝜌 ∗ (as done for JJe models in CMPZ21). The standard Satoh 𝑘 -decompositionin equation (16) for 𝜌 ∗ , and the generalised 𝑘 -decomposition in equation(18) for 𝜌 ∗ , are adopted. In the JHD model, an ellipsoidal Jaffe distributionis coupled to a massive and quite flat ( 𝑞 = 𝑎 / 𝑏 =
10) MN disc; in bothcomponents a generalised 𝑘 -decomposition is adopted. In the JLD model, theellipsoidal Jaffe component has the same flattening and size as in the JHDmodel, but the disc is significantly smaller, and counter-rotates in the innerregions, with the position-dependent Satoh parameter in equation (47), whilea constant Satoh parameter is applied to the Jaffe component. In all models,the DM halo has a spherical NFW profile with 𝜉 h = . 𝑐 = R h = R BH = . even in this case 𝜌 ∗ is not a Jaffe ellipsoid, unless 𝜉 = 𝜉 . As exten-sively discussed in CMPZ21, 𝜌 ∗ could be negative somewhere (andso unphysical) for some choices of R , 𝜉 and 𝑞 . Remarkably, theconditions required to assure 𝜌 ∗ ≥ 𝜌 ∗ is embedded in a NFW DM halo (spher-ically symmetric for simplicity), of mass 𝑀 h ( 𝑟 t ) = R h 𝑀 ∗ enclosedwithin a truncation radius 𝑟 t , scale-length 𝑟 h = 𝜉 h 𝑟 ∗ , and concentra-tion 𝑐 ≡ 𝑟 t / 𝑟 h : 𝜌 h ( 𝑟 ) = 𝜌 n R h 𝑓 ( 𝑐 ) 𝑠 ( 𝜉 h + 𝑠 ) , 𝜙 h ( 𝑟 ) = − 𝜙 n R h ln ( + 𝑠 / 𝜉 h ) 𝑓 ( 𝑐 ) 𝑠 (44)where 𝑠 ≡ 𝑟 / 𝑟 ∗ , and 𝑓 ( 𝑐 ) = ln ( + 𝑐 ) − 𝑐 /( + 𝑐 ) . We complete themodel with a central BH of mass 𝑀 BH = R BH 𝑀 ∗ .Summarising, JJE models are determined, besides the total stellarmass and scale-length, 𝑀 ∗ and 𝑟 ∗ , by the two parameters 𝜉 and 𝑞 for 𝜌 ∗ , the five parameters 𝑞 , 𝜉 , R , Υ ∗ , 𝑘 for 𝜌 ∗ , the twoparameters Υ ∗ , 𝑘 for 𝜌 ∗ , the three DM parameters 𝜉 h , 𝑐 , R h , andthe BH mass weight R BH (see Table 2 for a specific JJE model).JASMINE2 further generalises JJE models, with the addition of asecond DM component given by a shallow and very extended quasi-isothermal halo, as useful in simulations of gas flows in galaxiesresiding in groups or clusters (Ciotti et al. 2021 in preparation). MNRAS , 1–14 (2019)
C. Caravita, L. Ciotti and S. Pellegrini
According to the previous discussion, we can use JASMINE2 in twodifferent tests: we can give in input the homoeoidal expansion of thedensity-potential pairs, and compare the numerical solution with thatof CMPZ21, to check the importance of quadratic flattening terms inthe solution of the Jeans equations for JJe models; and we can give ininput the true ellipsoidal models, to check how well the homoeoidalexpansion reproduces, as a function of the flattenings, the internaldynamics of ellipsoidal systems. At the same time, of course, theanalytical results of JJe models are used to test the accuracy ofJASMINE2. To compare the numerical results of JASMINE2 withthe analytical formulae in CMPZ21, we construct a JJe system with R h =
0, and we focus on the component 𝜌 ∗ , for which the analyticalsolution of the Jeans equations is available.In the first experiment, we feed JASMINE2 with the homoeoidalexpansion for 𝜌 ∗ and 𝜌 ∗ , and we compare the numerical resultsof integration of the Jeans equations with the analytical results: weobtain excellent agreement for all the kinematical fields, better thana fraction of percent over the whole numerical grid, ranging from ≈ − to ≈
70 (in units of 𝑟 ∗ ). The results tend to be slightly morediscrepant at increasing flattenings, as expected since the analyticalresults in CMPZ21 are limited to the linear order in the flattenings,while JASMINE2 takes automatically into account also the second order terms when integrating the Jeans equations (due to the productbetween the density and the gradient of the potential). By increas-ing the numerical resolution in the central regions, and moving tosmaller and smaller distances from the center, we also verified thatthe asymptotic formulae for the projected kinematical fields werealso perfectly recovered numerically. This first experiments addedconfidence that JASMINE2 is working properly and with high accu-racy, but also provided a further support that the (quite cumbersome)analytical formulae in CMPZ21 are actually correct .In a second experiment, we compared the numerical results of JAS-MINE2 for the true ellipsoidal JJe models with the analytical resultsin CMPZ21 obtained from homoeoidal expansion, therefore movingbeyond the effect of second order approximation in the flatteningsexplored in the first experiment. We found that for relatively smallflattenings (corresponding to a ≈ E2-E3 galaxies) the homoeoidalexpansion truncated at the linear order provides quite good results,even when compared with true ellipsoidal models, with the mostsignificant discrepancy in the intermediate regions.In a third experiment, we finally compared the JASMINE2 solu-tions for the true ellipsoidal JJe models with the solutions obtainedby evaluating the gravitational potential with the standard ellipsoidalformula (e.g. BT08), and then integrating the Jeans equations numer-ically with
Mathematica , reaching perfect agreement.
We move now to illustrate the main properties of a specific JJE model(see Table 2). The total stellar distribution 𝜌 ∗ has a spherical Jaffeprofile obtained from equation (41) with 𝜉 = 𝑞 =
1; this quiteartificial case allows us to discuss some subtleties that can occurto the kinematical decomposition in multi-component systems. Thestellar component 𝜌 ∗ is obtained from equation (42) with 𝜉 = . 𝑞 = . R is 0 .
04, and Υ ∗ =
2, i.e. it is a quite small ellipsoidal These experiments are similar in the approach to those already performedwith JASMINE by using the analytical results of Smet et al. (2015): we recallthat the numerical integration of the potential in JASMINE2 uses the sameroutines of JASMINE. distribution at the center of the galaxy; note that, from equation(A15), the maximum possible value of R is 0 .
08. The component 𝜌 ∗ accounts for the remaining 96% of the total stellar mass of thegalaxy, and Υ ∗ =
6, so that 𝜌 ∗ could represent an elliptical galaxywith a central and younger stellar system. We add the spherical NFWDM halo, given by equation (44), with 𝜉 h = . 𝑐 =
10, and R h = 𝑅 e (see Footnote 3) is ≈ .
45 of the totalmass. Finally, in agreement with BH-galaxy scaling relations (seee.g. Kormendy & Ho 2013), the mass of the central BH is fixed to R BH = . 𝜌 ∗ and ˜ 𝜌 ∗ , and of the total stellar density˜ 𝜌 ∗ = 𝜌 ∗ / 𝜌 n . Being this last spherical, and ˜ 𝜌 ∗ oblate, ˜ 𝜌 ∗ in itscentral regions is slightly prolate, and this will affect its kinematicalfields, as anticipated in Section 2.3. Additional information on themodel structure is provided in the first column of Figure 2: the toppanel shows the radial profiles in the equatorial plane of R ˜ 𝜌 ∗ , R ˜ 𝜌 ∗ , ˜ 𝜌 ∗ , and R h ˜ 𝜌 h . The total 𝜌 ∗ is almost coincident with 𝜌 ∗ ,except for the central regions, where 𝜌 ∗ and 𝜌 ∗ are comparable. TheDM density 𝜌 h overcomes 𝜌 ∗ outside ≈ . 𝑅 e . The bottom panelshows the radial profiles in the equatorial plane of the contributionsto the circular velocity due to the various mass components: the BHcontribution is dominant in the inner regions, the DM in the outerregions, while at intermediate distances from the centre the resultingcircular velocity is quite flat.Similar trends can be seen in the radial profiles of the velocity fieldsin the equatorial plane of Figure 3, where, in the first column fromtop to bottom, we show the rotational velocity, the vertical velocitydispersion, and the azimuthal velocity dispersion, of 𝜌 ∗ and 𝜌 ∗ , andthe total mass-weighted and luminosity-weighted fields. Note that inthe three panels the vertical scale is the same, and the resulting sys-tem is clearly a slow rotator. This JJE model offers the opportunityto apply the generalised 𝑘 -decomposition in equation (18); becausethe field Δ , associated to the slightly prolate 𝜌 ∗ , is negative in thecentral regions. As discussed in Section 2.3, we verified that v 𝜑 isnowhere negative, and then we adopted the generalised decomposi-tion, with a quite small 𝑘 = .
2. The field Δ instead is everywherepositive, as expected, and so we adopted the standard Satoh formulawith 𝑘 = .
5. In the velocity profiles, the effect of the central BH isclearly visible; for example, the velocity dispersion profile of a Jaffemodel with R BH = ≈ 𝑅 e , are almost coincidentwith the profiles of the more massive component 𝜌 ∗ , in both themass-weighted and luminosity-weighted cases; this is not surprising,because in these regions 𝜌 ∗ coincides with 𝜌 ∗ (see Figures 1 and 2).The situation is different in the inner regions, where 𝜌 ∗ and 𝜌 ∗ arecomparable: here the total velocities have intermediate values, withthe luminosity-weighted profiles closer to the profiles of 𝜌 ∗ with thesmaller mass-to-light ratio.As an illustration of the projection procedure, in the first row ofFigure 4, we show the EO projected luminosity-weighted fields 𝜐 los L and 𝜎 los L , with the superimposed dotted contours representing thegalaxy isophotes of the surface brightness 𝐼 ∗ . The slow rotationof the model is apparent from the colorbar values, indeed 𝜐 los L is everywhere lower than 𝜎 los L . A curious feature is the slightlyvertically elongated shape of 𝜎 los L : this is not due to the prolateshape of 𝜌 ∗ in the central regions, but it is an effect of the generalised 𝑘 -decomposition, coupled with the fact that Δ is almost null in theexternal regions, and so here 𝜐 𝜑 ∼ 𝑘 𝜎 , as introduced in Section2.3. For example, an increase in 𝑘 would lead to an increase ofthe rotation in the external regions, with correspondent decrease of MNRAS000
5. In the velocity profiles, the effect of the central BH isclearly visible; for example, the velocity dispersion profile of a Jaffemodel with R BH = ≈ 𝑅 e , are almost coincidentwith the profiles of the more massive component 𝜌 ∗ , in both themass-weighted and luminosity-weighted cases; this is not surprising,because in these regions 𝜌 ∗ coincides with 𝜌 ∗ (see Figures 1 and 2).The situation is different in the inner regions, where 𝜌 ∗ and 𝜌 ∗ arecomparable: here the total velocities have intermediate values, withthe luminosity-weighted profiles closer to the profiles of 𝜌 ∗ with thesmaller mass-to-light ratio.As an illustration of the projection procedure, in the first row ofFigure 4, we show the EO projected luminosity-weighted fields 𝜐 los L and 𝜎 los L , with the superimposed dotted contours representing thegalaxy isophotes of the surface brightness 𝐼 ∗ . The slow rotationof the model is apparent from the colorbar values, indeed 𝜐 los L is everywhere lower than 𝜎 los L . A curious feature is the slightlyvertically elongated shape of 𝜎 los L : this is not due to the prolateshape of 𝜌 ∗ in the central regions, but it is an effect of the generalised 𝑘 -decomposition, coupled with the fact that Δ is almost null in theexternal regions, and so here 𝜐 𝜑 ∼ 𝑘 𝜎 , as introduced in Section2.3. For example, an increase in 𝑘 would lead to an increase ofthe rotation in the external regions, with correspondent decrease of MNRAS000 , 1–14 (2019) eans modeling with multiple stellar populations Figure 1.
The scaled stellar distributions ˜ 𝜌 ∗ , ˜ 𝜌 ∗ , and the dimensionless total stellar distribution ˜ 𝜌 ∗ = 𝜌 ∗ / 𝜌 n , of the three models of Table 2. The dotted contourscorrespond to the isodensities, with values spaced by 1 dex. 𝜎 los L , and with the net result of a more elongation of 𝜎 los L in thecentral regions. JHD and JLD models consist of a stellar profile 𝜌 ∗ given again bythe ellipsoidal Jaffe model in equation (42), coupled with a stellar MN disc 𝜌 ∗ , of total mass 𝑀 ∗ = R 𝑀 ∗ , and scale-lengths 𝑎 = ˜ 𝑎𝑟 ∗ , 𝑏 = ˜ 𝑏𝑟 ∗ : 𝜌 ∗ ( 𝑅, 𝑧 ) = 𝜌 n R ˜ 𝑏 ˜ 𝑎 ˜ 𝑅 + ( 𝜁 + √︁ ˜ 𝑧 + ˜ 𝑏 ) 𝜁 ( ˜ 𝑅 + 𝜁 ) / ( ˜ 𝑧 + ˜ 𝑏 ) / , (45) MNRAS , 1–14 (2019) C. Caravita, L. Ciotti and S. Pellegrini − − − M a ss d e n s i t y / ρ n JJE
J2DM
JHD
JMNDM
JLD
JMNDM10 − − − R/r ∗ . . . . . C i r c u l a r v e l o c i t y / √ φ n DMJ 2BH 10 − − − R/r ∗ DMJ MNBH 10 − − − R/r ∗ DMJ MNBH
Figure 2.
Radial profiles in the equatorial plane ( ˜ 𝑧 =
0) of the mass densities R ˜ 𝜌 ∗ (dashed blue), R ˜ 𝜌 ∗ (dotted-dashed magenta), ˜ 𝜌 ∗ (heavy solid), and R h ˜ 𝜌 h (dotted), normalised to 𝜌 n (top row), for the three models of Table 2. In the bottom row, we show the corresponding contributions to the total circular velocity(heavy solid) in the equatorial plane of the mass components, with the additional contribution of the central BH (solid), all normalised to √ 𝜙 n . 𝜙 ∗ ( 𝑅, 𝑧 ) = − 𝜙 n R √︁ ˜ 𝑅 + 𝜁 , 𝜁 = ˜ 𝑎 + √︁ ˜ 𝑧 + ˜ 𝑏 , (46)where R = − R from equation (34). For 𝑎 = 𝑏 = 𝑞 = 𝑎 / 𝑏 the disc flattening parameter. As in JJE models, we add the sphericalNFW halo in equation (44), and a central BH, so that the resultingmulti-component models are completely determined once the valuesof 𝜉 , 𝑞 , R , Υ ∗ , 𝑘 for 𝜌 ∗ , ˜ 𝑏 , 𝑞 , Υ ∗ , 𝑘 for 𝜌 ∗ , 𝜉 h , 𝑐 , R h for 𝜌 h , and R BH for the BH, are assigned, in addition to the total stellarmass 𝑀 ∗ and the scale length 𝑟 ∗ (see Table 2). In Section 4.2.1 weconsider the case of a dominant MN disc, when the ellipsoidal Jaffecomponent can be interpreted as the stellar halo of a disc galaxy, andin Section 4.2.2 the case of a small and counter-rotating stellar discat the center of a dominant stellar spheroid, as sometimes observedin real ETGs (e.g. Morelli et al. 2004; Krajnović et al. 2015; Mitzkuset al. 2017; see also Cappellari 2016). The parameters of the DMhalo and of the central BH are the same as in the JJE model. In the ”Jaffe - Heavy Disc” JHD model (see Table 2), the ellipsoidalJaffe stellar halo 𝜌 ∗ is characterised by a scale-length 𝜉 =
1, aflattening 𝑞 = .
8, a stellar mass fraction of 30% (i.e. R = . Υ ∗ =
6. The dominant MN disc 𝜌 ∗ ( R = . 𝑞 = 𝑏 = .
1, and a lower Υ ∗ = 𝜌 ∗ , ˜ 𝜌 ∗ , and ˜ 𝜌 ∗ , are shown. The resulting isodensity contours of 𝜌 ∗ would be classified as ”disky” near the equatorial plane, and as”boxy” at large distance from the plane. The radial profiles of thedensity distributions (including the DM), in the equatorial plane, areshown in Figure 2. It is apparent how, inside ≈ . 𝑟 ∗ the Jaffe halodominates, around 𝑟 ∗ the MN disc dominates, and 𝜌 h overcomes thetotal 𝜌 ∗ outside ≈ 𝑟 ∗ . Note that, even if R < R , 𝜌 ∗ dominatesthe total density in the central regions, due to the cuspy profile ofthe Jaffe density compared with the flat core of the MN density. Thedensity decomposition reflects on the circular velocity profiles in thebottom panel of the same Figure: the total 𝜐 c at small radii is totallydominated by the BH, and at large radii by the DM halo; while the”bump” around 3 𝑟 ∗ is due to the MN and the DM potentials.The radial profiles in the equatorial plane of the velocity fields,obtained from the Jeans equations, are shown in the middle col-umn of Figure 3, where from top to bottom the total mass- andluminosity-weighted rotational velocity, vertical velocity dispersion,and azimuthal velocity dispersion, are plotted together with the cor-responding quantities for each stellar component separately. For theadopted values of the parameters in Table 2, Δ turns out to be nega-tive in a quite central region, while Δ + 𝜎 is everywhere positive; wedecided to apply the generalised 𝑘 -decomposition in equation (18)to both stellar components, with 𝑘 = . 𝑘 = .
8. The totalvelocity profiles, in the central regions, are completely determinedby the Jaffe profile, because here 𝜌 ∗ > 𝜌 ∗ , compensating also forthe higher Υ ∗ ; in the external regions, instead, the total profiles aredominated by the MN disc. Furthermore, 𝜐 𝜑 stays well below 𝜐 c both in the inner and outer regions (see 𝜐 c in Figure 2), as a clearmanifestation of asymmetric drift in the equatorial plane (e.g. BT08). MNRAS000
8. The totalvelocity profiles, in the central regions, are completely determinedby the Jaffe profile, because here 𝜌 ∗ > 𝜌 ∗ , compensating also forthe higher Υ ∗ ; in the external regions, instead, the total profiles aredominated by the MN disc. Furthermore, 𝜐 𝜑 stays well below 𝜐 c both in the inner and outer regions (see 𝜐 c in Figure 2), as a clearmanifestation of asymmetric drift in the equatorial plane (e.g. BT08). MNRAS000 , 1–14 (2019) eans modeling with multiple stellar populations . . . . . . R o t a t i o n a l v e l o c i t y / √ φ n JJE
J2 0 . . . . . JHD
JMN − . − . − . . . . JLD
JMN0 . . . . . . V e rt i c a l v e l o c i t y d i s p e r s i o n / √ φ n J2 0 . . . . . . . . . . − − − R/r ∗ . . . . . . A z i m u t h a l v e l o c i t y d i s p e r s i o n / √ φ n J2 10 − − − R/r ∗ . . . . . − − − R/r ∗ . . . . . Figure 3.
Radial profiles in the equatorial plane ( ˜ 𝑧 =
0) of the rotational velocities (top row), vertical velocity dispersions (middle row), and azimuthal velocitydispersions (bottom row), normalised to √ 𝜙 n , for the three models of Table 2. Each panel shows the total mass-weighted (heavy solid) and luminosity-weighted(heavy dotted) fields, together with the corresponding fields of 𝜌 ∗ (dashed blue) and 𝜌 ∗ (dotted-dashed magenta). Notice the different values on the verticalscales of the first column (JJE model) and of the top right panel showing the counter-rotation (JLD model). Note that 𝜎 , associated with a flat density profile at the centre, ismuch higher than 𝜎 , associated with 𝜌 ∗ ∼ 𝑅 − in the inner regions,as can be expected from the integration of the vertical Jeans equationfor a power law density distribution in the gravitational field of apoint-mass (i.e. the BH). In addition, Δ of the MN model with thecentral BH vanishes at the centre (see e.g. Chapter 13 in Ciotti 2021),thus in the generalised 𝑘 -decomposition, 𝜐 𝜑 ∼ 𝑘 𝜎 (at variancewith what would happen in the standard Satoh decomposition, i.e. 𝜐 𝜑 = 𝑘 √ Δ ).In the second row of Figure 4, the luminosity-weighted projectedfields 𝜐 los L and 𝜎 los L are shown, and the high rotation of the discis clearly visible. The drop of 𝜐 los L inside 𝑟 ∗ is due to a drop of theintrinsic rotational velocity (Figure 3). Also 𝜎 los L shows the highest values near the equatorial plane, with a toroidal distribution aroundthe centre, and a drop inside 𝑟 ∗ . At variance with the JHD model, in the ”Jaffe - Light Disc” JLDmodel (see Table 2), the ellipsoidal Jaffe distribution 𝜌 ∗ accountsfor almost the whole stellar mass of the galaxy ( R = . 𝜉 = 𝑞 = . Υ ∗ = ) are unchanged. The component 𝜌 ∗ is a small MN disc,with ˜ 𝑏 = .
01 and R = .
04, while 𝑞 =
10 and Υ ∗ = 𝜌 ∗ is MNRAS , 1–14 (2019) C. Caravita, L. Ciotti and S. Pellegrini
Figure 4.
Edge-on projected luminosity-weighted rotational velocity 𝜐 los L (left), and velocity dispersion 𝜎 los L (right), normalised to √ 𝜙 n , for the threemodels of Table 2. Notice the different ranges of values on the colorbars. For the JLD model, the region shown is limited to 𝑟 ∗ to appreciate the central features,in particular the inner counter-rotating disc. The dotted contours show the galaxy isophotes of the surface brightness 𝐼 ∗ , with values spaced by 1 dex. (structurally) identical to that of the JHD model, while ˜ 𝜌 ∗ is muchmore concentrated, so that the total stellar density is distributed inan extended halo with a very small disc. Indeed, the disc is almostinvisible in the last panel, and it would be apparent only with a zoomin, as in Figure 4. The last column of Figure 2 shows the radialprofiles in the equatorial plane of the density components, with theirmass weights, and the resulting decomposition of the galaxy circularvelocity profile. Notice that the central values of 𝜌 ∗ are higher thanthose in the JHD model, due to its smaller size, compensating forthe reduced mass. In the circular velocity plot, this reflects into alarger contribution from the Jaffe component, and a smaller andinner ”bump” of the MN component. As a result, 𝜐 c is almost flatbetween 10 − 𝑟 ∗ and 10 𝑟 ∗ . In the last column of Figure 3, the radial profiles of the velocityfields in the equatorial plane are shown. As in the previous models, ofcourse, the total luminosity-weighted profiles, when distinguishablefrom the mass-weighted ones, are always closer to the profiles of thecomponent with the smaller mass-to-light ratio. For the JLD model,both Δ and Δ are everywhere positive, so we apply the standardSatoh decomposition. The stellar halo is modeled as a slow rotatorwith 𝑘 = .
5, while the circumnuclear stellar disc as a faster andcounter-rotating light disc. In order to have counter-rotation limitedto a central region, we adopt a position-dependent Satoh parameter,defined as follows: 𝑘 ( 𝑅, 𝑧 ) = 𝑘 + ( 𝑘 ∞ − 𝑘 ) 𝑟𝑟 + . 𝑟 ∗ , 𝑟 = √︁ 𝑅 + 𝑧 (47) MNRAS000
5, while the circumnuclear stellar disc as a faster andcounter-rotating light disc. In order to have counter-rotation limitedto a central region, we adopt a position-dependent Satoh parameter,defined as follows: 𝑘 ( 𝑅, 𝑧 ) = 𝑘 + ( 𝑘 ∞ − 𝑘 ) 𝑟𝑟 + . 𝑟 ∗ , 𝑟 = √︁ 𝑅 + 𝑧 (47) MNRAS000 , 1–14 (2019) eans modeling with multiple stellar populations (e.g. Ciotti et al. 2021, in preparation; see also Negri et al. 2014a foran alternative parametrisation), with 𝑘 = − . 𝑘 ∞ = .
1, where thenegative sign of 𝑘 assures the counter-rotation of the disc, as canbe seen in the top right panel of Figure 3. At very small radii (inside10 − 𝑟 ∗ ), the total rotational velocity is positive because the densityis dominated by the Jaffe component. We stress that the module of 𝜐 𝜑 decreases towards the centre, at variance with the JHD model,because now 𝜐 𝜑 = 𝑘 Δ , and Δ →
0, as explained in the previousSection. The central total vertical velocity dispersion is higher thanthat of the JHD model, even if the Jaffe component is structurallyidentical, because of the higher R and of the more concentrated MNdisc.In the last row of Figure 4, the los luminosity-weighted velocitiesshow clearly the effect of the inner thin disc; the region shown islimited to 𝑟 ∗ to appreciate the central features. In particular, in the 𝜐 los L distribution we have counter-rotation at small radii (but not inthe very centre). The disc is also responsible for the highest valuesof the 𝜎 los L in the equatorial plane, and the extended surroundingtoroidal distribution is also present, in analogy with the JHD model. We presented the theoretical framework, and its numerical imple-mentation in the code JASMINE2, for an efficient Jeans modeling ofmulti-component axisymmetric galaxies. The models can include anarbitrary number of stellar components, with different structural, dy-namical and population properties, a DM halo, and a central BH. Thestructural and dynamical properties of the single stellar componentscan be mass- or luminosity-weighted, and projected on the plane ofthe sky. The internal dynamics of each stellar component is implicitlydescribed by a 2-integral distribution function (in general differentfor each component); a decomposition of the azimuthal velocity fieldmust then be chosen. For each component, the code can work withthe Satoh (1980) 𝑘 -decomposition, where 𝜐 𝜑𝑖 = 𝑘 𝑖 √ Δ 𝑖 , and/or witha generalised 𝑘 -decomposition, where 𝜐 𝜑𝑖 = 𝑘 𝑖 √︃ Δ 𝑖 + 𝜎 𝑖 ; further-more, the parameter 𝑘 𝑖 can be constant or position-dependent. Thegeneralised decomposition allows for the modeling of systems with Δ 𝑖 < family of models. The scaledJeans solutions are then combined in PP, with the desired mass andluminosity weights, and through the choice of appropriate kinemat-ical decompositions, and then projected. The PP procedure can beperformed several times, obtaining different specific models in thesame family. Finally, for each model, the two physical scales 𝑀 ∗ and 𝑟 ∗ can be assigned. A further benefit of the presented approach isthe possibility to gain the most transparent understanding of the roleof each density component in determining the resulting kinematicalfields of the galaxy. In order to illustrate the features of the code, we presented threegalaxy models, composed of two stellar components, a sphericalNFW DM halo and a central supermassive BH. In the JJE model, thetotal spherical stellar profile and one ellipsoidal stellar component,both Jaffe models, are assigned; the second stellar component isgiven by their difference. This model, when the DM halo is absent,has several properties available in analytical form, and thus has beenused to test the code. The JHD model consists of a large and massiveMN stellar disc, coupled with an ellipsoidal Jaffe stellar model, thatcan be seen as the stellar halo of a disc galaxy. In the JLD model, anellipsoidal Jaffe component dominates in mass, and a MN stellar discis small, inner and counter-rotating, as sometimes found in early-typegalaxies.Ongoing applications of JASMINE2 include the building of multi-component galaxy models for numerical simulations of gas flows ingalaxies (e.g. Negri et al. 2014b; Gan et al. 2019a,b, Ciotti et al.2021 in preparation); the study of circumnuclear stellar discs (alsowith counter-rotation, see e.g. Morelli et al. 2004; Krajnović et al.2015; Mitzkus et al. 2017; Sormani et al. 2020; see also Cappellari2016); a systematic exploration of galaxy models constrained to lie onthe major observed Scaling Laws, extending the statistical approachpioneered in Bertin et al. (2002) and Lanzoni & Ciotti (2003). ACKNOWLEDGEMENTS
We thank Antonio Mancino for independent checks of the results ofJJE models.
DATA AVAILABILITY
No new data were generated or analysed in support of this research.
REFERENCES
Bertin G., 2014, Dynamics of GalaxiesBertin G., Ciotti L., Del Principe M., 2002, A&A, 386, 149Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. PrincetonUniversity PressCappellari M., 2008, MNRAS, 390, 71Cappellari M., 2016, ARAA, 54, 597Ciotti L., 2000, Lecture notes on stellar dynamicsCiotti L., 2021, Introduction to Stellar DynamicsCiotti L., Bertin G., 2005, A&A, 437, 419Ciotti L., Pellegrini S., 1996, MNRAS, 279, 240Ciotti L., Mancino A., Pellegrini S., Ziaee Lorzad A., 2021, MNRAS, 500,1054Gan Z., Ciotti L., Ostriker J. P., Yuan F., 2019a, ApJ, 872, 167Gan Z., Choi E., Ostriker J. P., Ciotti L., Pellegrini S., 2019b, ApJ, 875, 109Jaffe W., 1983, MNRAS, 202, 995Kormendy J., Ho L. C., 2013, ARA&A, 51, 511Krajnović D., et al., 2015, MNRAS, 452, 2Kuzmin G. G., 1956, Azh, 33, 27Lanzoni B., Ciotti L., 2003, A&A, 404, 819Maraston C., 2005, MNRAS, 362, 799Mitzkus M., Cappellari M., Walcher C. J., 2017, MNRAS, 464, 4789Miyamoto M., Nagai R., 1975, PASJ, 27, 533Morelli L., et al., 2004, MNRAS, 354, 753Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563Negri A., Ciotti L., Pellegrini S., 2014a, MNRAS, 439, 823Negri A., Posacki S., Pellegrini S., Ciotti L., 2014b, MNRAS, 445, 1351Plummer H. C., 1911, MNRAS, 71, 460Posacki S., Pellegrini S., Ciotti L., 2013, MNRAS, 433, 2259MNRAS , 1–14 (2019) C. Caravita, L. Ciotti and S. Pellegrini
Renzini A., Buzzoni A., 1986, in Chiosi C., Renzini A., eds, Astrophysicsand Space Science Library Vol. 122, Spectral Evolution of Galaxies. pp195–231, doi:10.1007/978-94-009-4598-2_19Satoh C., 1980, PASJ, 32, 41Smet C. O., Posacki S., Ciotti L., 2015, MNRAS, 448, 2921Sormani M. C., Magorrian J., Nogueras-Lara F., Neumayer N., Schönrich R.,Klessen R. S., Mastrobuono-Battisti A., 2020, MNRAS, 499, 7Yoon D., Yuan F., Ostriker J. P., Ciotti L., Zhu B., 2019, ApJ, 885, 16
APPENDIX A: POSITIVITY CONDITION FOR JJEMODELS
The stellar component 𝜌 ∗ , for the 2-component ellipsoidal modelsdescribed in Section 4.1, is given by the difference of an assignedtotal 𝜌 ∗ and an assigned 𝜌 ∗ . This approach naturally leads to discussthe positivity of 𝜌 ∗ , with a treatment similar to that followed in theappendix of CMPZ21, and references therein. We generalise equation(43) as 𝜌 ∗ ( 𝑅, 𝑧 )( − 𝛾 ) 𝜌 n = 𝜉𝑞𝑚 𝛾 ( 𝜉 + 𝑚 ) − 𝛾 − R 𝜉 𝑞 𝑚 𝛾 ( 𝜉 + 𝑚 ) − 𝛾 , (A1)recovering the case of JJE models for 𝛾 =
2. In order to discuss thepositivity condition for 𝜌 ∗ , we use spherical coordinates, so that ( 𝑅, 𝑧 ) = 𝑟 ( sin 𝜃, cos 𝜃 ) and 𝑚 = 𝑠 Ω , 𝑚 = 𝑠 Ω , 𝑠 ≡ 𝑟𝑟 ∗ , (A2)where Ω ≡ sin 𝜃 + cos 𝜃𝑞 , Ω ≡ sin 𝜃 + cos 𝜃𝑞 . (A3)The positivity of 𝜌 ∗ reduces to a condition on R , given by R ≤ R M ≡ inf I (cid:34) 𝜉𝑞 𝜉 𝑞 (cid:18) Ω Ω (cid:19) 𝛾 (cid:18) 𝜉 + 𝑠 Ω 𝜉 + 𝑠 Ω (cid:19) − 𝛾 (cid:35) , (A4)over the rectangular region I ≡ { 𝑠 ≥ , ≤ 𝜃 ≤ 𝜋 / } in the ( 𝑠, 𝜃 ) plane. Following the discussion in CMPZ21, we determine R M = min (cid:16) R c , R ∞ , R , R 𝜋 / , R int (cid:17) , (A5)where the first four quantities refer to the minimum value of the r.h.s.of equation (A4) over the boundaries of I , and R int is the value of aminimum (if it exists) in the interior of I . When 𝑞 ≠ 𝑞 , it is simpleto show that no critical points can exist in the interior of I , and so thediscussion reduces to the boundaries of I : geometrically, R M can bereached only at the centre ( 𝑠 = R c ), at infinity ( 𝑠 → ∞ , R ∞ ), alongthe symmetry axis ( 𝜃 = R ), or on the equatorial plane ( 𝜃 = 𝜋 / R 𝜋 / ).We begin with R c and R ∞ , obtaining R c = 𝜉 − 𝛾 𝑞 𝜉 − 𝛾 𝑞 min ≤ 𝜃 ≤ 𝜋 / (cid:18) Ω Ω (cid:19) 𝛾 , (A6) R ∞ = 𝜉𝑞 𝜉 𝑞 min ≤ 𝜃 ≤ 𝜋 / (cid:18) Ω Ω (cid:19) . (A7)Now, from equation (A3), it is easy to show that for a generic 𝛼 ≥ ( Ω / Ω ) 𝛼 reaches its minimum at 𝜃 = 𝜋 / 𝑞 ≤ 𝑞 ,and at 𝜃 = 𝑞 ≤ 𝑞 , so thatmin ≤ 𝜃 ≤ 𝜋 / (cid:18) Ω Ω (cid:19) 𝛼 = , 𝑞 ≤ 𝑞, (cid:18) 𝑞𝑞 (cid:19) 𝛼 , 𝑞 ≤ 𝑞 , (A8) and the conditions in equations (A6) and (A7) can be finally sum-marised as R c = 𝜉 − 𝛾 𝑞 𝜉 − 𝛾 𝑞 min (cid:32) , 𝑞 𝛾 𝑞 𝛾 (cid:33) , (A9) R ∞ = 𝜉𝑞 𝜉 𝑞 min (cid:32) , 𝑞 𝑞 (cid:33) . (A10)Along the symmetry axis, and in the equatorial plane, condition (A4)becomes R = 𝜉𝑞 𝜉 𝑞 inf ≤ 𝑠< ∞ (cid:18) 𝜉 𝑞 + 𝑠𝜉𝑞 + 𝑠 (cid:19) − 𝛾 , (A11) R 𝜋 / = 𝜉𝑞 𝜉 𝑞 inf ≤ 𝑠< ∞ (cid:18) 𝜉 + 𝑠𝜉 + 𝑠 (cid:19) − 𝛾 , (A12)and simple algebra finally shows that the results can be summarisedas R = 𝜉𝑞 𝜉 𝑞 min (cid:34) , (cid:18) 𝜉 𝑞 𝜉𝑞 (cid:19) − 𝛾 (cid:35) , (A13) R 𝜋 / = 𝜉𝑞 𝜉 𝑞 min (cid:169)(cid:173)(cid:171) , 𝜉 − 𝛾 𝜉 − 𝛾 (cid:170)(cid:174)(cid:172) . (A14)For the JJE models in Section 4.1, with 𝜉 < 𝜉 , 𝑞 < 𝑞 , and 𝛾 = R ≤ R M = 𝜉 𝑞 𝜉𝑞 . (A15) This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000