Anisotropy ansatz for the axisymmetric Jeans equations

The Jeans equations do not form a closed system, and to solve them a parametrization relating the velocity moments is often adopted. For axisymmetric models, a phenomenological choice (the"$b$-ansatz") is widely used for the relation between the vertical ($\sigma_z^2$) and radial ($\sigma_R^2$) components of the velocity dispersion tensor, thus breaking their identity present in two-integral systems. However, the way in which the ansatz affects the resulting kinematical fields can be quite complicated, so that the analysis of these fields is usually performed only after numerically computing them. We present here a general procedure to study the properties of the ansatz-dependent fields $\overline{v_{\varphi}^2}$, $\Delta = \overline{v_{\varphi}^2} - \sigma_z^2$ and $\Delta_R = \overline{v_{\varphi}^2} - \sigma_R^2$. Specifically, the effects of the $b$-ansatz can be determined before solving the Jeans equations once the behaviour over the ($R,z$)-plane of three easy-to-build ansatz-independent functions is known. The procedure also constrains the ansatz to exclude unphysical results (as a negative $\overline{v_{\varphi}^2}$). The method is illustrated by discussing the cases of three well-known galaxy models: the Miyamoto&Nagai and Satoh disks, and the Binney logarithmic halo, for which the regions and the constraints on the ansatz values can be determined analytically; a two-component (Miyamoto&Nagai plus logarithmic halo) model is also discussed.


INTRODUCTION
The Jeans equations (hereafter JEs) are one of the standard tools for the modelling of stellar systems (e.g.Binney & Tremaine 2008, hereafter BT08; see also Ciotti 2021, hereafter C21) and to extract information from observations (e.g.Scott et al. 2015, Cappellari 2016, Li et al. 2016, Zhu et al. 2023).However, the JEs are moments of the more fundamental Boltzmann equation, and in the collisionless limit they suffer in general from a "closure problem".Their closure can be obtained by assuming a (more or less motivated) dependence of the phase-space distribution function (hereafter DF) on the available integrals of motion.Alternatively, some relation between the velocity moments (usually between the velocity dispersion tensor components) can be imposed through some phenomenological "ansatz".This latter approach is not as elegant and physically sound as deriving all the model properties from the DF, coupled with the Poisson equation and asking for self-consistency (e.g.King 1963, Stiavelli & Bertin 1985, Bertin & Varri 2008), or using a DF built from actions (e.g.Binney 2010; see also Posti et al. 2015, Binney & Vasiliev 2023), or reconstructing the DF numerically with the Schwarzschild orbital superposition method (Schwarzschild 1979; see also, e.g., Statler 1987, Cappellari et al. 2007, Thomas et al. 2009).However, the methods mentioned above are not always well suited for exploratory works, because in general the solution of the self-consistency problem requires non-trivial numerical work.
An important advantage of using the JEs is that the stellar density distribution of the model can be chosen at the beginning to be in ★ E-mail: leonardo.dedeo2@unibo.itgood agreement with the observational data (for example by using specific density profiles, or by multi-component modelling), and that the closure ansatz guarantees some direct control on the resulting kinematics.Of course, the problem of the existence of a non-negative DF for the obtained model remains in general open.
In this work, we address two aspects worth a thorough investigation, both related to how a widely used closure ansatz (the -ansatz, Cappellari 2008) affects the solutions of the JEs: the first is to formalize a general procedure, to be applied before solving the equations, to determine the constraints on this ansatz, as for example those coming from the request of positivity of  2  , thus avoiding repeated and time-consuming numerical tests.The second is to gain some qualitative understanding of its effects on the kinematical fields, before their construction.Moreover, as a widely used decomposition of  2  (e.g. that proposed by Satoh 1980, and its variants) requires the positivity of certain functions, we also determine the conditions for its applicability.
The paper is organized as follows: in Sections 2 and 3, after reviewing the JEs for two-integral axisymmetric systems, the general solutions with the -ansatz, relating the vertical ( 2  ) and the radial ( 2  ) velocity dispersions, are derived, and cast in a form suitable for the successive investigation.In Section 4, we determine the constraints on the ansatz to assure the positivity of  2  and to use the Satoh decomposition and its generalizations.In Sections 5 and 6, we apply the procedure to some well-known galaxy models, for which a fully analytical treatment is possible: the Miyamoto & Nagai and Satoh disks, and the Binney logarithmic halo.A more realistic two-component model, made by a stellar Miyamoto & Nagai disk embedded in a dark matter Binney halo, is also discussed.In Section 7, the main results are summarized.

