Rotational equilibria by Lagrangian variational principle: toward multi-dimensional stellar evolutions
aa r X i v : . [ a s t r o - ph . H E ] S e p Mon. Not. R. Astron. Soc. , 000–000 (0000) Printed 26 September 2018 (MN L A TEX style file v2.2)
Rotational equilibria by Lagrangian variational principle:toward multi-dimensional stellar evolutions
Nobutoshi Yasutake , Kotaro Fujisawa , Shoichi Yamada Physics Department, Chiba Institute of Technology, Shibazono 2-1-1, Narashino, Chiba, 275-0023, Japan Advanced Research Institute for Science and Engineering, Waseda University, Okubo 3-4-1, Shinjuku, Tokyo, 169-8555, Japan
26 September 2018
ABSTRACT
We have developed a new formulation to obtain self-gravitating, axisymmetric configu-rations in permanent rotation. The formulation is based on the Lagrangian variationalprinciple with a triangulated mesh. It treats not only barotropic but also baroclinicequations of state. We compare the various stellar equilibria obtained by our newscheme with those by Hachisu’s self-consistent field scheme for the barotropic case,and those by Fujisawa’s self-consistent field scheme for the baroclinic case. Includedin these rotational configurations are those with shellular-type rotations, which arecommonly assumed in the evolution calculation of rotating stars. Although radiationprocesses, convections and meridional flows have not been taken into account in thisstudy, we have in mind the application of this method to the two-dimensional evolutioncalculations of rotating stars, for which the Lagrangian formulation is best suited.
Key words: stars: rotation— stars: evolution — stars: protostars
Stellar evolution theory is well established especially forspherically symmetric stars, in which there is no rota-tion and no magnetic field. After the works of pioneersof the field (Burbidge et al. 1957; Cameron 1957; Hayashi1961ab), Henyey proposed in his seminal papers a method,which later became the defacto standard for stellar evolu-tion calculations (Henyey et al. 1964). The method is stableand capable of self-consistently incorporating various phys-ical processes that occur in the stellar interior. It has beenmodified and extended a lot approximately in the last halfcentury to accommodate rotation(Maeder & Meynet 2000;Woosley et al. 2002).In stellar evolution calculations, we need to obtain stel-lar structures consistent with, nuclear reactions (and/ormolecule formations), radiative transport of energy, con-vective and circular motions as well as a realistic equa-tion of state. Since the time-scale of stellar evolution isnormally much longer than the dynamical time-scale, it iswell approximated by series of time-independent configura-tions in hydrostatic equilibrium. In the presence of rota-tion, this is not simple. One of the difficulties is to obtainrotational equilibria for a given angular momentum distri-bution. What is more, however, we do not know a prioriwhat the distribution of angular momentum looks like inthe stellar interior although the problem has been studiedtheoretically since the 19th century by Carl Jacobi, RichardDedekind, Peter Lejeune Dirichlet, and Bernhard Riemann to mention a few (Chandrasekhar 1969). It was in the lastcentury that some methods to obtain rotating hydrostaticequilibria were proposed, and it was mathematically shownthat cylindrical distributions of specific angular momen-tum are realized for EOSs, in which matter pressure is afunction of density alone. This type of EOS is referred to“ barotrope ” (Ostriker & Mark 1968; Hachisu 1986).The EOS is not barotropic in the realistic stellar in-terior, but baroclinic in general, i.e., pressure dependson density, entropy and chemical compositions . Thenthe numerical construction of rotational equilibria becomesmuch more difficult compared with the case for barotropes,since the first integral of the hydrostatic equations thatis available in the barotropic case no longer exists in thebaroclinic case (Uryu & Eriguchi 1994; Uryu & Eriguchi1995; Roxburgh 2006; Espinosa Lara & Rieutord 2007;Espinosa Lara & Rieutord 2013; Rieutord et al. 2016). Evenif some rotational configurations are somehow obtained, itis a highly non-trivial issue how to make a sequence out ofthem that represents stellar evolutions appropriately. Notethat in the presence of rotation fluid elements composingrotating stars could change their positions non-radially incomplicated ways as the stars evolve in time even withoutconvection. It would be highly difficult to describe such dis-placements of fluid elements with fluxes on a fixed numerical For compact stars such as white dwarfs and neutron stars, theEOS may be approximately barotropic , since the temperature isnegligibly low.c (cid:13)
N.Yasutake, K. Fujisawa, and S. Yamada mesh from a Eulerian point of view, since the motions areextremely slow. The Lagrangian formulation will be hencemore appropriate just as in the spherically symmetric case.And this is exactly the idea in this paper. We introduced inthis paper a triangulated mesh with each node representinga fluid element approximately. Starting from an arbitraryreference configuration, we seek for the Lagrangian displace-ment for each node that gives a rotational equilibrium as aresult. In so doing, the variational principle is employed. Themethod has its own difficulty, though.In this paper, we construct a couple of configurationsin rotational equilibrium with different angular momentumdistributions in order to demonstrate the capability of ournew scheme. Included in them is the so-called shellular-typerotations. It has been argued that in the stellar interiorturbulence is anisotropic, operating more strongly withineach isobaric surface than in its normal direction, and thata uniform rotation should be realized in the isobar as aconsequence. Such a rotation law is referred to as “shellu-lar rotation” rather than cylindrical rotation (Zahn 1992;Meynet & Meynet 1997). It should be stressed, however,that no rotational equilibrium with the shellular rotationwas constructed numerically, not to mention analytically .Very recently Fujisawa, one of the authors of this paper, hassucceeded in producing some configurations with shellular-type rotation(Fujisawa 2015). It is a nice demonstration ofperformance of our new scheme to reconstruct these config-urations.The organization of the paper is as follows. In Section2, we describe the formulation in detail: the Lagrangianvariational principle, on which our scheme (referred to asthe YFY scheme hereafter) is based, will be reviewed first;then the finite discretization on the triangulated grid willbe explained and the handling of self-gravity is also men-tioned; finally the minimization of the energy function withthe Monte Carlo technique, which is analogous to those uti-lized, e.g., in nuclear physics, will be described. In Section3, we demonstrate the nice performance of YFY scheme,constructing some rotational configurations and comparingthem with the those obtained by other numerical schemes,including Fujisawa’s self-consistent field scheme mentionedabove. In the last section, we summarize this study and dis-cuss possible extensions of our scheme as future works.
Rotational equilibria should satisfy the following equations:(i) the force balance equation: ∇ P + ρ ∇ φ − ρj ̟ e ̟ = 0 , (1)where P , ρ , φ , j , ̟ , e ̟ are pressure, density, gravitationalpotential, specific angular momentum, the distance from therotation axis, and the unit vector perpendicular to the ro-tation axis, respectively; The obvious exception is uniformly rotating stars. By ”shellular-type” we mean that the configuration is not per-fectly but almost shellular with iso-angular-velocity surfaces beingnearly aligned with isobars. (ii) the Poisson equation:∆ φ = 4 πGρ, (2)where G is the gravitational constant; these equations aresupplemented with the equations of state (EOS).Although equations (1), (2) look simple, especially for avery realistic EOS, numerically solving them may not be soeasy. Most of previous works assumed a barotropic EOS, inwhich pressure is a function of density alone, and for which itis mathematically shown that rotation is uniform on concen-tric cylinders and, as a consequence, there is a first integral ofequation (1). The realistic equation of state, however, is notbarotropic but baroclinic, in which case pressure dependson density, entropy and chemical compositions. Then thereis no first integral available any longer and the rotation lawbecomes unknown a priori. In this paper, we do not use theforce balance equation, equation (1). Instead the variationalprinciple is utilized as described in detail below. We first describe the Lagrangian variational principle, whichis obeyed by rotational equilibria. We obtain the energyfunctional, which is minimized by a rotational equilibriumand will be discretized on a triangulated mesh later fornumerical calculations. The variational principle itself isnot our own invention but is actually given in a textbook(Tassoul 1978), although it is a core of the entire formulationproposed in this paper.Suppose an arbitrary reference configuration with den-sity, specific entropy and specific angular momentum dis-tributions ( ρ ( ¯ r ) , s ( ¯ r ) , j ( ¯ r )), where ¯ r is a spatial coordinateand will also serve as a Lagrange coordinate as will becomeclear shortly. We, then, introduce a diffeomorphic mapping: r = f ( ¯ r ), which is identified with the Lagrangian displace-ment that deforms the reference configuration to anotherconfiguration. Then the following relation should hold if themass were to be conserved in the deformation: ρ ( ¯ r ) dV ( ¯ r ) = ρ ( r ) dV ( r ) , (3)where ρ ( ¯ r ) is a density at ¯ r and dV ( ¯ r ) is an infinitesimalvolume around ¯ r in the reference configuration whereas ρ ( r )and dV ( r ) are the density and volume at the correspondingposition in the mapped configuration.The volume element istransformed with the Jacobian J of the mapping as dV ( ¯ r ) = J ( f − ) dV ( r ) , (4)the density after the displacement is expressed as ρ ( r ) = ρ ( ¯ r ) J ( f − ) . (5)Introducing another diffeomorphism g , we consider theconsecutive mapping g · f , which we also refer to as theproduct of maps g and f and interpret as a Lagrangiandisplacement: ¯ r f → r g → r ′ , (6)where r ′ is the coordinate after the mapping g .Using the formula for the product of Jacobians, J (( g · f ) − ) = J ( f − ) · J ( g − ) , (7) c (cid:13) , 000–000 otational equilibria by Lagrangian variational principle we obtain the density after the mapping shown in equation(6) as ρ ( r ′ ) = ρ ( ¯ r ) J (( g · f ) − ) = ρ ( ¯ r ) J ( f − ) J ( g − ) = ρ ( r ) J ( g − ) . (8)This equation implies that the intermediate configurationobtained after the mapping f served as the reference con-figuration for the configuration after the mapping g andone can forget about the original reference configuration.Employing this fact, we consider infinitesimal displacementsalone in the following: δ f : r → r ′ = r + ξ (9)where ξ is the so-called Lagrangian displacement vector, aninfinitesimal vector that connects the original point with itsimage by the map. The Jacobian for the (inverse) infinites-imal mapping is given as J ( δ f − ) = J − ( δ f ) = 1 − div( ξ ) . (10)Then the Lagrangian variation of density ∆ ρ is obtained as ∆ ρ ≡ ρ ( r ′ ) − ρ ( r ) = ( J ( δ f − ) − ρ ( r ) = − div( ξ ) ρ ( r ) . (11)In this study, we consider the variational principle, inwhich rotational equilibria are obtained. Supposing thatthere is no Lagrangian variation in specific entropy and an-gular momentum for the above infinitesimal displacement,∆ s = 0 , (12)∆ j = 0 , (13)we consider the variation of the following functional, towhich we refer as the energy functional hereafter: E ( ξ ) = Z ερ dV + 12 Z ρφ dV + Z ρ (cid:16) j̟ (cid:17) dV. (14)In this expression the integration domain is the stellar in-terior, and ε and φ denote the specific internal energy andgravitational potential, respectively; ̟ is the distance fromthe rotational axis. We consider the variation of each termon the right hand side of equation (14) in turn. The variationof the first term is evaluated as follows: δ (cid:18)Z ερ dV (cid:19) = Z ∆ ε ρ dV = Z (cid:18) ∂ε∂ρ (cid:19) s ∆ ρ ρ dV = Z (cid:18) Pρ (cid:19) ∆ ρ ρ dV = − Z P div( ξ ) dV = Z ∇ P · ξ dV. (15)Here ∆ ε is the Lagrangian variation of internal energy, whichis given by the Lagrangian variation of density upon usingthe assumption that there is no entropy variation, equation(12). The final expression is obtained with the use of equa-tion (11).The second term is evaluated as follows: δ (cid:18) Z ρφ dV (cid:19) = 12 Z ( δρ · φ + ρ · δφ ) dV = 12 Z ( −∇ ( ρ ξ ) · φ + ρ · δφ ) dV This is regarded as a functional of the map f . Its variation is afunctional of δ f and hence of ξ as shown shortly in the following. = 12 Z ( ρ ξ · ∇ φ + ρ · δφ ) dV, (16)where the second equality is obtained with an employmentof the following relation: δρ = ∆ ρ − ξ ∇ ρ. (17) δφ in the last term of the final expression is the Eulerianvariation of the gravitational potential, which is evaluatedfrom the Poisson equation, ∆ φ = 4 πGρ, and its solutionwith the Green’s function, φ = − Z dV ′ ρ ( r ′ ) / | r − r ′ | . (18)Then R ρ · δφ dV is given as Z ρ · δφ dV = − G Z δρ ( r ′ ) ρ ( r ) | r − r ′ | dV ′ dV = Z δρ ( r ′ ) φ ( r ′ ) dV ′ = Z ρ ξ · ∇ φ dV. (19)Equation (16) is hence expressed as δ (cid:18) Z ρφ dV (cid:19) = Z ρ ξ · ∇ φ dV. (20)The variation in the last term of equation (14) is calcu-lated as δ (cid:18)Z ρ (cid:16) j̟ (cid:17) dV (cid:19) = Z ρ ∆ (cid:18) j ̟ (cid:19) dV = Z j ξ · ∂∂ r (cid:16) ̟ (cid:17) ρ dV = − Z j ̟ ξ · e ̟ ρ dV. (21)In the second equality we employed the assumption thatthere is no Lagrangian variation in the specific angular mo-mentum, equation (13).We are now ready to derive the Euler-Lagrange equa-tion for the energy functional in equation (14). All the vari-ations given by equations (15), (20) and (21) are expressedwith the Lagrangian displacement vector ξ and the Euler-Lagrange equation corresponds to its coefficient in δE . Theresultant equation is obviously equation (1). This impliesthat the rotational equilibrium is a stationary point of thefunctional and that it may be obtained not by solving equa-tion (1) but by minimizing the energy functional somehow .In the above formulation of the variational principle, wetreated the gravitational potential φ as a functional of den-sity, formally solving the Poisson equation with the Green’sfunction. It is possible, however, to treat the gravitationalpotential also as an independent variational variable. Theenergy functional should be modified then to E ( ξ ) = Z ερ dV + Z ρφ dV + Z ρ (cid:16) j̟ (cid:17) dV + 14 πG Z
12 ( ∇ φ ) dV. (22) The rotational equilibrium may not be a minimum point of theenergy functional, which will be indeed the case if the configu-ration is convectively unstable. Possible treatments of such caseswill be discussed elsewhere.c (cid:13) , 000–000
N.Yasutake, K. Fujisawa, and S. Yamada
In this case the Poisson equation for the gravitational po-tential itself is obtained as a Euler-Lagrange equation withrespect to φ . Note that the factor of 1/2 disappears fromthe second term, reflecting the fact that φ is an independentvariable. The problem with this formula from a practicalpoint of view is that the fourth term has a non-compact sup-port, i.e., the integral region extends to infinity. This prob-lem may be alleviated, however, by imposing an appropriateboundary condition at the stellar surface . In the following,we employ the first formulation, in which the gravitationalpotential is treated as a functional of density through thesolution of the Poisson equation. The second form of theenergy functional, equation (22), is employed, however, toderive a discretized Poisson equation in section 2.3. The Lagrangian variational principle presented above is notour own invention but has been known for many years(Tassoul 1978). The following implementation, however, isour original idea. We adopt a finite element method: themeridian section of an axisymmetric star is covered by a tri-angulated mesh, and approximate the energy functional onthis mesh. Various quantities are assigned to each grid pointor node. They are moved to artificially deform the star. Weexplain the method of these procedures in this subsection.The triangulated mesh consists of nodes and edges, withadjacent nodes being connected by an edge. The so-calledadjacency matrix, which is commonly used in graph theoryand gives the relation between neighboring nodes with its( i , j ) component A ij being unity (zero) if the i -th and j -thnodes are (not) connected with an edge, is convenient andwill be used in evaluating various quantities numerically.Note that this matrix is symmetric. To each node we assigncoordinates, mass, specific angular momentum, entropy andvolume fractions of various elements . The latter four areconserved quantities as long as we do not consider nuclearreactions and transport of energy and angular momentumand are fixed in the node shifts, ρ i = m i /V i = m i / ( 13 X cell ∈ i V ∆ ) . (23)In this expression, m i and ρ i are the mass and density as-signed to the i -th node, respectively, and the summationruns over the cells surrounding the i -th node whose volumesare denoted by V ∆ .The energy functional, equation (14) is approximatedon this mesh as: E FEM ( r i ) = X i ε i m i + 12 X i φ i m i + X i (cid:16) j i ̟ i (cid:17) m i , (24)which is now actually a function of the coordinates of allnodes, r i , with the subscript specifying the individual node.This corresponds to an approximate evaluation of the Jaco-bian J ( r i ) at each node. The minimization of this energy This is nothing but the evaluation somehow of the integral out-side the star. See section 2.4. In this paper we ignore the fractions of elements for simplicity. gives the coordinates of nodes in the rotational equilibriumfor given distributions of mass, specific entropy and angularmomentum on the triangulated mesh. We then get the resul-tant density at each node through equation (23). It shouldbe obvious that the specific internal energy ε i is also re-garded as a function of the coordinates of nodes, since it isa function of density and specific entropy s i , the former ofwhich is obtained once the node positions are determinedas just mentioned and the latter is fixed. As for the gravita-tional potential φ i , equation (2) is solved numerically on thesame triangulated mesh for the density obtained this way. Itcan be hence regarded also as a function of the coordinatesof nodes. As mentioned just now, in our formulation, the gravitationalpotential is regarded as a functional, or a function in the dis-cretized version, of density by solving the Poisson equation.As also noted earlier in section 2.1, we employ the energyfunctional given by equation (22) to derive the discretizedversion of Poisson equation.We first approximate the second and fourth terms onthe right hand side of the equation on the triangulated meshas X k Z cell (cid:16) ( 14 πG ) 12 ∇ φ · ∇ φ + φρ (cid:17) dV, (25)where the volume integral is done for each triangular celland summed over all cells. The integrand is approximatedlinearly in each cell as φ = α φ z + α φ ̟ + α φ , (26)where the coefficients, α ’s, are given in terms of the valuesof the gravitational potential at three nodes of the cell as α φa = T − ab φ b . In this expression, the subscripts a and b specify the node ofthe cell and the matrix T ab is determined by the geometry ofthe mesh alone and given in Appendix A. Then the deriva-tive of φ in the integrand of equation (25) is approximatedas ∇ φ = ( α φ , α φ , . (27)Using this approximation, we evaluate the first term of equa-tion (25) as X k Z cell ∇ φ · ∇ φ dV = X k (cid:0) α φ α φ + α φ α φ (cid:1) Z cell dV = X k X a,b ∈ k (cid:0) T − a T − b + T − a T − b (cid:1) φ a φ b V ∆ . (28)The second term of equation (25) is, on the other hand,expressed as X k Z cell φρ dV = X i φ i m i , (29)in which the sum on the right hand side is taken over all c (cid:13) , 000–000 otational equilibria by Lagrangian variational principle nodes. Collecting the terms, we obtain the expression forequation (25) as X k X a,b ∈ k (cid:0) T − a T − b + T − a T − b (cid:1) V ∆ · φ a φ b = X i φ i m i . (30)The derivative of this expression with respective to φ i givesthe discretized Poisson equation as X j M ij φ j = 4 πGm i . Here we introduce a symmetric matrix M ij defined as X i,j M ij φ i φ j = X k X a,b ∈ k (cid:0) T − a T − b + T − a T − b (cid:1) V ∆ φ a φ b , (31)in which the double sum on the left hand side runs overall nodes. Note that the matrix M ij is determined by themesh geometry alone. This approximate Poisson equation issolved algebraically as φ i = 4 πG M − ij m j . (32)In so doing the boundary condition is imposed Dirichlet-typeat the stellar surface with the values of the gravitationalpotential being calculated with the Green’s function (seebelow).In order to validate the formulation, we compare it withanother method. In fact, the Green’s function method, whichis employed to calculate the boundary values of the poten-tial, is actually applicable everywhere. Adopting the multi-pole expansion of the Green’s function, we write the poten-tial as φ i = − πG ∞ X l =0 P l (cos θ i ) Z π sin θ ′ dθ ′ P l (cos θ ′ ) × (cid:18)Z r i ( r ′ ) l +2 r l +1 dr ′ + Z ∞ r i r ′ l ( r ′ ) l − dr ′ (cid:19) ρ i ( r ′ , θ ′ ) , (33)where P l ’s are the Legendre polynomials. The integrals areevaluated numerically.The comparison is done for a non-spherical density dis-tribution shown in the left panel of Fig. 1. The triangulatedmesh we deploy is displayed with the nodes and edges in theright panel. The mass m i associated with node i is obtainedfrom the density ρ i at the node through equation (23). Weevaluate the gravitational potential either via equation (32)or equation (33). The left and middle panels of Fig. 2 corre-spond to the numerical results for equations (32) and (33),respectively, which look almost identical to each other. As amatter of fact, the right panel of Fig. 2 demonstrates thatthe both results are in agreement within less than 1% errors.Here the error is defined as the absolute value of the ratio ofthe difference between the two results to the one obtainedwith equation (33). We employ a Monte Carlo technique to search for minimumenergy function. The reason is the following. It is well knownfrom the perturbation theory based on the Lagrangian dis-placement vector that a class of trivial
Lagrangian dis- placements make the distributions of physical quantities un-changed (Friedman & Schutz 1978). An example is the isen-tropic case, i.e. the specific entropy is uniform in the stellarinterior. This is most easily understood from the fact thatthe Lagrangian displacement vector ξ = 1 /ρ ∇ × η for anarbitrary axisymmetric vector field η , which would be ob-tained as a Euler-Lagrange equations for the discretized en-ergy function, equation (24), does not change the density asa function of spatial position. This residual gauge degree offreedom renders the system of equations for r i indefinite andwould make singular the matrix obtained by linearizationif one were to employ the Newton-Raphson method. Thismeans that, without an appropriate gauge fixing, the for-mulation could not be applied to isentropic configurations,which are normally supposed to be the simplest case. Even ifthe specific entropy is not homogeneous, the energy functionis rather insensitive to the displacements of nodes close tothe rotation axis and/or to the surface in general, which im-plies that there are always small (i.e. almost singular) eigen-values in the matrix of the linearized equations, which couldhamper the employment of the Newton-Raphson method.The Monte Carlo method described below could natu-rally alleviate these difficulties. In this technique we displacerandomly the grid points and search literally for the con-figuration that makes the energy function minimum. Theproblem of the gauge degree of freedom does not arise inthe Monte Carlo method, since one solution is automati-cally picked up out of infinitely many possibilities by a dice.And that is sufficient indeed for our purpose because physi-cal quantities are identical for any choice. Such an approachis not special and is indeed employed in other fields such asnuclear physics, in which deformed nuclei are constructed ina way that is analogous to the one we adopt in this paper.The actual minimization procedure is the following.Given an initial configuration, we move one of the nodesin the radial direction slightly and see if this shift lowersthe value of the energy function or not. If not, we cancelthe shift. The amount of dislocation is a random variablebut is limited to ∼
1% of the lengths of the edges attachedto the node. We do the same things in the lateral directionand sweep the entire mesh. We repeat this process until thevalue of the energy function is no longer lowered comparedwith the one of 100 sweeps before. This condition is, how-ever, not enough to get hydro-static equilibria, and we alsoneed to check virial residuals, which will be described later.In the current version of the formulation, the connectionsof nodes by edges represented by the adjacency matrix A ij introduced in subsection 2.2 are fixed through the minimiza-tion process. It is hence important to perform the initial tri-angulation properly so that the iteration process would leadto convergence. This is particularly the case when rotation isvery rapid and the resultant rotational equilibrium is differ-ent and some triangular cells may become highly deformed.It may be necessary then to implement the re-gridding some-how.The nodes on the outer boundary, which are referredto as the anchor nodes and have a vanishing mass, are notmoved during this minimization process, since an expansionin the radial direction would continue for ever otherwise. Insome cases, however, the outer boundary comes too close toor goes too far from the nearby active nodes. We then shiftthe anchor nodes appropriately and fix them again there- c (cid:13) , 000–000 N.Yasutake, K. Fujisawa, and S. Yamada ρ
0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] Z [ k m ] X [10 km] Figure 1. (colour on-line). An example of density distribution (left panel), and the edges and nodes (right panel). The colour bar showsthe density in cgs unit. ∆φ
0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] Figure 2. (colour on-line). Left, middle panels show the gravitational potential by the expressions in equation (32) and equation (33)for the density profile given in Fig. 1. The unit is in cgs. Right panel shows the difference ratio between the results of left and middlepanels. The ratio is given by the difference divided by the result of middle one. after. We find that large deformations of triangular cells de-grade accuracy of the energy function and in some casesgenerate a fictitious local minimum. To avoid such artifactsthe re-gridding may be a solution as mentioned above but inthe current version we apply a smoothing to the deformedportion of the grid, which has experienced more than 30 %of change in area by a single Monte Carlo sweep. We knowempirically that such a large variation in area is commonlyaccompanied by the appearance of a very acute angle in thecell. We introduce the residual in the virial relation definedas V C = (cid:12)(cid:12)(cid:12)(cid:12) T + W + 3 R P dVW (cid:12)(cid:12)(cid:12)(cid:12) , (34)in which T and W are the rotational and gravitational ener-gies, respectively, and the third term is the volume integralof pressure. V C is equal to zero in the exact theory and isa measure of numerical accuracy commonly used in the lit-erature (e.g. Eriguchi & Mueller 1985), and we also adoptit together with energy convergence. In this paper, we havetypically achieved V C < − with ∼
500 grid points, whichis comparable previous studies where V C ∼ O (10 − ) wasachieved (Eriguchi & Mueller 1985). It is possible to improveaccuracy by deploying a larger number of nodes. In this section we apply our new method to some test cases:barotropic configurations in subsection 3.1 and baroclinicones in subsection 3.2. We compare the resultant configu-rations with ones obtained by other Eulerian schemes. Insubsection 3.3, we mention some technical issues such as thesmoothing procedure.
We begin with the barotropic case, in which pressure de-pends on density alone, P = P ( ρ ). We consider the isen-tropic configurations as a representative of barotropes. Herewe compare our solutions with those given by Hachisu’sself-consistent field (HSCF) scheme (Hachisu 1986), awell-established numerical method for rapidly rotatingbarotropic configuration.In this section, we consider two types of rotation laws.As repeatedly mentioned, the barotropic configuration ro-tates cylindrically, i.e., the iso-angular velocity surfaces areconcentric cylinders.(i) rigid rotation:Ω ( r ) = Ω , (35) c (cid:13) , 000–000 otational equilibria by Lagrangian variational principle (ii) differential rotation:Ω ( r ) = Ω ( ̟/R e ) + d , (36)where Ω( r ) is an angular velocity at the position of r . Inequation (36) R e is the equatorial radius, Ω is a parameterthat specifies the overall magnitude of angular velocity, and d is a dimensionless parameter that determines the a degreeof differential rotation: as d increases, the angular velocitybecomes uniform whereas in the opposite limit the rotationalvelocity becomes constant (Eriguchi & Mueller 1985). Weset d = 0 .
90 in the following.It should be reminded that in our YFY scheme the dis-tribution of angular velocity is not known a priori, since weassign a conserved specific angular momentum to each nodeof the triangulated mesh in the initial reference configura-tion, which is arbitrary, and the node position is changedlater. This is not the case for the HSCF scheme, which is aEulerian method and derives the density distribution on thefixed mesh, assuming from the beginning that the rotationis cylindrical. Note that even the cylindrical rotation is notguaranteed in the YFY scheme. It is hence a good diagnosisfor our method. In order to make comparison between thetwo schemes, we take the following strategy: we first employthe HSCF scheme to obtain the rotational equilibrium fora given rotation law and cover the resultant configurationwith a triangulated mesh and obtain the mass and specificangular momentum that should be assigned to each node;we then expand by 20 % in the radial direction and use itas the initial reference configuration for our YFY scheme.If it works properly, the mesh should return to the originalshape after the minimization of the energy function. Inci-dentally we have confirmed that if the artificial expansion ofthe mesh is not administered, the mesh does not change bythe minimization as it should.One of such comparisons is shown in Fig. 3, where thedensity distributions obtained with the two schemes are pre-sented for a rigidly rotation. Although HSCF scheme utilizesnon-dimensional quantities are normalized with the max-imum density ρ and equatorial radius R e , we give hererather arbitrarily specific values ρ and R e as ρ = 1 . × g cm − and R e = 2 . × km in order to facilitate com-parison. The ratio of poler radius R p to the equatorial set to R p /R e = 0 .
80 in this model, which corresponds to the ratioof rotational energy to gravitational one of
T / | W | = 3 . s = 14 . B withk B being the Boltzmann constant if the matter is composedof hydrogens alone.The two configurations look quite similar to each other,which is also confirmed in the right panel in which relativedifferences are displayed as a contour color map: the densitydistributions agree with each other within an error of 5 %in most places. A closer inspection reveals, however, thatthe deviation gets larger near the surface, where the densityitself is quite low. The vicinity of the rotation axis is anotherregion, where the accuracy is compromised. This is due tothe fact that our scheme is based on the variational principleand it is intrinsically difficult to determine the configurationof the regions that do not contribute much to the energyfunctional. Some multi-layer treatment may be needed tohandle this problem. In Fig. 4 we show how the energy function (left panel)and the virial residual (right panel) change as the Monte-Carlo sweep proceeds for the same model. The iteration pro-cess is terminated after 248 sweeps, at which point the valueof the energy function is no longer lowered and, importantly,the virial residual is sufficiently small ( < ∼ − ) simulta-neously. One recognizes that there are occasional glitches(three major ones and many minor ones) in the left panel, atwhich the value of the energy function rises discontinuously.There are corresponding discontinuities in the right panel.These irregularities are produced by what we call smoothingwhich is a reconfiguration of the mesh we conduct when themesh gets distorted too much. Without such a measure, thesearch is trapped in local energy minimum as we describelater.The corresponding change of the configuration is pre-sented in Fig. 5. One can see a gradual shrink of the entiregrid. Note that although the energy function looks settled toa minimum around 130 sweeps as shown in the left panel ofFig. 4, the configuration is far from rotational equilibrium.This is an example of the false local minimum just men-tioned. The fact that this is a fake is reflected in the valueof the virial residual given in the right panel of Fig. 3. Thisis the reason why we need to conduct the smoothing and thevirial residual is an important diagnostic quantity to judgeconvergence.Next, we present in Fig. 6 an example of the differentialrotation given by equation (36). As in Fig. 3, the left panelshows the result by the YFY scheme, whereas the middlepanel is for the HSCF scheme. We set again the central den-sity and equatorial radius to ρ = 1 . × g cm − and R e = 2 . × km. The specific entropy is constant in spaceand the same as for the rigid rotation in the previous case.Then the total mass is also unchanged. The density distri-butions agree with each other within 5 % in most regions,but not near the surface again.The behavior of the energy function (left panel) andvirial residual (right panel) is qualitatively the same as whatwe saw in Fig. 4 for the rigid rotation. The iteration processis terminated after 248 sweeps. These panels are very similarto Fig. 7, though the rotation law is different.In Fig. 8, the logs of the stellar structures for a differen-tial rotation given by equation (36) are illustrated again. Thestructure changes from initial model—the results of HSCFscheme expanded by 20 % to the radial direction— to thehydro-static solution. The star continues to transform at 150sweeps, although the energy function is converged in Fig. 7.We also show the distributions of angular momentumin this model as Ω = 1 . × − rad s − , and T / | W | is 4.3%. As for the numerical accuracy, as we see some differences( ∼ not assumed butobtained as a result of computations. This is a strong vindi-cation that our formulation works well indeed. We assumethe other rotational laws including a rigid rotation in thereference configuration, and find in all cases cylindrical ro-tations in the outcomes. c (cid:13) , 000–000 N.Yasutake, K. Fujisawa, and S. Yamada
YFY scheme ρ
0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] ∆ρ
0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] Figure 3. (colour on-line). The density distributions for a rigidly rotating barotropic configuration obtained with the YFY scheme(leftpanel), and the HSCF scheme(middle panel) as well as the relative differences thereof (right panel). The unit of density is given in cgs.The central density, equatorial radius, (uniform) specific entropy and ratio of rotational energy to gravitational one are ρ = 1 . × g cm − , R e = 2 . × km, s = 14 . B and T/ | W | = 3 . -0.352-0.35-0.348-0.346-0.344-0.342-0.34 0 50 100 150 200 250 E ne r g y F un c t i on Iteration Number 0.0001 0.001 0.01 0.1 1 0 50 100 150 200 250 V i r a l R e s i dua l Iteration number
Figure 4. (colour on-line). The values of the energy function (left panel) and of the virial residual (right panel) for the rigid rotationgiven in Fig. 3. The configuration settled to an rotational equilibrium iteration process is converged after 248 sweeps, where the energyfunction is longer lowered and the virial residual reaches < − , and the iteration is terminated there. The energy function is normalizedby G M ⊙ / R ⊙ , in which R ⊙ is the solar radius. The spikes that are seen in the energy function from time to time are due to smoothingsperformed to avoid false local minima. It is obvious that the EOS is generally baroclinic, i.e.,pressure depends not only on density but also entropyand chemical abundance in stars. It is, then, necessaryto handle such baroclinic equilibria for the study on thestructures of rotating stars. Such investigations are notmany in the literature, though (Uryu & Eriguchi 1994;Espinosa Lara & Rieutord 2007; Espinosa Lara & Rieutord2013; Fujisawa 2015; Rieutord et al. 2016). In the stellarevolution calculation, the shellular rotation is commonly as-sumed as a rotation law. It is not obvious, however, whetherthere exist such hydrostatic equilibria or not, especially forrapid rotation. Recently, Fujisawa (2015), one of the authorsof this paper, formulated a new self-consistent field schemethat is suitable for the construction of baroclinic equilibria.It is a Eulerian method and may be regarded as a (nontriv-ial) extension of the HSCF scheme we used in the previoussubsection. It is probably superior to the formulation pro-posed in this paper in terms of accuracy. The advantage ofthe YFY scheme, however, is that it is a Lagrangian methodand suitable to construct an evolutionary sequence out ofan ensemble of rotational equilibria. Anyway, with two in-dependent numerical schemes at hand to obtain baroclinic equilibria, it is a good opportunity for us to make compar-isons and validate both of them simultaneously.The baroclinicity is most easily introduced by an arti-ficial modification of the entropy distribution. We employthe polytrope-like EOS, p = Kρ γ with K being position-dependent,(i) spherical iso- K surfaces: K ( r ) = K (cid:18) e r R e (cid:19) , (37)(ii) oblate iso- K surfaces: K ( r ) = K (cid:26) e (cid:18) e P (cos θ ) (cid:19) r R e (cid:27) , (38)where, the Legendre polynomials are denoted by P n with θ being colattitude and employed to give non-sphericity to K distribution. Note that K corresponds to a (constant) spe-cific entropy in the original polytropic EOS and, hence, themodifications above may be regarded as an introduction ofentropy-dependence to pressure. Two dimensionless param-eters e and e specify radial non-uniformity whereas e givescolatitudinal dependency. We set e = 0 . e = 0 .
45 and e = 0 .
80 here. These models are highly artificial admittedly,and more realistic models, in which the entropy distributionsare determined self-consistently with energy generations and c (cid:13) , 000–000 otational equilibria by Lagrangian variational principle Z [ k m ] X [10 km]iteration number = 0 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 50 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 100 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 150 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 200 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 248 Figure 5. (colour on-line). The change of configuration in the searching process for the solution given in Fig.3.
YFY scheme ρ
0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] ∆ρ
0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] Figure 6. (colour on-line). Same as Fig. 3 but for the differential rotation obeying the law given in equation (36). The values ρ , R e , s ,and R p /R e are the same as those in Fig. 3. -0.4-0.398-0.396-0.394-0.392-0.39-0.388 0 50 100 150 200 250 E ne r g y F un c t i on Iteration Number 0.0001 0.001 0.01 0.1 1 0 50 100 150 200 250 V i r a l R e s i dua l Iteration number
Figure 7. (colour on-line). Same as Fig. 4 but for the differential rotation in Fig. 6. The iteration is terminated after 246 sweeps, atwhich point the energy functional is no longer lowered and the virial residual is sufficiently small < − .c (cid:13) , 000–000 N.Yasutake, K. Fujisawa, and S. Yamada Z [ k m ] X [10 km]iteration number = 0 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 50 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 100 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 150 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 200 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 246 Figure 8. (colour on-line). Same as Fig.5 but for the same differential rotation as in Figs. 6, 7. jYFY scheme 0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] ∆ j 0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] Figure 9. (colour on-line). Distribution of specific angular momentum for the same differential rotation as in Figs. 6-8. The angularmomentum is in unit of 10 cm rad s − . For this model, Ω = 1 . × − rad s − , and T/ | W | = 4 . radiative transfer, will be presented in the forthcoming pa-pers.In the FSCF scheme, the rotation law is given on theequatorial plane and the angular velocity distribution issolved self-consistently. We hence employ the rotation lawgiven in equation (36) only on the equatorial plane. The re-sultant configurations are expanded radially by 20 % andused as reference configurations for the YFY scheme just asin the previous sub section.We first present the results for case (i) with the spheri-cal isentropic surfaces given in equation (37). In Fig. 10, weshow isopycnic color-contours, in which the left panel is forthe YFY scheme and the middle panel for the FSCF scheme.Just as in the HSCF scheme, the FSCF scheme deals withnon-dimensional quantities normalized with the maximumdensity ρ and equatorial radius R e . To facilitate the com-parison with our scheme, we set ρ = 1 . × g cm − and R e = 2 . × km. The values of K and Ω chosen to be5 . × in cgs unit and 1 . × − rad s − , respectively.This value of Ω is a few hundred times the solar angularvelocity. We find again a reasonable agreement in most re-gions as is evident in the right panel of the figure, in whichthe relative difference is shown. Note that the number of thegrid points employed in the FSCF method, N r = 513 and N θ = 257, is much larger than that deployed in the YFYmethod, N = 489. It is also apparent from the same figurethat the deviation gets larger in the vicinity of the surfacefor the reason mentioned earlier.In Figs. 11 and 12 we present the values of the energyfunction and virial residual as a function of the number ofsweeps and the corresponding change of the configuration,respectively, just as in the barotropic case. It is observedthat the energy function is settled to a certain value after130 sweeps. It is apparent from the value of virial residual c (cid:13) , 000–000 otational equilibria by Lagrangian variational principle ∼ − that this is a false minimum. The true one is reachedat 261 sweeps after a couple of shuffles by the smoothings,which are marked by both major and minor spikes in theenergy function shown in the left panel of Fig. 12. Thisdemonstrates again that the virial residual is an importantdiagnostic measure for true minimum and the smoothing isan indispensable process to escape the trap of false minima.It is well known that the stability of rotational equilib-ria against axisymmetric perturbations can be judged by theSolberg and Høiland criterion which is summarized as fol-lows: stable configurations satisfy (i) the specific entropy s ,and (ii) the specific angular momentum j increases from thepoles to the equator on the isentropic surface. The rotationalconfiguration we obtained above satisfies these conditionsas can be understood from Fig. 13, in which we display theiso-angular momentum surfaces on top of the density color-contours; the arrows in the figure indicate the direction, inwhich the angular momentum increases.There is another well known rule we should care for: theBjerknes–Rosseland rule claims that if d Ω /dz is positive ata given point, the isobaric surface is more oblate than theisopycnic one there (Tassoul 1978). This means in particu-lar that if d Ω /dz is positive everywhere, matter temperaturesare lower at the poles than one on the equator on each iso-baric surface (Tassoul 1978). The rule is also confirmed forthe results presented above. The isobaric surfaces and iso- K (and, hence, iso-entropic) surfaces are shown in Fig. 14. Asis evident from the figure, matter on each isobaric surfacehas indeed a smaller value of K at poles than on the equator.On the other hand, the gradient of specific angular moment dj/dz , and hence d Ω /dz is also slightly positive as can beseen in Fig. 13.We move on to case (ii), i.e., with oblate isentropicsurfaces as specified in equation (38). What is importantwith this entropy distribution is that it leads to rota-tional equilibria with a shellular-type rotation demonstratedby (Fujisawa 2015). Figure 15 compares the isopycnic sur-faces obtained with the YFY scheme (left panel) and forthe FSCF scheme (middle panel). As quantified in the rightpanel, in which the relative difference is presented, the bothresults are in agreement within an error of ∼ V C < − is satisfied.We show in Fig. 18 the angular velocity distributionsobtained with the two methods for the same model. It is ob-vious that they are not cylindrical but of shellular-type . Therelative difference between the results obtained with the twoschemes is ∼ d Ω /dz is nega-tive everywhere in the current case with the oblate isentropicsurfaces as found in Fig. 18. As expected and confirmed in-deed in Fig. 20, the values of K are larger at poles than onthe equator this time. This model is hence also consistentwith the Bjerkeness-Rosseland rule. c (cid:13) , 000–000 N.Yasutake, K. Fujisawa, and S. Yamada
YFY scheme ρ
0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] ∆ρ
0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] Figure 10. (colour on-line). Same as Fig. 3 but for a baroclinic configuration with the spherical isentropic surfaces given in equation(37). The maximum density and equatorial radius normalization factors ρ and R e on FSCF scheme are set to ρ = 1 . × g cm − and R e = 2 . × km, respectively, and the constant K in equation (37) is chosen to be 5 . × in cgs unit. -0.287-0.286-0.285-0.284-0.283-0.282-0.281-0.28-0.279-0.278 0 50 100 150 200 250 300 E ne r g y F un c t i on Iteration Number 0.0001 0.001 0.01 0.1 1 0 50 100 150 200 250 300 V i r a l R e s i dua l Iteration number
Figure 11. (colour on-line). Same as Fig. 4 but for the baroclinic configuration presented in Fig. 10. The iteration is terminated at 261sweeps. Z [ k m ] X [10 km]iteration number = 0 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 52 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 104 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 156 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 208 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 261 Figure 12. (colour on-line). Same as Fig. 5 but for the baroclinic configuration presented in Figs. 10 and 11.c (cid:13) , 000–000 otational equilibria by Lagrangian variational principle YFY scheme s j 0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] ∆ j 0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] ∆ s 0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] Figure 13. (colour on-line). (Upper panels) iso-angular momentum contours on top of the entropy color maps for the baroclinicconfiguration shown in Figs. 10-11. The results for the YFY scheme and the FSCF scheme are displayed in the left and right panels,respectively. The arrows indicate the direction, in which the angular momentum increases. (Lower panels) the relative differences ofspecific angular momentum (left panel) and entropy (right panel) between the results obtained with the schemes.
YFY scheme PK 0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] ∆ P 0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] ∆ K 0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] Figure 14. (colour on-line). (Upper panels) iso- K surfaces on top of the pressure color-contours for the baroclinic configuration in Figs.10-13. The arrows indicate the direction, in which K increases. The left and right panels correspond to the results for the YFY andFSCF schemes, respectively. (Lower panels) the relative differences of pressure (left panel) and K (right panel) the results for the twoschemes.c (cid:13) , 000–000 N.Yasutake, K. Fujisawa, and S. Yamada
YFY scheme ρ
0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] ∆ρ
0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] Figure 15. (colour on-line). Same as Fig. 10 but for the oblate isentropic surfaces given in equation (38). The maximum density andequatorial radius are set to ρ = 1 . × g cm − and R e = 2 . × km, respectively. -0.288-0.287-0.286-0.285-0.284-0.283-0.282-0.281-0.28-0.279-0.278 0 50 100 150 200 250 300 E ne r g y F un c t i on Iteration Number 0.0001 0.001 0.01 0.1 1 0 50 100 150 200 250 300 V i r a l R e s i dua l Iteration number
Figure 16. (colour on-line). Same as Fig. 11 but for the configuration shown in Fig. 15. The iteration is terminated at 275 sweeps. Z [ k m ] X [10 km]iteration number = 0 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 55 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 120 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 165 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 220 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 275 Figure 17. (colour on-line). Same as Fig. 12 but for configuration displayed in Figs. 15 and 16.c (cid:13) , 000–000 otational equilibria by Lagrangian variational principle YFY scheme Ω
0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] ∆Ω
0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] Figure 18. (colour on-line). The left and middle panels show the angular velocity distributions of shellular-type obtained for theconfiguration shown as in Figs. 15 - 17 with the YFY and FSCF schemes, respectively. The right panel exhibits the relative difference.
YFY scheme s j 0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] ∆ j 0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] ∆ s 0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] Figure 19. (colour on-line). Same as Fig. 13 but for the configuration presented in Figs. 15-18.
In the application of the YFY scheme, we have always con-structed the reference configurations from the configurationsobtained with the HSCF or FSCF scheme, so far, expandingthem radially by a certain factor. This may raise a concernthat the reference configuration must be prepared in thisparticular way initially for our scheme. In order to demon-strate that this is not the case in fact, we try three moredifferent constructions and see if the final configurations arethe same or not.To ensure that all models have the same angular mo-mentum distribution in the final configuration, we still needto prepare the reference configurations by deforming this fi-nal configuration. Although it does not matter what config-uration is employed for the latter, we adopt the model with the shellular-type presented in Figs. 15, 18, 19 and 20. Thistime we do not expand but shrink the final configurationfor this model by 30 % in three different ways: (A) radially,(B) horizontally, and (C) vertically. To the three referenceconfigurations thus obtained, we apply the YFY scheme andsee if the same equilibrium, which should be also identicalto the configuration shown in Fig.17, is reached.We show the changes of the configurations with theMonte Carlo sweeps in Figs. 21-23, which correspond to de-formations (A)-(C) respectively. The upper left panels inthese figures exhibit the reference configurations, which areclearly different from the final configurations shown in thelower right panels. In particular, the reference configurationis prolate for case (B) whereas the rotational equilibriumis oblate owing to centrifugal forces. It is evident from thecomparison of the lower right panels in the figures that the c (cid:13) , 000–000 N.Yasutake, K. Fujisawa, and S. Yamada
YFY scheme PK 0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] ∆ P 0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] ∆ K 0 1 2 3X [10 km] 0 1 2 3 Z [ k m ] Figure 20. (colour on-line). Same as Fig. 14 but for the the oblate isentropic surfaces given in equation (38). resultant configurations after the optimizations look almostidentical with each other as well as with the original con-figuration given in Fig.17. We, hence, conclude that ourscheme is robust to the variation of reference configuration.It should be also noted that in the application to the stel-lar evolution calculation, we will be able to use the stellarstructure obtained at the previous time step as the referenceconfiguration, which should be hence pretty close to the trueequilibrium.It is one of the great advantages with the YFY schemebased on the Lagrange coordinates that we can record thehistory of the fluid element assigned to each node. This istrue not only for the Monte Carlo sweep as has been demon-strated so far but also for the true evolution in time when thescheme is applied to the 2-dimensional calculation of stellarevolution. Since nuclear burnings occur locally in each fluidelement, the chemical evolution is most easily treated onthe Lagrange coordinates and hence by the YFY scheme.It should be also mentioned that since the specific entropywould be conserved for each fluid element in the absence ofthe generation and transport of heat, radiative transport ofenergy will be formulated most unambiguously on the La-grange coordinates. This will be also the case for the trans-port of angular momentum. Convective transport of energyand angular momentum may be handled as well if it is jus-tified to treat them approximately as diffusions.The left three panels of Fig. 24 present the historiesof the values of the energy function for the cases shown inFigs. 21-23. The values of the energy function decrease inmost of the time in the Monte Carlo sweep. Sporadical smallglitches are due to the smoothings, which are necessary, as repeatedly mentioned, to ensure the convergence to the trueminimum, avoiding the trap of false minima.Shown in the right panels of Fig. 24 are histories of thevirial residuals for the same models. In the previous sections,we used the condition V C < − as one of the diagnostics tojudge if the equilibrium is reached as described in section 2.This condition works well also in the current cases as can beseen from the comparison of the right panels in Fig. 24 withFigs. 21-23. The importance of the criterion is also under-stood from Fig. 16, in which the virial residual is still large, V C > − , after ∼
120 iterations although the energy func-tional appears to have hit the minimum. The configurationis still out of equilibrium indeed, which can be confirmed inthe corresponding structure shown in Fig. 17.It should be noted, however, that the condition on thevirial residual alone is not sufficient for the judgement, sincethe virial residual is an integrated quantity and may becomesmall accidentally. For example, the configuration presentedin the upper-middle panel in Fig. 23 for the 190th itera-tion has not yet reached the equilibrium although the virialresidual satisfies V C < − as shown in Fig. 24. It is hencemandatory to take into account the convergence of energyfunction in addition to the condition on the virial residual.Considering the smoothing administered every ∼
100 sweeps(see the last paragraph of this section), we impose the con-dition that the energy function at the i − th iteration, E ( i ),should not lower than E ( i − i − i th step. Al-though this simple criterion seems to work well at least forthe models considered in this paper, it may need improve-ment.In this paper it is not our purpose to construct a se- c (cid:13) , 000–000 otational equilibria by Lagrangian variational principle quence of rotational configurations with different angularmomenta up to the mass shedding limit. As a matter offact, the Monte Carlo method employed in this paper tofind a configuration that gives the minimum to the energyfunctional is not suited very much for finding the limitingconfiguration for two reasons. Firstly, in the Monte Carlomethod many (non-equilibrium) configurations are tried toreach the true equilibrium and in so doing some trial con-figurations give energies lower than the true minimum er-roneously owing to numerical errors more often than not,which leads to a “numerical” mass shedding and no conver-gence. Secondly, the equilibrium configuration may not beunique for a given distribution of specific angular momen-tum on the Lagrangian coordinates. It is well known in factthat the angular momentum is not a monotonic function ofthe ratio of the polar to equatorial radii, R p /R e , near themass shedding limit (see, for example, Table 3 in Fujisawa(2015), in which the angular momentum increases initiallyas the ratio decreases but it decreases after the peak thatappears at R p /R e ∼ . R p /R e = 0 . R p /R e = 0 .
60, which we find is es-sentially the same as the result obtained with the FSCFmethod. Incidentally, this model is in the regime, where thesolution is unique in the above sense.At the end of section, we would like to discuss the im-portance of the smoothing. We employ the variational prin-ciple to obtain the hydrostatic equilibria. The idea is sim-ple, but the implementation is not indeed. Without an ap-propriate smoothing, Monte Carlo sweeps tend to make toodeformed triangles, which lead the generation of a false localenergy minimum. In Fig. 26 we show the result of the calcu-lation without the smoothing of acute-angled triangulations.There is no glitch in the panels of the energy function norof the virial residual, since we do not take into account thesmoothing. The energy function is settled after 1000 sweeps,although the virial residual does not satisfy V C < − . Thisresult suggests that the search of minimum energy is trappedin a local minimum. The resultant structure is far from thefinal results shown in Figs. 21-23. This is the reason why thesmoothing process is necessary.There will be many possible ways to alleviate highlyacute triangles, which could be the origin of false local min-ima. Followings are the simple measures we employed in thisstudy:(1) if there is a triangular cell that experiences more than30 % of change in area after a shift of a node duringthe Monte Carlo sweep, we move this particular node ofthe triangle to the center of the polygon formed with itsnearest neighboring nodes.(2) since the local shifts of some nodes in the first step some-times produce zigzag distributions of the nodes that areinitially located on one of the concentric surfaces (seethe construction of the triangulated mesh in the refer-ence configuration), we make them somewhat smootherby moving the grid on the top or bottom of the zigzag dis- tribution toward the middle of the two neighboring nodes;the first step is always followed by the second step.(3) in rare cases, the order of the grid points on the axisor on the equator is changed by their large shifts in theMonte Carlo sweep; then we restore the original orderby pulling back the grids in trouble toward the originalpositions appropriately; this measure is take only whenit is needed.(4) even if we do not find any large areal-change, we admin-ister the smoothing prescribed in (1) to all nodes period-ically; from our experience it should be done every ∼ ∼ the number of grid points / 5) sweeps; although thisis mainly meant to be a preventive measure, it is alsoeffective for the escape from a local minimum.The above prescriptions are admittedly empirical and de-pend on the construction of the triangulated mesh as well asthe direction of Monte Carlo sweep. Although there shouldbe more sophisticated and systematic ways to treat the falselocal minimum, we defer a further improvement to the fu-ture work, since the measures given above appear to workwell at least for the models considered in this paper. c (cid:13) , 000–000 N.Yasutake, K. Fujisawa, and S. Yamada Z [ k m ] X [10 km]iteration number = 0 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 105 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 210 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 315 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 420 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 526 Figure 21. (colour on-line). Same as Fig. 17 but for the reference configuration given by the deformation for case (A). Z [ k m ] X [10 km]iteration number = 0 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 109 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 218 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 327 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 436 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 548 Figure 22. (colour on-line). Same as Fig. 21 but for case (B).c (cid:13) , 000–000 otational equilibria by Lagrangian variational principle Z [ k m ] X [10 km]iteration number = 0 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 190 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 422 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 633 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 844 0 1 2 3 0 1 2 3 Z [ k m ] X [10 km]iteration number = 1057 Figure 23. (colour on-line). Same as Fig. 21 but for case (C).
We have developed a new formulation to obtain rotationalequilibria numerically. The scheme can handle not onlybarotropic but also baroclinic configurations, which arecritically important for the application to realistic stars.Such an achievement is itself a major break through tothe status quo, in which previous works are not manyand, more importantly, they are all based on the Euleriandescription (Uryu & Eriguchi 1994; Uryu & Eriguchi 1995;Espinosa Lara & Rieutord 2007; Espinosa Lara & Rieutord2013; Rieutord et al. 2016). Our method, on the otherhand, employs the Lagrangian description just as in one-dimensional counterparts for non-rotating stars and hencestellar evolution calculations, since we can trace the poten-tially complex movements of each fluid element. Our formu-lation is based on the variational principle, in which rota-tional equilibria are obtained as the configurations that op-timize the energy functional for given distributions of mass,specific entropy, and angular momentum on the Lagrangiancoordinates. In this paper all physical quantities are dis-cretized on triangulated grids.In order to validate our scheme, we compare the configu-rations obtained by our scheme with those by other Eulerianschemes: the HSCF scheme developed by Hachisu (Hachisu1986) for the barotropic configuration, and the FSCF schemeconceived by Fujisawa (Fujisawa 2015) more recently for thebaroclinic one. We have confirmed that all the configurationsincluding the ones with shellular-type rotations obtained inthis paper are linearly stable against axisymmetric pertur-bations and comply with the Bjorkness-Rosseland rule.In this comparison, we have reproduced the equilibriaa shellular-type rotation, which were obtained only recentlyby Fujisawa although the shellular rotation is commonly as- sumed in the one-dimensional stellar evolution calculationwith rotation being taken into account only approximately.It is found that the result obtained with the YFY schemeagree with those derived by other schemes with an error of5 ∼ %, which we think is reasonable, considering the rela-tively small number of nodes (489) we used in this paper.It is then a legitimate question to ask how the accuracy isimproved with the number of nodes. In order to address thisissue, we have conducted additional computations, changingthe node number from ∼
200 to ∼ N in our currentlyunparallelized code, where N is the node number.We find that the higher resolutions with ∼
800 and ∼ c (cid:13) , 000–000 N.Yasutake, K. Fujisawa, and S. Yamada -0.29-0.28-0.27-0.26-0.25-0.24-0.23 0 100 200 300 400 500 600 E ne r g y F un c t i on Iteration Number (A) 0.0001 0.001 0.01 0.1 1 0 100 200 300 400 500 600 V i r a l R e s i dua l Iteration number (A)-0.29-0.285-0.28-0.275-0.27-0.265-0.26-0.255 0 100 200 300 400 500 600 E ne r g y F un c t i on Iteration Number (B) 1e-05 0.0001 0.001 0.01 0.1 1 0 100 200 300 400 500 600 V i r a l R e s i dua l Iteration number (B)-0.288-0.286-0.284-0.282-0.28-0.278-0.276-0.274 0 200 400 600 800 1000 1200 E ne r g y F un c t i on Iteration Number (C) 1e-06 1e-05 0.0001 0.001 0.01 0.1 1 0 200 400 600 800 1000 1200 V i r a l R e s i dua l Iteration number (C)
Figure 24. (colour on-line). Same as Fig. 16 but for the reference configurations given by the deformation for cases (A)-(C). Z [ k m ] X [10 km] Figure 25. (colour on-line). The density distribution (left panel), and the edges and nodes (right panel) for shellular-type configurationfor R p /R e = 0 .
60 obtained with the YFY scheme. The result is essentially the same as the one derived with the FSCF scheme. The unitof density is given in cgs. The central density, equatorial radius are ρ = 1 . × g cm − , R e = 2 . × km, respectively.c (cid:13) , 000–000 otational equilibria by Lagrangian variational principle -0.29-0.28-0.27-0.26-0.25-0.24-0.23 0 200 400 600 800 1000 E ne r g y F un c t i on Iteration Number (A) 0.01 0.1 1 0 200 400 600 800 1000 V i r a l R e s i dua l Iteration number (A) Z [ k m ] X [10 km]iteration number = 1000 Figure 26. (colour on-line). The values of the energy function (left panel) and of the virial residual (middle panel) with initial model (A),but without the smoothing against acute-angled triangles. The right panel corresponds to the resultant structure at 1000-th iterationnumber. V i r a l R e s i dua l [ x - ] Total Node Number -0.287-0.286-0.285-0.284-0.283-0.282-0.281-0.28-0.279-0.278 0 500 1000 1500 2000 2500 3000 3500 E ne r g y F un c t i on Iteration Number 0.0001 0.001 0.01 0.1 1 0 500 1000 1500 2000 2500 3000 3500 V i r a l R e s i dua l Iteration number
Figure 27. (colour on-line). The range of converged virial residual for each resolution (left panel). The middle and the right panels showan example of the energy functional and the virial residual with 1073 nodes. the surface or to the rotational axis are the least contribu-tors, since they have either small densities or volumes; thismeans in turn that their positions are very difficult to ob-tain accurately; this may be understood if one scrutinizesthe right panels of Figs. 2, 3, 6, 9, 10, 15, 18 as well as thelower panels of Figs. 14, 19, 20; on the hand, these nodesgive contributions of the order of 10 − to the virial residual.In principle, if the number of nodes is sufficiently largeparticularly in the regions of our concern, then the virialequality should be improved. The current resolutions we canafford are way too small, though, as can be understood fromthe figures showing the node configurations. This may beremedied if we implement the multi-layer treatment.Multi-layer treatment may be useful to improve the ac-curacy , which will be also required in applying the newscheme to rotating stars in advanced evolutionary stages.As explained earlier, the numerical error comes mainly fromthe regions, whose contributions to the energy functionalare minor. If the multi-layer treatment is employed, theseminor regions can be treated separately from the major re-gions. Then the minor contributions are no longer minor.The important thing is that the variational principle canbe applied to each region individually, with other regionsserving as fixed gravitational potentials. In order to obtainthe global equilibrium, we need to find an equilibrium con- Such a treatment is already employed in their scheme by Es-pinosa Lara and Rieutord (Espinosa Lara & Rieutord 2013) and,combined with the spectral method, achieves very high accuracy.See also Kiuchi et al. (Kiuchi et al. 2010) figuration in each region consecutively and iteratively untilall these configurations do not change any longer simulta-neously. In addition to the minor regions mentioned above,the central region may also need to be treated with a specialcare, since the nodes are rather sparse there as can be seen,e.g., in the right panel of Fig. 1.We have demonstrated that our scheme is robust, ob-taining the same final configuration in equilibrium irrespec-tive of the reference configurations assumed initially. Wehave also observed, on the other hand, that the Monte Carlosweeps to get to the minimum of the energy functional areprone to be trapped by false local minima, which are gen-erated by deformations of the numerical grid and that thesmoothing should be properly administered to escape them.In order to distinguish the true minimum from a false one,the virial residual is found to be a good diagnostic. Althoughwe have not made any attempt in this paper to reconstructa triangulated mesh when it becomes too deformed, such re-gridding will be necessary when applying the YFY schemeto the evolution of rotating stars.The next step toward the application to the stel-lar evolution calculation will be to incorporate ra-diation transport in rotational equilibria. Unlike theprevious works (Uryu & Eriguchi 1994; Uryu & Eriguchi1995; Roxburgh 2006; Espinosa Lara & Rieutord 2007;Espinosa Lara & Rieutord 2013; Rieutord et al. 2016), wewill solve time-dependent diffusion equation on the same tri-angulated grid, which will not be a difficult task. The merid-ional circulation is another important ingredient (Zahn1992) when considering the application of our scheme to ac-tual rotating stars, since it is supposed to play a major role c (cid:13) , 000–000 N.Yasutake, K. Fujisawa, and S. Yamada in the re-distribution of angular momentum and elements(e.g. Mathis 2009). Although such motions cannot be han-dled directly with our scheme, they may be incorporatedeither as advections or as diffusions as just as in current 1Dcalculation.Convection is even more difficult to treat. The redis-tributions of angular momentum and elements as well asentropy in this case may be approximated as diffusions oradvections like the meridional circulation but the real diffi-culty with the convection is the fact that the convectivelyunstable configurations do not correspond to the minimumof the energy functional but to a saddle point. A couple ofideas are currently being tested and will be reported else-where.Since our scheme is applicable only to axisymmetricconfigurations, in which the specific angular momentum isconserved in each fluid element, extension of our formal-ism to 3D configurations such as triaxial equilibria (Tassoul1978) is much beyond the scope of this paper.Possible applications of our scheme are not limited tothe ordinary stars, though. They will be extended to e.g.,compact stars, proto-stars and proto-planets to mention afew. Incorporation of magnetic fields and/or general relativ-ity should be considered in due course.
ACKNOWLEDGEMENTS
We are grateful to K. Suzuki, T. Yamasaki, M. Okamoto,and T. Maruyama for fruitful discussions. This work waspartially supported by JSPS KAKENHI Grant Numbers25105510, 24244036 and Grant-in-Aid for Scientific Researchon Innovative Areas, No.24103006.
APPENDIX A: TRANSFORMATION MATRIX
In this appendix we give the transformation matrix T lm thatis used in the expression of the discretized energy functional,equation (24). It is actually introduced in association withthe interpolation of physical quantities in each triangularcell. Suppose that a quantity X is given by on all nodes andwe interpolate them in each triangular cell as X = α X z + α X ̟ + α X , (A1)in which ̟ and z are the cylindrical coordinates of an arbi-trary point in the cell and α ’s are numerical coefficients to bedetermined. Evaluating the above equation at three apexes(labeled with l ) of the n -th triangular cell, one obtains thefollowing equations: X l ( n ) = T lm ( n ) α Xm ( n ) , (A2)where l, m = 1 , ,
3. The matrix T ( n ) in this expression isreferred to the transformation matrix in this paper, in whichdefined as T ( n ) = z ( n ) ̟ ( n ) 1 z ( n ) ̟ ( n ) 1 z ( n ) ̟ ( n ) 1 ! . (A3)This matrix is useful in the finite element method employedhere. For example, the volume of the triangular cell V ∆ isevaluated as, Figure A1.
The diffeomorphic mapping that gives the Jacobian. V ∆ = Z cell π̟ d̟ dz = 2 π · | J | · ̟ + ̟ + ̟ . (A4)Here | J | is the Jacobian for the diffeomorphic mapping f of an isosceles right triangle constructed by the unit fun-damental vectors e x , e y of the two-dimensional Cartesiancoordinates to the triangle with three apexes X X X asshown in Fig. A1. The Jacobian is then equal to the deter-minant | T ( n ) | . This is obtained from the following relationbetween two forms: Here, the areal element under the map-ping is shown as du ∧ dv = J ( f ) dx ∧ dy, (A5)where u and v are the vectors −−−→ X X and −−−→ X X , respectively. | T ( n ) | / ̟ + ̟ + ̟ ) / ̟ -coordinate of thegeometrical centre of the triangule. REFERENCES
Burbidge K. M., Burbidge G. R., Fowler W. A., Hoyle F.1957, Review of Modern Physics, 29, 4Cameron, A. G. W., 1957, Chalk River Report, CRL-41Chandrasekhar, S. 1969, The Silliman Foundation Lec-tures, New Haven: :
Ellipsoidal figures of equilibrium , YaleUniversity PressEndal, A. S., & Sofia, S. 1976, ApJ, 210, 184Espinosa Lara F., Rieutord M., 2007, A&A, 470, 1013—, 2013, A&A, 552, A35Eriguchi Y., Mueller E., 1985, A&A, 146, 260Friedman, J. L., & Schutz, B. F. 1978, ApJ, 221, 937Fukuda, I. 1982, PASP, 94, 271Fujisawa, K. 2015, MNRAS 454, 3060Hachisu, I. 1986, ApJS, 61, 479Hayashi, C., 1961, PASJ, 13, 442: Hayashi, C., 1961,PASJ, 13, 450.Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528,368Henyey, L. G., Forbes, J. E., & Gould, N. L. 1964, ApJ,139, 306Kiuchi, K., Nagakura, H., & Shoichi, Y. 2010, ApJ, 717,666Maeder, A., & Meynet, G. 2000, Astron. Astrophys., Suppl.Ser., 38, 143Maeder, A., Meynet, G., Lagarde, N., & Charbonnel, C.2013, A&A, 553, A1 c (cid:13) , 000–000 otational equilibria by Lagrangian variational principle Mathis, S. 2009, A&A, 506, 811Mathis, S. 2013, in Goupil M., Belkacem K., Neiner C.,Lignieres F., Green J. J., eds, Lecture Notes in Physics,Vol. 865, Seismology for Studies of Stellar Rotation andConvection. Springer-Verlag, Berlin, p. 49Mathis S., Palacios A., Zahn J.-P., 2004, A&A, 425, 243Mathis, S., Decressin, T., Eggenberger, P., & Charbonnel,C. 2013, A&A, 558, A11Meynet G. & Maeder A. 1997, A&A, 321, 465Ostriker, J. P., & Mark, J. W.-K. 1968, ApJ, 151, 1075Potter A. T., Tout C. A., Brott I., 2012a, MNRAS, 423,1221Potter A. T., Tout C. A., Eldridge J. J., 2012b, MNRAS,419, 748Rieutord, M., Espinosa Lara, F., Putigny, B. 2016, J. Com-putational Phys., 318, 277Roxburgh, I. W., 2006, A&A, 454, 883Takahashi K., Umeda H., Yoshida T., 2014, ApJ, 794, 40Talon, S., Zahn J.-P., Maeder A., Meynet G., 1997, A&A,332, 209Tassoul, J.-L. 1978, Princeton Series in Astrophysics,Princeton:
Theory of rotating stars , University Press, andreferences thereinUryu, K., & Eriguchi, Y. 1994, MNRAS, 269, 24—, 1995, MNRAS, 277, 1411Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Reviewsof Modern Physics, 74, 1015, and references thereinYasutake, N., Fujisawa K., Yamada S. 2015, MNRAS let-ter, 446, L56Zahn, J.-P., 1992, A&A, 265, 115 c (cid:13)000