TWO-INTEGRAL SYSTEMS
As well known, the JEs for an axisymmetric stellar system described by a two-integral DF  =  (,   ), where  and   are respectively the orbital energy and the axial angular momentum, are (see e.g.BT08; C21):  * (, ) and Φ(, ) are the stellar density distribution and the total (e.g., stars, plus dark matter) gravitational potential.Standard cylindrical coordinates are used,  2  is the vertical velocity dispersion, and a bar over a symbol indicates the average over velocity in phase-space.In particular,  2  =   2 +  2  , where   is the streaming velocity field in the azimuthal direction, and  2  is the azimuthal velocity dispersion.Finally, the only non-vanishing ordered velocity field is   , while   =   = 0; moreover,  2  =  2  , and all the mixed components of the velocity dispersion tensor vanish.
The solutions of equation ( 1) with null boundary conditions at infinity are where is a commutator-like operator (see e.g.Barnabè et al. 2006; see also C21 and references therein).From the second of equation ( 2), Δ can be obtained or by differentiation of  *  2  , or by integration of the commutator, as in equation ( 3).This latter approach is to be preferred, as it usually reveals important properties of the solutions that are not apparent in the first approach based on differentiation: for example, it is immediate to show that the commutator vanishes for a spherical density  * () in a spherical total potential Φ() and that, in the case of ellipsoidal densities in ellipsoidal potentials, the commutator is everywhere positive (negative) when the density shape is flatter (rounder) than the potential (see e.g.C21).Of course, when using the commutator for the computation of Δ, the radial derivative of  *  2  is obtained as a byproduct: an identity we will use in the following.Once  2  and Δ are known, one has where the second identity involving again a commutator can be easily proved from equations ( 2) and (3).
Notice that a model with  2  < 0 somewhere is certainly physically inconsistent, but  2  ≥ 0 everywhere is not a sufficient condition for consistency: as well known, there are models with "acceptable" solutions of the JEs and a negative (unphysical) DF (see e.g.Ciotti & Pellegrini 1992; see also Chapter 14 in C21, and references therein).

The Satoh 𝑘-decomposition
For axisymmetric models with Δ ≥ 0 everywhere, Satoh (1980) introduced the widely used -decomposition of  2  : where  is constant with 0 ≤  2 ≤ 1.If  = 0, then   = 0 and no net rotation is present, while, if  2 = 1, then  2  =  2  and the system is an isotropic rotator.Moreover, if Δ = 0 (as for spherical models), then no ordered rotation is possible, and the system is isotropic independently of .
As shown in Ciotti & Pellegrini (1996), the original Satoh decomposition can be easily generalized to assume a spatially dependent  (, ), provided that where the upper limit  2 M corresponds to maximally rotating models with no net velocity dispersion in the azimuthal direction, and is then obtained from equation (6) imposing  2  = 0 everywhere; in this case, the density flattening is fully supported by the streaming velocity field   .Of course,  is not required to be positive, thus allowing the modelling of counter-rotating structures with negative   (see e.g.Negri et al. 2013).
In Appendix A1, we discuss the most general decomposition for  2  , which holds also for models with Δ < 0.

MORE GENERAL SYSTEMS
Having assessed the classical case of two-integral axisymmetric systems, we now turn to the focus of the paper, i.e. the study of the properties of the kinematical fields associated with more general ansatz distinguishing between  2  and  2  , and so implicitly based on a DF with a third integral in addition to  and   .If the third integral is an even function of   and   , the first of equation ( 1) remains unchanged, while the second becomes (Cappellari 2008, C21): Notice that, in addition to  2  , all quantities depending on  * and Φ, such as the commutator [ * , Φ] and the function  in equation ( 4), remain the same as in the two-integral case.
Different ansatz can be introduced to solve equation ( 8), which is independent of the vertical velocity dispersion, and contains the two unknown functions  2  and  2  .In the following, we study in detail the "-ansatz", relating  2  with  2  .Introduced by Cappellari (2008), it was adopted for the Jeans Anisotropic Modelling method (JAM; see also Cappellari 2020), that is widely used to reproduce the properties of observed galaxies (e.g., Cappellari et al. 2013, Zhu et al. 2016, Loubser et al. 2020, Nitschai et al. 2021, Surti et al. 2024).This solution of the JEs, implying the alignment of the velocity ellipsoid with the cylindrical coordinates, was presented as an efficient modelling able to capture the main properties of the velocity ellipsoid inferred from extensive three-integrals Schwarzschild's modelling of integral-field stellar kinematics, under the mass-follows-light hypothesis (Cappellari 2008).The main motivation for the adoption of the "-ansatz", in fact, is that it allows to model adequately the observations, in particular the integral-field spectroscopy of axisymmetric galaxies classified as regular rotators and stellar disks (see Sections 2.3 and 2.4 of Cappellari 2008;Cappellari 2016).Of course, on the theoretical side, the -ansatz is not the only one possible, and in Appendix A2 we present the "-ansatz", leading to a nice closure of equation ( 8).

The b-ansatz
In the "-ansatz", the unknown  2  is linked to  2  through the choice of the function (, ) ≥ 0, as When  = 0, the system has no radial velocity dispersion, while  = 1 gives the two-integral case.
Inserting equation ( 9) in equation ( 8), recalling the definition of  in equation ( 4), and solving for Δ  , we obtain: where in the last expression we have used again equation (4).Note how  multiplies functions that are independent of the adopted ansatz.Once Δ  is known, if we restrict to a  function that depends only on  (or to a constant , a special case commonly used), / = 0 on the l.h.s. of equation ( 10), and where from equation ( 4) Only () functions leading to  2  ≥ 0 everywhere are physically acceptable, and Section 4.1 deals with this request.
Finally, from equation (11), we obtain: where from equations ( 4) and ( 12) From equation ( 13), if Δ  ≥ 0 and () ≥ 1, then Δ ≥ 0, and if Δ  < 0 and () ≤ 1, then Δ < 0. A complete discussion on the conditions for the positivity of Δ  and Δ is given in Sections 4.2 and 4.3.Notice that the values of Δ in the -ansatz are not the same of the two-integral case, because, with the introduction of the ansatz,  2  remains unaltered, but  2  changes.Also important, the functions ,  and  are ansatz-independent.
A different decomposition for the azimuthal motions can be given by choosing the parameter , introduced in Cappellari et al. (2007): where clearly with the first inequality required for   2 ≥ 0. Positive  values correspond to radial anisotropy, while  < 0 gives tangential anisotropy; in two-integral systems,  = 0 corresponds to the isotropic case.Notice that, from equation ( 17),  can be either positive or negative if Δ  > 0; instead, if Δ  < 0,  cannot be taken negative, i.e. the velocity dispersion tensor can only be radially anisotropic.In Section 7, some findings based on the use of this decomposition to interpret recent observations (Wang et al. 2021) are discussed in light of our analysis.

PHYSICAL CONSTRAINTS ON THE ANSATZ
Equation ( 8) can be solved only with the introduction of some closure, as those presented in Section 3.1 or Appendix A2.However, arbitrary ansatz functions can lead to unphysical solutions of the JEs, such as negative values of  2  .Therefore, it would be useful to know in advance (i.e., before solving the equations) what constraints must be imposed on the ansatz function in order to avoid unphysical solutions, and also how specific choices of the ansatz functions affect the properties of the solutions in different regions of space.This Section is dedicated to these problems.
In practice, being  2  ≥ 0 independent of the ansatz, the positivity of  2  is guaranteed by equation ( 9); therefore, we focus on the positivity of  2  only.Information on the sign of Δ and Δ  is also relevant for the modelling because, although not directly related to model consistency, their positivity is needed to apply a decomposition of  2  as those in equations ( 6) and ( 15).In the following, we restrict for simplicity to models with Φ/ ≥ 0 everywhere (the most common situation), and we limit to consider  = (); finally, we indicate with P = (, ) a generic point in the (, )-plane.

Positivity of 𝑣 2
: the B region We discuss here the positivity of  2  , a condition necessary to model consistency.Introducing the sets equation ( 11) shows that, independently of (),  2  ≥ 0 over the region B +0 = B + ∪ B 0 .Instead, the behaviour of  over B − constrains the possible choices of (), as we now describe.We indicate with Pr(B ± ) the projection of B ± on the -axis, and with B ±  =  :  ≷ 0,  ∈ Pr(B ± ) the radial section at fixed  of B ± (see Figure 1 for a qualitative illustration).As we restrict to systems with a reflection symmetry about the equatorial plane, without loss of generality, in the following we limit the discussion to  ≥ 0. From equation ( 11), the condition where the last equality has been obtained using equation ( 14).Notice that the independence of  from  implies that the upper limit  M (), determined over B −  , applies also to points in the complementary section B +  : therefore, the function  M () provides a constraint over the whole rectangular strip defined by  ≥ 0 and  ∈ Pr(B − ).From now on, we will refer to this region as to "the strip Pr(B − )".In the special case of a constant , the condition in equation ( 19) reduces to  ≤  M , where  M is the minimum of  M ().

Positivity of Δ 𝑅 : the D region
As seen in Section 3.1, the necessary condition to apply the decomposition to  2  is to have Δ  ≥ 0 or Δ ≥ 0. Introducing the sets where  is defined in equation ( 4), the first of equation ( 10) (where now / = 0) shows that, independently of (), Δ  ≥ 0 over the region where Pr(D ± ) is the projection of D ± on the -axis, D ±  =  :  ≷ 0,  ∈ Pr(D ± ) , and the last equality follows from equation (4).Due to the independence of  from , the condition in equation ( 21) must be verified over the whole strip Pr(D − ).In the special case of a constant , then Δ  ≥ 0 for  ≤  0 , where  0 is the minimum of  0 ().
Notice that from Theorem B.1 the strip Pr(B − ) is contained in the strip Pr(D − ), and in this common region  0 () ≤  M (), as proved in Theorem B.4.

Positivity of Δ: the C region
The determination of the positivity of Δ as a function of () is more complicated, depending on the sign of both the  and  functions in equation ( 13).As done above, we introduce the sets and The sign of Δ at a point P depends on its position in the (, )-plane, as follows: and over For a spatially constant , the previous inequalities hold replacing  1 () and  2 () with  1 and  2 , which are respectively the maximum of  1 (), and the minimum of  2 ().

Some general considerations
The properties of the kinematical fields resulting from a given choice for () are determined by the relations between the ansatzindependent regions B ± , C ± and D ± and the functions defined over them.
In this case, from the considerations above, B − and C − are disjoint, but this does not mean that the two projections in equation ( 25) are disjoint too: actually, they may coincide with the -axis.In the special case of [ * , Φ] = 0 everywhere, in equation ( 25) we would have  1 =  0 =  2 = 1.We conclude with an example of the qualitative effects of an increase of  (for simplicity assumed as spatially constant) on  2  ,  2  , Δ  and Δ, in the various regions.Suppose a given model is assigned, and the regions B ± , C ± and D ± have been determined; furthermore, suppose [ * , Φ] ≥ 0 everywhere, so that equation ( 25) holds.We start considering the case of  = 0 (i.e. 2  = 0 everywhere): as 0 is smaller than  0 and  M , then  2  and Δ  are positive everywhere, while Δ < 0 somewhere in C − (because 0 <  1 ), and Δ ≥ 0 everywhere in C + (because 0 <  2 ).For increasing ,  2  increases everywhere,  2  decreases over B − and increases over B + , while Δ  decreases over D − and increases over D + .As B − ⊆ D − , over B − one has that Δ  decreases because  2  decreases and  2  increases; over the remaining part of D − (i.e.B + ∩D − ), instead, Δ  decreases, and so  2  increases more than  2  .Increasing further  above  1 , Δ becomes positive over C − , then when  >  0 , Δ  becomes negative somewhere over D − , and for  >  2 , Δ becomes negative somewhere over B − ; finally, for  >  M ,  2  < 0 somewhere in B − , and the model becomes unphysical.These behaviours are summarized in the table in Figure 2.

ONE-COMPONENT MODELS
In light of the results in Section 4, we now explore the behaviour of some well-known galaxy models that allow for an almost complete analytical treatment, i.e. the Miyamoto & Nagai (1975) (hereafter MN) disk, the Satoh (1980) disk, and the Binney logarithmic halo (BT08).In Section 6, we then consider the case of the two-component model made by a stellar MN disk embedded in a dark matter Binney logarithmic halo.

The Miyamoto & Nagai disk
The potential-density pair of the MN disk of total mass  * , and scale lengths  * and  * , can be written as: where  ≡ √︁ 1 +  2 , and  ≡  * / * measures the flattening of the disc.For  * = 0, the MN model reduces to the Plummer (1911) sphere, and for  * = 0 to the razor-thin Kuzmin disc (Kuzmin 1956, Toomre 1963).In the formulae above and in this Section,  and  are assumed to be normalized to  * ≠ 0.
From equations ( 2) and (3), we have: In particular, [ * , Φ * ] ≥ 0 everywhere, vanishing in the spherical case ( = 0): therefore, equation ( 25) and the considerations made at the end of Section 4 apply.For reference, in the top panels of Figure 3, we show the 2D maps of the ansatz-independent fields  2  ,  [ * , Φ * ] / * (i.e. the Δ field in the two-integral case), and Φ/, for a disc flattening of  = 1.The first task to explore the model behaviour is the identification of the B, C and D regions.Quite remarkably, equation ( 27) allows for an analytical expression of the ,  and  functions: i.e.B − is the portion of the (, )-plane at the right of the hyperbola B 0 with vertex  = (1 + )/ √ 5 on the equatorial plane, and asymptotes  = ± √ 5, so that Pr(B − ) coincides with the -axis.In the top panels of Figure 4, the red line shows B 0 for three representative values of .Moreover, the radial derivative of the first of equation ( 27) shows that  < 0 over the whole space, and  = 0 for  = 0, i.e.D 0 coincides with the -axis; therefore, Pr(D − ) coincides with the -axis, and D + is empty.Finally, from the second of equation ( 13), we get The biquadratic in  above has a positive discriminant, and a permanence and a variation of signs in the coefficients: it follows that C − is the portion of the (, )-plane contained between the -axis and the (positive) square root of the (positive) solution of the biquadratic.We do not report here the expression for C 0 , which however can be determined without difficulty; Pr(C − ) coincides with the -axis.
In the top panels of Figure 4, the green line shows C 0 for three representative values of .Again from Figure 4, it is apparent that B − ⊆ C + and C − ⊆ B + , as expected from the general discussion in Section 4; notice also how the separation between B − and C − (i.e. the region B + ∩ C + ) becomes larger as  increases.Moreover, in the spherical case ( = 0) we would obtain B 0 = C 0 , B − = C + and B + = C − .
We now determine the constraints on () for the positivity of  2  , Δ  and Δ. Remarkably, all the computations in equations ( 19), ( 21), ( 23) and ( 24) can be carried out explicitly, and we obtain In the top-left panel of Figure 5, we show the four functions for  = 1, where the chain of inequalities in equation ( 25) is apparent, considering that for the MN disk Pr(B − ) and Pr(C − ) both coincide with the -axis.If one restricts to  functions independent of , we have independent of : in particular, this holds for the spherical case (the Plummer sphere).
Figure 6 illustrates the behaviour of  2  , Δ  and Δ for the MN disk with  = 1, and for three values of , i.e.  = 0.5, 1, 2. The adopted values map three different kinematical configurations, as can be seen from equation (33).For  = 0.5,  2  and Δ  are everywhere positive, while Δ is expected to be negative over some region in B + ∩ C − = C − , being 0.5 <  1 , and positive over B − ∩ C + = B − , being 0.5 <  2 .This is nicely confirmed by the three left panels in Figure 6, where white regions indicate negative values of the fields.The three central panels correspond to  = 1, i.e. to the twointegral solutions (see also the three top panels in Figure 3).The fields in this case are everywhere positive, in agreement with the general discussion and equation ( 33).Note how the increase of  from 0.5 to 1 produces the expected changes summarized in Figure 2, i.e.  2  decreases over B − and increases over B + , Δ  decreases everywhere being D − for this model coincident with the (, )-plane, and Δ behaves qualitatively as  2  , becoming non-negative over C − because 1 =  1 .Finally, the three right panels of Figure 6 correspond to the unphysical model with  >  M , and a large region in the  2  map becomes white: increasing further  would increase the size of the white region up to the whole B − .However,  2  continues to increase in B + .As 2 >  0 , Δ  is now negative over a large portion of the (, )-plane and, increasing further , Δ  would become negative everywhere.Finally, the negative region of Δ has now switched to the B − region, being 2 >  2 .Interestingly, the shape and the position of the white region of  2  resemble those where the numerically evaluated  2  for the similar model of the Satoh disk with  = 4/3 (Cappellari 2020, Figure 10) goes to zero.

The Satoh disk
The potential-density pair of the Satoh (1980) disk of total mass  * , and scale lengths  * and  * , can be written as: where again  ≡ √︁ 1 +  2 ,  ≡  * / * measures the flattening of the disk, and all lengths are assumed to be normalized to  * ≠ 0 as throughout the Section.Notice that, at variance with the MN model, the spherical limit of the Satoh model ( = 0) reduces to the point mass case of no practical interest.From equations (2) and (3), one has: Similarly to the MN disk, [ * , Φ * ] ≥ 0 everywhere, and also the maps of  2  , [ * , Φ] and Φ/ are very similar to those of the MN model in Figure 3, and thus are not shown.
The B, C and D regions can be determined analytically.In particular, we have so that B − is the portion of the (, )-plane to the right of the curve B 0 with vertex coordinates  = 0 and  = √︁ ( + 2)/5, and asymptotes  = ± √ 5; Pr(B − ) coincides with the -axis.The radial derivative of the first of equation ( 35) shows that  < 0 everywhere, except for the -axis, where  = 0; therefore, as for the MN disk, Pr(D − ) coincides with the -axis, and the set D + is empty.Finally, from the second of equation ( 13),  (37) and again, similarly to the MN disk, the resulting expression is a biquadratic with two real solutions.From the fact that  ≥ 1, the coefficients present a permanence and a variation of their sign: therefore, C − is the region between the -axis and the (positive) square root of the (positive) solution of the biquadratic.The explicit expression of C 0 can be obtained without difficulty; we just notice that, as for the MN disk, Pr(C − ) coincides with the -axis.Again, all the results summarized in Section 4 about the relative position of the regions apply, i.e.B − ⊆ D − , and The limitations on () for the positivity of  2  , Δ  and Δ can all be obtained analytically.From equations ( 19), ( 21), ( 23) and ( 24), where  * √︁ ( + ) 2 − 1.In the top-right panel of Figure 5, we show the four functions above for  = 1; again, they fulfil the inequalities in equation ( 25), considering that Pr(B − ) and Pr(C − ) both coincide with the -axis.For a spatially constant , we have independently of .
We do not show the analogous of Figure 6, due to the close similarity, for sufficiently flat models, with the MN case, so that all comments made for the MN disk apply.The only noticeable difference is that now  1 ≠  0 ≠  2 , while for the MN disk they are the same [see equation ( 33

The Binney logarithmic halo
The potential-density pair of the Binney logarithmic halo (BT08) of asymptotic circular velocity  h , scale length  h , and potential flattening  is: where throughout this Section the lengths are assumed normalized to  h ≠ 0. As well known, for  < 1/ √ 2,  h becomes negative on the -axis.The dynamical properties of this model are given by Evans (1993) (see also Evans 1994 andEvans &de Zeeuw 1994 for the properties of the larger family of the so-called power-law models); here we just use equations ( 2) and ( 3) to obtain it vanishes for  = 1 (the spherical limit), and it is ≤ 0 for  > 1.We finally notice that, in the two-integral case, showing that 2  becomes negative for  > √ 2, and the models are hence inconsistent.Actually, Evans (1993) proved that the twointegral DF becomes negative for  ≳ 1.08 so that models with 1.08 ≲  ≤ √ 2 do not exist 2 although they have non-negative  2  and  2  .From now on, we consider models with 1/ √ 2 ≤  ≤ 1, so that [ h , Φ h ] ≥ 0 everywhere: in the central panels of Figure 3, the 2D maps of  2  , [ h , Φ h ], and Φ/ are shown for  = 0.75.As for the two previous models, the B, C and D regions can be determined in a fully analytical way.We notice however that for  ≃ 1 the Binney halo is qualitatively different from a disk, resembling a spheroid; therefore, we expect some important differences compared to the previously discussed cases.From equation ( 43), we have With some work, it can be proved that the discriminant of the biquadratic is strictly positive independently of ; moreover, for  ≥ 1/ √ 2, the last term in equation ( 45) is negative, so that a permanence and a variation of the sign of its coefficients are present.We conclude that B 0 is the (positive) square root of the (positive) solution, and B − is the region of the (, )-plane to the right of B 0 .As for the previous models, Pr(B − ) coincides with the -axis; for simplicity, we do not report the explicit expression of B 0 here, but we show it as the red line in the bottom panels of Figure 4, for three representative values of .The properties of the D region are more complicated compared to the previous cases: in fact, while for  ≥ √︁ 2/3,  < 0 everywhere (vanishing on the -axis) and thus D − fully covers the (, )-plane, for  < √︁ 2/3 and the D + region now exists: for example, in the bottom right panel of Figure 4, the dashed line shows D 0 for  = 0.75, and D + is the part of the plane above it.Notice how, independently of , Pr(D − ) coincides with the -axis.Finally, from equation ( 13), where, in the considered range of 1/ √ 2 ≤  ≤ 1, the discriminant of the biquadratic is positive, and the coefficients present a permanence and a variation of the sign.Therefore, C 0 is the (positive) square root of the (positive) solution, and C − is the region between the -axis and C 0 ; Pr(C − ) coincides with the -axis.In the bottom panels of Figure 4, the green line shows C 0 for three representative values of .It is also apparent how B − ⊆ D − , B − ∩C + = B − and B + ∩C − = C − , as expected from the general discussion; similarly to the MN model, notice that the B + ∩ C + region becomes larger as  decreases, i.e. for more flattened systems.Again, in the spherical case ( = 1) we would obtain We now determine the constraints on () from the request of positivity for  2  , Δ  and Δ.From equation ( 19), with 1/ while, for  ≤  M ,  M () is a complicated function monotonically decreasing from 2/ 2 to the minimum reached on the equatorial plane; in the spherical case ( = 1),  M () = 2 for all .From equation ( 21), and, for  ≤  0 ,  0 () is a complicated function monotonically decreasing from 1/ 2 to the minimum reached again on the equatorial plane; in the spherical case,  0 () = 1 for all , as expected from the general discussion.Finally, from equations ( 23) and ( 24),  1 () = 1, and while, for  ≤  2 ,  2 () is again a complicated function which decreases monotonically, reaching its minimum on the equatorial plane; in the spherical case,  2 () = 1.In the bottom-left panel of Figure 5, we show the four () functions described above for the representative case of  = 0.75: the chain of inequalities in equation ( 25) is apparent, even though the profiles qualitatively differ from those of the two disk cases.In particular, if one restricts to the constant  case, we have Figure 7 shows the 2D maps of  2  , Δ  and Δ for the Binney logarithmic halo with  = 0.75 and for three values of , i.e.  = 0.5, 1, 3.75.The adopted values of  allow to illustrate three representative kinematical behaviours.When  = 0.5, we expect  2  and Δ  to be everywhere positive, while a negative Δ is expected over some region in C − , being 0.5 <  1 = 1, and positive over B − , being 0.5 <  2 ≃ 2.16.This is nicely confirmed by the three left panels, where white regions indicate negative values.The three central panels correspond to  = 1, i.e. to the two-integral solutions, and complement the three middle panels in Figure 3.The fields for  = 1 case are everywhere positive, in agreement with the general discussion and equations ( 49), ( 51) and ( 53).The increase of  from 0.5 to 1 produces, as expected, that  2  decreases over B − and increases over B + , Δ  decreases over D − and increases over D + , and Δ behaves qualitatively as  2  , becoming non-negative over C − because 1 =  1 .Finally, the three right panels correspond to the unphysical model with  >  M ≃ 3.33, and a portion of B − is already white; compared to the disk models, this region is more confined near the equatorial plane.Increasing further  would extend the white region to cover the whole B − , while  2  would continue increasing in B + .As 3.75 >  0 ≃ 1.26, Δ  is now negative over a large portion of D − and, increasing further , Δ  would become negative over the whole D − region, i.e. below the yellow line.At variance with disks, however, Δ  would remain positive and increase above the line.Finally, the negative region of Δ has now switched to the B − region, being 3.75 >  2 .

A TWO-COMPONENT MODEL: THE MIYAMOTO-NAGAI DISK IN A BINNEY LOGARITHMIC HALO
We now study whether and how the presence of a dark matter (DM) halo changes the conclusions obtained for one-component models.In particular, we consider the two-component model made by a stellar MN disk embedded in a DM Binney logarithmic halo.We adopt this model because  2  and the commutator for the stellar distribution can be given in closed form for arbitrary values of  and  (Smet et al. 2015, hereafter S15); the resulting expressions are however sufficiently complicated to exclude the possibility of a simple analytical study of the B, C and D regions, which are determined in a numerical way starting from the analytical formulae in S15.
By using the normalizations adopted for the MN model, i.e.  * for lengths and   * / * for potentials, the dimensionless total potential is in which Φ * = −1/, as given in equation ( 26), and from equation ( 42) the dimensionless halo potential is where  h is the scale length of the halo now in units of  * , i.e. a dimensionless parameter.For R = 0, the model reduces to the self-gravitating MN disk, while for R ≫ 1 we obtain the "halodominated" case, in practice equivalent to considering Φ = Φ h in the JEs.Notice that, at sufficiently large  = √︁  2 +  2 , the model is always halo-dominated, independently of the value of R.
In the two-integral case,  2  of the stellar distribution can be obtained by adding equations ( 10) and ( 17) of S15, while [ * , Φ] can be obtained from equation ( 27) of S15.It can be shown that, while  2  ≥ 0 regardless of the model parameters, the sign of [ * , Φ] depends on the value of , and if  < 1, there are regions at large  where [ * , Φ] < 0. Instead, if  ≥ 1, then [ * , Φ] ≥ 0 everywhere.We notice that an asymptotic analysis of the solutions of the JEs shows that, independently of  and R,  2 becomes negative at large distances from the origin for  ≲ 0.85; instead, for  ≳ 0.85,  2  ≥ 0 everywhere.
As we expect that the halo-dominated models are the most different compared to the one-component MN disk, in the following we focus on these models, and for simplicity, we consider a spherical ( = 1) halo.In the three bottom panels of Figure 3, we show the maps of  2 , [ * , Φ] and Φ/ for a halo-dominated MN disk with  = 1, and  h = 2.By comparison with the MN panels in the figure, it is apparent that the main effect of the halo is to produce an increase of  2  , and of the commutator at large radii near the equatorial plane, as expected given the flat rotation curve of the halo at large radii.Instead, in the more central regions, the self-gravitating effects of the disk are not visible because the gravitational field of the disk is neglected everywhere in the presented model.
However, even though the gravitational potential is dominated by the halo, the disk density distribution enters the ,  and  functions, so we expect that the B, C and D regions, whose determination now must be carried out numerically, are not coincident with those of the Binney halo (see Section 5.3).Indeed, now the B and C regions are more similar to those of the MN disk, and D − coincides with the (, )-plane, with  = 0 on the -axis.
Due to the numerical determination of the B, C and D regions, also the various constraints on the  values had to be found numerically.In the bottom-right panel of Figure 5, we show the four () constraints for the halo-dominated model with  = 1,  = 1 and  h = 2. Reassuringly, they obey the inequality in equation ( 25), being the commutator positive everywhere.For the constant  case, they can be obtained analytically through an asymptotic analysis, and we found that  1 =  0 =  2 = 1, and  M = 5/4.
Figure 8 shows the 2D maps of  2  , Δ  and Δ for the same  = 1,  = 1 and  h = 2, and for  = 0.5, 1, 2 (as in Figure 6 for the MN disk).The main result is that the fields are more similar to those of the one-component MN disk in Figure 6 than to those of the onecomponent Binney model in Figure 7, even though the MN disk is halo-dominated.Therefore, we conclude that under the -ansatz a DM halo does not alter appreciably the qualitative behaviour of the kinematical fields of the disk, even when the halo is very massive.

SUMMARY AND CONCLUSIONS
Unless a specific functional form of the phase-space DF is assumed, the JEs do not form a closed system of equations, and to solve them a choice must be made for the link between the velocity moments.In this work, we considered the JEs of axisymmetric systems with three unknowns ( 2  ,  2  and  2  ), and focused on the choice of  2  = (, ) 2  (the "b-ansatz", where  is a positive function), introduced in Cappellari ( 2008) and widely used in the dynamical modelling of stellar systems.After recalling the JEs in cylindrical coordinates, and the expressions for the ansatz-independent quantities  2 and [ * , Φ], and the associated functions ,  and  defined over the (, )plane, we investigated the behaviour of the kinematical quantities that depend on  as follows: • We first gave the general expressions, as a function of  and of the ansatz-independent quantities, for  2  and the related functions Δ =  2  −  2  and Δ  =  2  −  2  , which enter the Satoh (1980) decomposition of  2  to derive the ordered rotational velocity   .These general expressions are fundamental to investigate the positivity of  2  , required for the model consistency, and of Δ and Δ  , necessary for the Satoh (1980) decomposition and its variants.These general expressions also allow to predict, before solving the JEs, the trends of  2  , Δ and Δ  as a function of  throughout the galaxy, thus avoiding a time-consuming numerical exploration of the parameter space in the modelling.
• In particular, for systems with [ * , Φ] ≥ 0, we showed that as  increases, in addition to the trivial fact that  2  =  2  increases,  2  and Δ decrease over B − and increase over B + , and Δ  decreases over D − and increases over D + .
After the general analysis of the problem and the setup of the procedure to investigate the kinematical properties of a model as a function of , we illustrated the method with three widely used galaxy models: the Miyamoto & Nagai and Satoh disks, and the Binney logarithmic halo.In doing so, we obtained the following results: • The shape of the B ± , C ± and D ± regions and the limits on () can all be obtained analytically.For the two disks, D − (the region where   *  2  / < 0) coincides with the (, )-plane, independently of the disk flattening; for the Binney halo, instead, there is a critical value of the potential flattening such that for flatter potentials the region D + appears, bounded by a curve containing the -axis.For all the models, B − and C − are also bounded by hyperbola-like curves but contain the -axis.Concerning the limits on (), for all the models Pr(B − ) = Pr(C − ) = Pr(D − ) coincide with the -axis, and  1 () ≤ 1 ≤  0 () ≤  2 () ≤  M ().
• With the aid of two-dimensional maps, we illustrated how the kinematical fields change over the galaxy by changing , assumed spatially constant.The results are in agreement with the general considerations, both for their trend with , and for the regions where they become negative, confirming the utility of the preliminary analysis.
• We finally applied the procedure to the two-component model made by a stellar MN disk embedded in a dominant, spherical dark matter Binney halo; in this case, the investigation has been carried out numerically.The dark halo does not alter appreciably the qualitative trend of the kinematical fields of the stellar disk as a function of , even when the halo is very massive.
A final comment is worth on the relation between  2  and  2  as a function of  and , the other parameter that measures the orbital anisotropy, with positive values corresponding to radial anisotropy, and negative ones to tangential anisotropy [see equations ( 16) and ( 17)].As noticed in Section 3.1,  is necessarily positive where Δ  < 0; where Δ  > 0, instead,  can be positive ( 2 > 1 in the Satoh decomposition) or negative ( 2 < 1).Allowing for  and  as free parameters in the JAM method, the stellar kinematics of regular rotators in the ATLAS 3D survey could be successfully modelled (e.g., Cappellari 2016).It turned out that on average  > 1 and  ≈ 0, within about one effective radius; a limit was also found in the orbital anisotropy, of  = 1 − 1/ < ∼ 0.7 intr , where  intr is the intrinsic flattening of the galaxy.This limit had first been found with the Schwarzschild orbit superposition method for the galaxies of the SAURON survey (Cappellari et al. 2007), and thereafter has been often shown in the (/, ), and (  , ) diagrams of larger samples of regular rotators, until those of the MaNGA survey (e.g., Wang et al. 2020).In a recent work, Wang et al. (2021) investigated the origin of this limit: with the JAM modeling, they built mock samples of axisymmetric galaxies with a velocity dispersion ellipsoid appropriate to reproduce the observed fast rotators ( ≈ 0), and with the mass-follows-light hypothesis.They found an upper limit on  by excluding values producing a negative   2 over a non-negligible region of the galaxy, clearly an unphysical solution.Our analysis can provide an interpretation for the existence of this limit.In fact, from equation ( 16), the request of positivity for   2 with  ≈ 0 translates into the request of Δ  ≳ 0, and then of  ≲  0 .Interestingly, for the mass models considered in this work,  0 increases with the flattening [see equations ( 31), (39), ( 50) and ( 51)], in agreement with the trend shown by the maximum values estimated for  (of  max ∼ 0.7 intr ).
Each region is made by two disjoint parts fully covering the (, )-plane, thus the determination of just B − , C − and D − suffices for investigating a given model.In Appendix B, some general results about the relative positions of the regions and the constraints on () are presented.In particular, it is shown that B − is always contained in D − .Moreover, if [ * , Φ] > 0 over the whole (, )-plane (as for oblate self-gravitating models), then B − ⊆ C + and C − ⊆ B + ; if [ * , Φ] = 0 everywhere (as in spherically symmetric models), then B − = C + , B + = C − , and D −0 coincides with the whole (, )plane; if [ * , Φ] < 0 (as for prolate self-gravitating models), then B + ⊆ C − , C + ⊆ B − , and D − coincides with the whole (, )plane.Finally, from Theorems B.5 and B.7, it follows that for models with [ * , Φ] ≥ 0 everywhere (a quite common situation), one has:

Figure 2 .
Figure 2. Scheme representing the behaviour, as a function of , of the kinematical fields  2  , Δ  and Δ over the regions B ± , C ± and D ± , for a model with [, Φ] > 0. Each row refers to the indicated region; arrows point upwards if the field is increasing, and downwards if the field is decreasing for increasing ; black arrows indicate a positive value, and red arrows a negative one.For  >  M , the model is unphysical.

Figure 3 .
Figure 3. Maps in the (, )-plane of  2  (left column),  [, Φ] / (central column), and Φ/ (right column) for the MN disk with  = 1 (top row), for the Binney logarithmic halo with  = 0.75 (central row), and for the MN stellar disk ( = 1) in a dominant (Φ = Φ h ) and spherical ( = 1) Binney logarithmic dark matter halo with  h = 2.The fields are in units of  2 h for the Binney halo, and of   * / * for the other two models.Black solid lines show equally spaced isodensity contours.

Figure 4 .
Figure 4.The ansatz-independent regions B, C , and D of the MN disk (top panels) for different values of , and of the Binney logarithmic halo (bottom panels) for different values of ; both models become flatter from left to right.Red and green areas correspond to the B − and C − regions, respectively, while heavy red and green lines are the B 0 and C 0 sets.Except for the Binney model with  = 0.75, where D − lies below the D 0 line (dashed), for all the other cases D − fully covers the (, )-plane, with D 0 coinciding with the -axis.For all models, Pr(B − ) and Pr(C − ) coincide with the -axis, and as [, Φ] ≥ 0 everywhere, from Theorem B.2 we have B − ∩ C + = B − , and C − ∩ B + = C − . )].

Figure 6 .
Figure 6.Maps in the (, )-plane of  2  (top row), Δ  (central row) and Δ (bottom row) in units of   * / * for the MN disk with  = 1 and constant  = 0.5 (left column),  = 1 (central column) and  = 2 (right column).Black solid lines show equally spaced isodensity contours of the stellar distribution.The red and green lines show, respectively, the B 0 and C 0 curves, while D 0 coincides with the -axis.White regions correspond to negative values of the fields.

Figure 7 .
Figure 7. Maps in the (, )-plane of  2  (top row), Δ  (central row) and Δ (bottom row) in units of  2 h , for the Binney logarithmic halo with  = 0.75 and constant  = 0.5 (left column),  = 1 (central column) and  = 2 (right column).Black solid lines show equally spaced isodensity contours of the stellar distribution.The red, green and black dashed lines show, respectively, the B 0 , C 0 and D 0 curves.White regions correspond to negative values of the fields.

Figure 8 .
Figure 8. Maps in the (, )-plane of  2  (top row), Δ  (central row) and Δ (bottom row) in units of   * / * for the MN disk with  = 1 in a dominant (Φ = Φ h ), spherical ( = 1) Binney logarithmic halo with  h = 2, and constant  = 0.5 (left column),  = 1 (central column) and  = 2 (right column).Black solid lines show equally spaced isodensity contours of the stellar distribution.The red and green lines show, respectively, the B 0 and C 0 curves.White regions correspond to negative values of the fields.