# Toroidal figures of equilibrium from a second-order accurate, accelerated SCF method with subgrid approach

J.-M. Huré, F. Hersant

► To cite this version:

J.-M. Huré, F. Hersant. Toroidal figures of equilibrium from a second-order accurate, accelerated SCF method with subgrid approach. Monthly Notices of the Royal Astronomical Society, 2017, 464 (4), pp.4761 - 4776. 10.1093/mnras/stw2569. hal-01384481

HAL Id: hal-01384481

https://hal.science/hal-01384481v1

Submitted on 24 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Toroidal figures of equilibrium from a second-order accurate, accelerated SCF method with subgrid approach

J.-M. Huré<fnref n="1" />,<fnref n="2" />* and F. Hersant<fnref n="1" />,<fnref n="2" />

Accepted 2016 October 3. Received 2016 July 28; in original form 2016 June 23

## ABSTRACT

We compute the structure of a self-gravitating torus with polytropic equation of state (EOS) rotating in an imposed centrifugal potential. The Poisson solver is based on isotropic multigrid with optimal covering factor (fluid section-to-grid area ratio). We work at second order in the grid resolution for both finite difference and quadrature schemes. For soft EOS (i.e. polytropic index $n \ge 1$), the underlying second order is naturally recovered for boundary values and any other integrated quantity sensitive to the mass density (mass, angular momentum, volume, virial parameter, etc.), i.e. errors vary with the number $N$ of nodes per direction as $\sim 1/N^2$. This is, however, not observed for purely geometrical quantities (surface area, meridional section area, volume), unless a subgrid approach is considered (i.e. boundary detection). Equilibrium sequences are also much better described, especially close to critical rotation. Yet another technical effort is required for hard EOS ($n < 1$), due to infinite mass density gradients at the fluid surface. We fix the problem by using kernel splitting. Finally, we propose an accelerated version of the self-consistent field (SCF) algorithm based on a node-by-node pre-conditioning of the mass density at each step. The computing time is reduced by a factor of 2 typically, regardless of the polytropic index. There is a priori no obstacle to applying these results and techniques to ellipsoidal configurations and even to 3D configurations.

**Key words:** equation of state – gravitation – methods: analytical – methods: numerical.

## 1 INTRODUCTION

The theory of figures of equilibrium is a broad and ancient astrophysical subject, aiming at understanding the flattening of the Earth and stars due to proper spin and tides (Ivory 1831; Chandrasekhar & Lebovitz 1962; James & Kopal 1963; Clement 1967; Chandrasekhar 1969). While very complex structures and shapes are permitted, matter can basically attain two configurations, depending on angular momentum and rotation profile (Chandrasekhar 1969; Hachisu 1986; Horedt 2004): an ellipsoidal shape or a toroidal shape at high momentum. This latter class is of great interest since discs and rings are seen in various astrophysical environments (Barker 2001), around planets and satellites (Goertz & Ip 1982), around young and evolved stars (Washimi, Shibata & Mori 1996; Dutrey et al. 2014), and around stellar and supermassive black holes (Rouan et al. 1998; Kishimoto et al. 2011; Elvis 2012). The subject is less documented (Kowalewsky 1885; Dyson 1893; Wong 1974; Hachisu 1986; Odrzywołek 2003; Petroff & Horatschek 2008; Slaný et al. 2013), probably for technical reasons (i.e. the strong deviation with respect to sphericity). For compact systems subject to rapid

evolutions, the framework of general relativity is pertinent (Kuwahara 1988; Tohline & Hachisu 1990; Nishida, Eriguchi & Lanza 1992; Nishida & Eriguchi 1994; Abramowicz, Karas & Lanza 1998; Cai & Shu 2003; Stuchlík, Slaný & Kovář 2009; Font, Montero & Shibata 2010; Hamerský & Karas 2013). Most investigations of self-gravitating tori are, however, performed in classical gravity. Configurations have been computed in detail by Hachisu (1986) and collaborators (see also Eriguchi 1978; Eriguchi & Mueller 1993; Hashimoto et al. 1993; Lai, Rasio & Shapiro 1993; Kley 1996; Horedt 2004). Despite restrictive assumptions (polytropic gas, symmetries, etc.), there is a wide and rich collection of possible equilibria (Ansorg, Kleinwächter & Meinel 2003). This paper revisits the determination of single-body, toroidal configurations from the self-consistent field (SCF) method. Special attention is paid to the convergence process leading to the physical solution and to the preservation of the accuracy order, which are rarely discussed in the literature (Odrzywołek 2003). Here, we work at second order in the grid resolution. Following the pioneering work by Lanza (1987), the Poisson solver is based on multigrid. Indirectly, the polytropic equation of state (EOS) plays a critical role. For instance, we show that the fluid volume – a key quantity when building equilibrium sequences – appears quite uncertain due to curvature effects at the inner and outer edges. This is fixed by a subgrid

<footnote n="1">¹Univ. Bordeaux, LAB, UMR 5804, F-33615 Pessac, France</footnote>

<footnote n="2">²CNRS, LAB, UMR 5804, F-33615 Pessac, France</footnote>

<footnote>* E-mail: jean-marc.hure@u-bordeaux.fr</footnote>
approach. There is an additional difficulty for $n < 1$: the presence of infinite gradients of the mass density at the fluid surface makes standard quadrature schemes imprecise. This is solved by treating the underlying power-law behaviour separately. Globally, the fluid boundary needs to be accounted for in order to fully recover the order of underlying schemes, i.e. second order. The virial parameter is then a good indicator of the precision of the solution and integrated values. Finally, we propose a semi-empirical recipe to speed up the SCF loop by a factor up to 2 typically. It is a mass density pre-conditioner.

The paper is organized as follows. In Section 2, we list the assumptions, review the equation set and define the enthalpy associated with the polytropic EOS. We briefly comment on the existence of solutions, and in particular how the index $n$ shapes the mass density profile at the fluid surface and below. We put the problem into dimensionless form and write down the virial equation, associated kernels and the virial parameter. In Section 3, we express the three constants of the problem and recall the principle of the SCF method, including self-normalization (which permits us to isolate a particular solution). We briefly outline the multi-grid method used to solve the Poisson equation, and define the associated nested grids, from the coarsest level to the finest one. The formula for boundary values (BVs) is given. Then comes the discretization of the Poisson equation and a short error analysis, which gives the optimal multigrid level. We discuss how the polytropic hypothesis impacts on the accuracy of quadratures at second order. Section 4 is devoted to examples obtained under rigid rotation for both hard and soft EOS, with index $n \in [0.5, 1, 1.5]$, at a low resolution of 128 $\times$ 128. We show how the gravitational, centrifugal and internal energy densities are distributed within the fluid section. The effect of various parameters (number of nested grids, covering factor) on the convergence process is discussed. In Section 5, we describe the subgrid approach and kernel splitting, which are necessary to go beyond the raw model. A few examples demonstrate its efficiency. In Section 6, we show that the number of SCF iterations can be reduced by an appropriate pre-conditioning of the mass density field inside the loop. This gain (a factor of 2 typically) directly concerns the computing time and is therefore important when performing high-resolution models. The last section is devoted to concluding remarks.

# 2 PHYSICAL MODEL

## 2.1 Hypothesis

We study the equilibrium of an isolated, compressible fluid with pressure $P$ and mass density $\rho$ rotating in its own gravitational potential $\Psi$ at rotation rate $\Omega$. We assume no accretion, no internal circulation and no viscosity. The EOS is of the form $f(P, \rho) = 0$. Barotropes (and in particular polytropes) enable us to cover various situations, e.g. perfect gas, gas with relativistic electrons, degenerate gas, convective media, etc. (Cox 1968). The analysis is performed in cylindrical coordinates $(r, \phi, Z)$ with vector basis $(e_R, e_\phi, e_Z)$, and the rotation axis is along $e_Z$, i.e. $\Omega = \Omega e_Z$.

The meridional section of the fluid is denoted by $S$, and it intercepts the fluid surface through a closed curve $\Gamma(r_\Gamma)$. We define a unit vector $e_\Gamma$ locally normal to $\Gamma$ and oriented outwards and the associated coordinate is

$$q - q_{\Gamma} = (r_{\Gamma} - \mathbf{r}) \cdot e_{\Gamma}, \qquad (1)$$

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Configuration for the self-gravitating torus (limited to the half-plane $Z \ge 0$). Points A and B are the radii at the inner and outer edges, respectively, while point M is an interior point. The meridional section is denoted by $S$, and $e_\Gamma$ is the unit vector normal to $\Gamma$, the fluid surface.</figcaption></figure>

with $q = q_\Gamma$ on to $\Gamma$ and $q > q_\Gamma$ inside the fluid (i.e. below the surface). Under axial and equatorial symmetries, we have

$$\left\{ \begin{array}{l} \partial_{\phi} f = 0, \\ f(Z) - f(-Z) = 0, \end{array} \right. \qquad (2)$$

for any scalar quantity $f$. The system and notations are summarized in Fig. 1. We focus on toroidal configurations, where matter has deserted the rotation axis.

## 2.2 Equation set

According to the Schwarz condition, the angular momentum and subsequently the rotation rate are independent of the Z-coordinate (Amendt, Lanza & Abramowicz 1989). The circular velocity is of the form $v = \Omega(R)Re_\phi$ and the equilibrium of the self-gravitating fluid is described by the equation set

$$\left\{ \begin{array}{l} -\frac{v^2}{R} e_R = -\frac{1}{\rho} \nabla P - \nabla \Psi, \\ f(P, \rho) = 0, \\ \Delta \Psi = 4\pi G \rho. \end{array} \right. \qquad (3)$$

The Poisson equation, i.e. equation (3 c), is subject to boundary conditions. BVs are accessible uniquely from the generalized Newton's law

$$\Psi = -G \iiint_{\text{fluid}} \frac{\rho(r')d^3r'}{|r - r'|}, \qquad (4)$$

where $r$ and $r'$ are spherical vectors. If we define the centrifugal potential $\Phi$ as

$$\Phi = - \int \Omega^2(R)RdR < 0, \qquad (5)$$

then equation (3 a) takes the form

$$\frac{1}{\rho} \nabla P + \nabla(\Psi + \Phi) = 0. \qquad (6)$$

The coupling between the rotation rate and the spatial distribution of matter is a tricky problem, which exceeds the present context (e.g. Binney & Tremaine 1987). In Balmforth et al. (1992), the link is made through conservation of potential vorticity. In general however, $\Omega$ remains an input profile. Various options are classically
considered (Bodenheimer & Ostriker 1973; Eriguchi & Mueller 1985, 1993; Hachisu 1986; Reese 2006), including power laws of the radius (e.g. Slaný et al. 2013). We will not favour a particular profile, except for illustrative purposes.

## 2.3 Enthalpy field

Equation (6) can be fully integrated provided the pressure term is the gradient of a scalar field $H$, the enthalpy, i.e.

$$ \frac{1}{\rho} \nabla P \equiv \nabla H, \qquad (7) $$

which is possible for a polytropic EOS

$$ P = K \rho^{\gamma}, \qquad (8) $$

where $K$ and $\gamma$ are positive constants. We thus have

$$ H = K \frac{\gamma}{\gamma - 1} \rho^{\gamma - 1} + \text{cst}, \qquad (9) $$

where the constant is generally forgotten. So, the equation set becomes

$$ \left\{ \begin{array}{l} \nabla (\Phi + H + \Psi) = 0, \\ H = K(n+1)\rho^{1/n}, \\ \Delta\Psi = 4\pi G\rho, \end{array} \right. \qquad (10) $$

where the polytropic index $n \ge 0$ is defined by

$$ \gamma = 1 + \frac{1}{n}. \qquad (11) $$

For isothermal structures where $\gamma = 1$ and possibly for layered systems made of different polytrope s where continuity is required (Singh 1970; Murphy & Fiedler 1985; Caimmi 1986; Beech 1988; Rucinski 1988; Curry & McKee 2000; Horedt 2004), equation (9) still holds but the integration constant must be accounted for. Actually, with the definition (Kimura 1981)

$$ H = K \frac{\gamma}{\gamma - 1} \left( \rho^{\gamma - 1} - \rho_0^{\gamma - 1} \right), \qquad (12) $$

where $\rho_0$ is a constant, we have

$$ \lim_{\gamma \to 1} H = K \ln \frac{\rho}{\rho_0}. \qquad (13) $$

The incompressible approximation $n = 0$ is also of great interest (Chandrasekhar 1969). In this case, equation (10 b) is obsolete because $\rho$ is uniform in the system throughout and disconnected from the enthalpy field.

## 2.4 Comment on the existence of solutions

We see from equation (10 a) that $H + \Psi + \Phi$ is a space invariant, namely

$$ H = C - \Psi - \Phi, \qquad (14) $$

where $C$ is a constant. Because $\Psi$ and $\Phi$ are both negative, and due to the correlation between $\Psi$ and $\rho$ (the larger the mass density, the deeper the potential well), $H$ is expected to rise inside matter, with one maximum (or more). The regions occupied by the fluid (where $\rho > 0$) are dictated by the value of constant $C$ (to be compared with Jacobi's constant in the restricted three-body problem). When $C > 0$, matter fills the whole physical space. So, finite-size systems require $C < 0$.

Clearly, equation (10 b) allows for a large variety of mathematical solutions, depending on the polytropic index. For some indices,

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Expected density profile (bold) and derivative (red arrow) when crossing the fluid boundary $\Gamma$ normally (coordinate $q$), depending on the polytropic index $n$ (see also Fig. 1). The enthalpy is shown (blue).</figcaption></figure>

$\rho$ can be negative or complex. We retain the physical solutions associated with $H > 0$ for any $n$. This condition is explicitly stated by replacing equation (10 b) by

$$ K(n+1)\rho^{1/n} = \sup(H, 0). \qquad (15) $$

This cut-off is a major difference with the approach based on the Lane–Emden equation, apart from the rotation term. Besides, we focus on single bodies while $H > 0$ can occur in two disconnected regions (or more), leading to multibody configurations (see e.g. Hachisu, Eriguchi & Nomaoto 1986).

## 2.5 Fluid surface and interior

If the fluid is isolated, the surface is defined by $\rho = 0$. This corresponds to the zero-enthalpy level. In the isothermal limit (and possibly for other particular indices), this surface may not exist at all, or be rejected at infinity, which case is well known for spherical configurations (Schmitz 1990; Horedt 2004). Just below the boundary, $\rho$ rises according to equation (15). Perpendicularly to $\Gamma$ (see Section 2.1 and Fig. 1), the mass density gradient is

$$ \nabla_q \rho \sim n \nabla_q H \times \rho^{\frac{n-1}{n}}, \qquad (16) $$

which, by linearly expanding $H$ around $q = q_\Gamma$, is given by

$$ \nabla_q \rho \sim n(q - q_\Gamma)^{n-1}. \qquad (17) $$

As a result, the mass density gradient is finite at the fluid boundary for $n \ge 1$ (soft EOS), but infinite for $n < 1$ (hard EOS). This latter case introduces a difficulty from a numerical point of view (see below). Since the gravitational potential is a continuous function of space (Kellogg 1929; Durand 1953) as well as the rotation law, $H$ is also continuous, which property is transmitted to $\rho$ inside the fluid, from equation (14). Fig. 2 summarizes four different situations: the larger the index, the more peaked the mass density, and the wider the wings. Conversely, the closer to zero the index, the sharper the transition from the fluid to the external medium, with a bare jump for $n = 0$.
## 2.6 Dimensionless problem

It is convenient to work with adimensioned quantities, especially for a numerical approach. We use the generic notation $f_0 \hat{f} \equiv f$, where $f_0$ is the magnitude of $f$, which leaves $\hat{f}$ of the order of unity. Space coordinates are normalized by the physical length-scale $L$, i.e.

$$ \left\{ \begin{aligned} \hat{R} &\equiv \frac{R}{L}, \\ \hat{Z} &\equiv \frac{Z}{L}, \end{aligned} \right. \qquad (18) $$

and we set

$$ \left\{ \begin{aligned} \hat{\rho} &\equiv \frac{\rho}{\rho_0}, \\ \hat{P} &\equiv \frac{P}{P_0}, \\ \hat{\Omega} &\equiv \frac{\Omega}{\Omega_0}, \\ \hat{\Psi} &\equiv \frac{\Psi}{\Psi_0}, \\ \hat{H} &\equiv \frac{H}{H_0}, \\ \hat{\Phi} &\equiv \frac{\Phi}{\Phi_0}, \end{aligned} \right. \qquad (19) $$

for the mass density, pressure, rotation rate, gravitational potential, enthalpy and centrifugal potential, respectively, where

$$ \left\{ \begin{aligned} P_0 &= K \rho_0^\gamma, \\ \Psi_0 &= G \rho_0 L^2, \\ H_0 &= K(n+1)\rho_0^{1/n}, \\ \Phi_0 &= \Omega_0 L^2. \end{aligned} \right. \qquad (20) $$

The dimensionless equation set then writes

$$ \left\{ \begin{aligned} \hat{\Psi} + C_1 \hat{H} + C_2 \hat{\Phi} &= C_3, \\ \hat{\rho}^{1/n} &= \sup(\hat{H}, 0), \\ \Delta \hat{\Psi} &= 4\pi \hat{\rho}, \end{aligned} \right. \qquad (21) $$

where

$$ \hat{\Phi} = - \int \hat{\Omega}^2(\hat{R}) \hat{R} d\hat{R}, \qquad (22) $$

$$ \hat{\Psi} = - \iiint_{\text{fluid}} \frac{\hat{\rho}(\hat{r}')}{|\hat{r} - \hat{r}'|} d^3 \hat{r}', \qquad (23) $$

for BVs, and the constants are

$$ \left\{ \begin{aligned} C_1 &= \frac{K(n+1)\rho_0^{\gamma-2}}{G L^2}, \\ C_2 &= \frac{\Omega_0^2}{G \rho_0}, \end{aligned} \right. \qquad (24) $$

for finite indices. For $\gamma = 1$, we have

$$ \left\{ \begin{aligned} H_0 &= K, \\ \hat{\rho} &= e^{\hat{H}}, \\ C_1 &= \frac{K}{G \rho_0 L^2}. \end{aligned} \right. \qquad (25) $$

While notations and definitions can vary from one author to another, these formulae are standard. Note that constant $C_1$ does not exist in Hachisu (1986), since it is included into the definition of $\hat{H}$.

## 2.7 Virial theorem and associated test

The virial theorem is established by projecting the vectorial Euler equation, i.e. equation (3 a), on the position vector $\mathbf{r}'$ and integrating over the mass distribution. It is an energy equation which shows how the gravitational $W$, internal $U$ and kinetic $T$ energies are spatially

distributed (e.g. Cox 1968). For the actual problem where only the gravitational and pressure forces are present, the theorem reads

$$ W + 2T + U = 0, \qquad (26) $$

where

$$ \left\{ \begin{aligned} W &= \frac{1}{2} \iint \Psi d\mathbf{m}, \\ T &= \frac{1}{2} \iiint v^2 d\mathbf{m}, \\ U &= 3\Pi, \\ \Pi &= \iint \frac{P}{\rho} d\mathbf{m}. \end{aligned} \right. \qquad (27) $$

When other forces are acting in volume or at the surface, the virial equation must be modified. This is for instance the case for loaded polytropes (Fang, Xiang & Wang 1983), with magnetic fields (Yoshida & Eriguchi 2006), when the fluid is overpressurized (Umemura & Ikeuchi 1986), for multibody configurations (Taniguchi & Nakamura 2000), etc. Given axial symmetry and equation (24), we easily show that equation (26) becomes

$$ \hat{W} + \frac{C_1}{n+1} \hat{U} + 2C_2 \hat{T} = 0, \qquad (28) $$

for finite indices $n$ (the second term is just $C_1 \hat{U}$ for $\gamma = 1$), where

$$ \left\{ \begin{aligned} W &= GL^5 \rho_0^2 \hat{W}, \\ \Pi &= P_0 L^3 \hat{\Pi}, \\ T &= L^5 \rho_0 \hat{\Omega}_0^2 \hat{T}, \end{aligned} \right. \qquad (29) $$

and

$$ \left\{ \begin{aligned} \hat{W} &= \pi \iint_S \hat{\Psi} \hat{\rho} \hat{a} d\hat{a} d\hat{z}, \\ \hat{U} &= 3 \hat{\Pi}, \\ \hat{\Pi} &= 2\pi \iint_S \hat{P} \hat{a} d\hat{a} d\hat{z}, \\ \hat{T} &= \pi \iint_S \hat{\Omega}^2 \hat{\rho} \hat{a}^3 d\hat{a} d\hat{z}. \end{aligned} \right. \qquad (30) $$

Interestingly enough, equation (28) can be written in the form

$$ 2\pi \iint_S \kappa_{\text{Virial}} d\hat{a} d\hat{z} = 0, \qquad (31) $$

where

$$ \kappa_{\text{Virial}} = \kappa_W + 2\kappa_T + \kappa_U \qquad (32) $$

is the energy per unit meridional section area and

$$ \left\{ \begin{aligned} \kappa_W &= \frac{1}{2} \hat{\rho} \hat{\Psi} \hat{a}, \\ \kappa_U &= \frac{3C_1}{n+1} \hat{P} \hat{a}, \\ \kappa_T &= \frac{1}{2} C_2 \hat{\Omega}^2 \hat{\rho} \hat{a}^3, \end{aligned} \right. \qquad (33) $$

with $\kappa_U = 3C_1 \hat{P} \hat{a}$ in the isothermal case. This kernel varies in space: it is positive where dispersion forces dominate, and negative when cohesive (i.e. gravitational) forces dominate.

When equation (21) is solved numerically, equation (28) is expected to be satisfied too. Since numerical methods have not infinite precision but produce small errors, the virial equation is generally not strictly zero. The self-consistency of a solution can then be checked a posteriori by comparing the left-hand side of equation (28) to its largest term. Relative to the gravitational term, the virial test writes

$$ \mathrm{VP} \stackrel{?}{=} 0, \qquad (34) $$
where the virial parameter is

$$VP = \frac{1}{|\hat{W}|} \left( \frac{C_1}{n+1} \hat{U} + 2C_2 \hat{T} \right) - 1. \quad (35)$$

This test does not prove the exactness of the computed solution, but only its self-consistency. It must also be interpreted with caution since, in numerical calculus, small quantities are difficult to determine with precision.

### 3 SOLVING THE NON-LINEAR PROBLEM AT SECOND ORDER

We see that $\hat{\rho}$ (or $\hat{H}$) is the solution of a non-linear, integro-algebraic equation. Given $L$, $\rho_0$, $\Omega_0$, $K$ and $n$, one can compute $C_1$ and $C_2$ from equation (24), but constant $C_3$ remains undetermined unless $\hat{\Psi}$, $\hat{H}$ and $\hat{\rho}$ are known at some point in space. In contrast, if the three constants are set arbitrarily, the solution in the form of $\hat{\Psi}$, $\hat{H}$ and $\hat{\rho}$ may not exist. Actually, for three reference points A, B and M of space, these constants must satisfy

$$\left\{ \begin{aligned} \hat{\Psi}_A + C_1 \hat{H}_A + C_2 \hat{\Phi}_A &= C_3, \\ \hat{\Psi}_B + C_1 \hat{H}_B + C_2 \hat{\Phi}_B &= C_3, \\ \hat{\Psi}_M + C_1 \hat{H}_M + C_2 \hat{\Phi}_M &= C_3, \end{aligned} \right. \quad (36)$$

which is not automatic. By selecting points A and B on to $\Gamma$ where the enthalpy is zero, we get

$$\left\{ \begin{aligned} C_1 &= -\frac{\hat{\Psi}_M \Delta \hat{\Phi}_{AB} + \hat{\Psi}_A \Delta \hat{\Phi}_{BM} + \hat{\Psi}_B \Delta \hat{\Phi}_{MA}}{\hat{H}_M \Delta \hat{\Phi}_{AB}}, \\ C_2 &= -\frac{\Delta \hat{\Psi}_{AB}}{\Delta \hat{\Phi}_{AB}}, \\ C_3 &= \frac{\hat{\Psi}_B \hat{\Phi}_A - \hat{\Phi}_B \hat{\Psi}_A}{\Delta \hat{\Phi}_{AB}}, \end{aligned} \right. \quad (37)$$

with $\Delta \hat{\Phi}_{AB} = \hat{\Phi}_A - \hat{\Phi}_B$ and so on for $\Delta \hat{\Phi}_{BM}$ and $\Delta \hat{\Phi}_{MA}$.

#### 3.1 Principle of the SCF method

To get the equilibrium, the three constants are often regarded as unknowns and solved together with the enthalpy, mass density and gravitational potential through an iterative procedure. A common method consists in computing successively (i) the mass density from the enthalpy through equation (21 b), (ii) the gravitation potential from equation (21 c) and the centrifugal potential from equation (22), (iii) constants $C_1$, $C_2$ and $C_3$ from equations (37), and (iv) the enthalpy from equation (21 a), and so on until input and output match. This is the so-called self-consistent field (SCF) method (Ostriker & Mark 1968; Hachisu 1986). It requires a starting guess $\hat{H}^{(0)}$ and a convergence criterion to detect, at some step $t$ within the cycle, if the enthalpy $\hat{H}^{(t)}$ has converged in the numerical sense. The capacities of convergence of the SCF method are seemingly misunderstood (Odrzywołek 2003), but it happens to work in most cases if (i) the three reference points are conveniently chosen, and (ii) the enthalpy is appropriately rescaled within the cycle. Traditionally, points A and B are radii at the inner and outer edge of the torus, respectively. Since $\hat{\Omega}$ is generally monotonic with the radius, this ensures that $\Delta \hat{\Phi}_{AB}$ is maximum.

#### 3.2 Enthalpy normalization

As outlined in Hachisu (1986), $\hat{H}$ does not necessarily remain under control. From equations (21 a) and (21 b), we have

$$\hat{\rho}^n = \frac{C_3 - C_2 \hat{\Phi} - \hat{\Psi}(\hat{\rho})}{C_1}, \quad (38)$$

and we see that both sides respond in the same way when $\hat{\rho}$ varies. The SCF cycle is therefore prone to runaway. However, since neither $C_2$ nor $C_3$ depends on point M, we have an extra degree of freedom to avoid any drift. We can impose that the enthalpy at some point in the fluid takes a particular value (see Section 3.2). By inserting $C_1^{(t)}$ into equation (43), we get

$$\hat{H}^{(t+1)} = \hat{H}_M^{(t)} \times \frac{C_3^{(t)} - C_2^{(t)} \hat{\Phi}^{(t)} - \hat{\Psi}^{(t)}}{C_3^{(t)} - C_2^{(t)} \hat{\Phi}_M^{(t)} - \hat{\Psi}_M^{(t)}}, \quad (39)$$

meaning that if M is held fixed, its enthalpy is constant during the SCF iterations. Then, at convergence (i) $\hat{H}_M$ has conserved its initial value, and (ii) the maximum enthalpy is not under control, and it does not necessarily occur at point M. This is not an obstacle to convergence, but there are clearly an infinity of solutions, each dictated by $\hat{H}_M^{(0)}$. We follow Hachisu (1986) by requiring

$$\left\{ \begin{aligned} \hat{H}_M &= \max(\hat{H}), \\ \hat{H}_M &= 1, \end{aligned} \right. \quad (40)$$

at convergence. Now, point M is not known in advance but changes during the SCF iterations. In fact, this ‘self-normalized’ version is obtained by computing constants $C_2$ and $C_3$ at time $t$ and constant $C_1$ at time $t+1$ from the new estimate $\hat{H}^{(t+1)}$. This is feasible because only $C_1$ depends upon $\hat{H}_M$. At step $t$, we successively determine

(i) $\hat{\rho}^{(t)}$ from equation (21 b),

(ii) $\hat{\Psi}^{(t)}$ from equation (21 c),

(iii) $\hat{\Phi}^{(t)}$ from equation (22),

(iv) constants $C_2^{(t)}$ and $C_3^{(t)}$ from equations (37b) and (37c),

(v) $(C_1 \hat{H})^{(t+1)}$ from equation (21a), namely

$$(C_1 \hat{H})^{(t+1)} = C_3^{(t)} - C_2^{(t)} \hat{\Phi}^{(t)} - \hat{\Psi}^{(t)}, \quad (41)$$

(vi) its maximum value, happening at some point M

$$\max[C_3^{(t)} - C_2^{(t)} \hat{\Phi}^{(t)} - \hat{\Psi}^{(t)}] = C_1^{(t+1)}, \quad (42)$$

(vii) and the new enthalpy field

$$\hat{H}^{(t+1)} = \frac{C_3^{(t)} - C_2^{(t)} \hat{\Phi}^{(t)} - \hat{\Psi}^{(t)}}{C_1^{(t+1)}}. \quad (43)$$

Since the computational grid is fixed, $\hat{\Phi}$ does not vary during the cycle and so step (iii) can be removed from the list and executed once for all. Convergence is checked after step (vii). A natural indicator of convergence is (see below)

$$\delta \hat{H}^{(t+1)} = \hat{H}^{(t+1)} - \hat{H}^{(t)}. \quad (44)$$

#### 3.3 Multigrid

The main difficulty arising in the actual problem is the determination of the potential in equation (21a), which varies during the SCF iterations due to the evolution of the enthalpy. We solve equation (21c) for $\hat{\Psi}$ on a uniform grid $\{\hat{R}_i, \hat{Z}_j\}$ from multigrid. This technique is very efficient for solving partial differential equations (PDEs) since it combines nested grids and relaxation (e.g. Briggs, Henson & McCormick 2000; Gheller, Sartoretto & Guidolin 2004; Guillet & Teyssier 2011). It has already been employed in the context of figures of equilibrium (Lanza 1987, 1992; Brandt & Lanza 1988). Short wavelengths are fastly smoothed by relaxation on a given grid, in contrast to long wavelengths. This is circumvented by using low-resolution versions of the finest grid, thereby
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Uniform grid used to solve the Poisson equation from multigrid. At the finest level, it contains $2^\ell + 1$ nodes per direction. BVs are computed for part of the grid where $Z \ge 0$ and then copied below the equatorial plane. There are $N_0 \ge 0$ nodes leftwards/rightwards to the inner/outer edge free of matter.</figcaption></figure>

converting long wavelengths into small ones. The communication between all levels, from the coarsest to the finest grid and vice versa, is achieved from restriction and prolongation operators, enabling us to write fields and solve the discretized PDE at all levels. The restriction operator is based on full weighting, while the prolongation operator employs bilinear interpolation. All the grids are visited, from finest grid to coarsest and back. This is called a 'V-cycle'. At each level, one red-black Gauss-Seidel algorithm (or more) is performed. Despite a good convergence rate (about one extra digit fixed by V-cycle), a unique V-cycle does not deliver the exact solution of the discretized problem, only a better approximation. However, by applying a V-cycle several times, multigrid becomes almost equivalent to matrix diagonalization methods. In single precision, about 7 successive V-cycles are necessary typically, and this is 15 in double precision.

Because multigrid has nominal performance when the spacing is the same in all directions, we use a square grid [$\hat{R}_0$, $\hat{R}_N$] $\times$ [$\hat{Z}_0$, $\hat{Z}_N$]. In order to vary the fluid section-to-grid area ratio (i.e. the covering factor), we put $N_0 \ge 0$ nodes in between the inner and outer edge of the fluid and the grid boundary, as indicated in Fig. 3. Since there are no forces susceptible to thicken the torus in the Z-direction (i.e. the section S is expected to be oblate), the grid size is imposed by the equatorial diameter of the torus $2\hat{R}_e = \hat{R}_B - \hat{R}_A$ (expected to be larger than the polar one $2\hat{R}_p$). In spite of the equatorial symmetry, we found no real advantage to work in the half-plane $Z \ge 0$ only. This would mean (i) Neumann conditions at the grid bottom $Z = 0$, (ii) much more Dirichlet BVs to compute and (iii) a smaller covering factor. The grid is therefore centred in the equatorial plane. Under these conditions, the number of subintervals between points A and B is

$$N' = N - 2N_0, \qquad (45)$$

and the numerical resolution of the finest grid is

$$\hat{h} = \frac{2\hat{R}_e}{N'}, \qquad (46)$$

both in the $\hat{R}$- and $\hat{Z}$-directions. The grid extends radially from $R_0$ to $R_N$ with

$$\begin{cases} \hat{R}_0 = \hat{R}_A - N_0 \hat{h} \ge 0, \\ \hat{R}_N = \hat{R}_B + N_0 \hat{h}, \end{cases} \qquad (47)$$

and we have

$$\begin{cases} \hat{Z}_N = \frac{1}{2}(\hat{R}_N - \hat{R}_0), \\ \hat{Z}_0 = -\hat{Z}_N, \end{cases} \qquad (48)$$

in the vertical direction. The sampling is

$$\begin{cases} \hat{R}_i = \hat{R}_0 + i\hat{h}, & i \in [0, N], \\ \hat{Z}_i = \hat{Z}_0 + j\hat{h}, & j \in [0, N]. \end{cases} \qquad (49)$$

This is for the finest level $\ell$ (i.e. highest resolution) and we have

$$N + 1 = 2^{\ell} + 1, \qquad (50)$$

per direction. The nested grids are easily generated. With mesh size doubling, the coarser grid at level $\ell - 1$ is obtained by keeping nodes $2i$ and $2j$. At level $\ell - 2$, we keep nodes $4i$ and $4j$. At level $1 < k < \ell$, we keep nodes $2^{(\ell-k)}i$ and $2^{(\ell-k)}j$. For the coarsest level, $k = 1$, there are just three nodes per direction, nine in total, where eight belong to the grid boundary.

## 3.4 Discretization of the Poisson equation. Error analysis

The discretization of equation (21 c) writes

$$\left( \delta_{\hat{R}\hat{R}}^2 + \frac{1}{\hat{R}_i} \delta_{\hat{R}} + \delta_{\hat{Z}\hat{Z}}^2 \right) \hat{\Psi}_{i,j} + \tau_{i,j} = 4\pi \hat{\rho}_{i,j}, \qquad (51)$$

where $\delta$ and $\delta^2$ denote the approximations for the first and second derivatives and $\tau = O(\hat{h}^2)$ is the truncation error. By using second-order centred schemes, the dimensionless potential at node i, j is given by

$$\hat{\Psi}_{i,j} = \frac{1}{4} \left( a_+ \hat{\Psi}_{i+1,j} + a_- \hat{\Psi}_{i-1,j} + \hat{\Psi}_{i,j+1} + \hat{\Psi}_{i,j-1} - 4\pi \hat{\rho}_{i,j} \hat{h}^2 \right), \qquad (52)$$

where

$$\begin{cases} a_+ = 1 + \frac{\hat{h}}{2\hat{R}_i}, \\ a_- = 1 - \frac{\hat{h}}{2\hat{R}_i}. \end{cases} \qquad (53)$$

The magnitude of the scheme error is

$$\tau \approx -\frac{\hat{h}^2}{12} \left( \partial_{\hat{R}\hat{R}\hat{R}}^4 + 2 \frac{1}{R} \partial_{\hat{R}\hat{R}\hat{R}}^3 + \partial_{\hat{Z}\hat{Z}\hat{Z}}^4 \right) \hat{\Psi}, \qquad (54)$$

and explicitly decreases with $\hat{h}$. As usual in error analysis (Mathews & Fink 2004), we assume that the third and fourth derivatives are bounded, i.e.

$$|\partial_{\hat{R}\hat{R}\hat{R}}^4 \hat{\Psi}| + 2 \frac{1}{R} |\partial_{\hat{R}\hat{R}\hat{R}}^3 \hat{\Psi}| + |\partial_{\hat{Z}\hat{Z}\hat{Z}}^4 \hat{\Psi}| \le D. \qquad (55)$$

At a node (i, j), the total error is then given by

$$E_{i,j}(\hat{h}) = \hat{\Psi}_{i,j} \left( \frac{4e}{\hat{h}^2} + \frac{e}{\hat{h}\hat{R}} \right) + \frac{\hat{h}^2}{12} D, \qquad (56)$$

where e is the computer precision. As is well known, finite differences become inaccurate if the spacing is too small, leaving round-off errors and loss of significance. The error is minimum for

$$\hat{h}^4 - \frac{6e \hat{\Psi}_{i,j}}{D \hat{R}} \hat{h} - \frac{48e \hat{\Psi}_{i,j}}{D} = 0. \qquad (57)$$
<table><caption>**Table 1.** Optimal spacing and associated multigrid level for second-order accuracy, depending on the computer precision *e*. The last rows give the absolute error expected in potential values for $\ell = 7$.</caption><thead><tr><th></th><th>Single precision</th><th>Double precision</th></tr></thead><tbody><tr><td>e</td><td>6 × 10<sup>-8</sup></td><td>1 × 10<sup>-16</sup></td></tr><tr><td>ħ<sub>opt</sub></td><td>0.006–0.04</td><td>3 × 10<sup>-4</sup>–9 × 10<sup>-6</sup></td></tr><tr><td>ℓ<sub>opt</sub></td><td>5–7</td><td>12–17</td></tr><tr><td>E/Ψ</td><td>0.007–2 × 10<sup>-4</sup></td><td>5 × 10<sup>-9</sup>–5 × 10<sup>-6</sup></td></tr><tr><td>ℓ</td><td></td><td>7</td></tr><tr><td>Absolute error</td><td></td><td>~ħ<sup>2</sup> ~ 2 × 10<sup>-5</sup></td></tr></tbody></table>

When the second derivative is the limiting quantity, the optimal spacing is $\hat{h}_{opt} \sim e^{1/4}$. The first derivative can bring a similar contribution if $\hat{h} \sim R$ typically (i.e. the computational grid is close to the Z-axis). It all depends on $\Psi/D$, which is difficult to estimate. The optimal spacing is in the range

$$ \left( \frac{6e\hat{\Psi}_{i,j}}{D\hat{R}} \right)^{1/3} \lesssim \hat{h}_{opt} \lesssim \left( \frac{48e\hat{\Psi}_{i,j}}{D} \right)^{1/4}, \quad (58) $$

which can be converted into an optimal level $\ell_{opt}$ of the multigrid from equation (50). With standard numbers, we get values listed in Table 1.

### 3.5 Boundary values (BVs)

Under axial symmetry and equatorial symmetries, equation (4) simplifies somehow. The adimensioned potential takes the form (Durand 1953; Huré 2005)

$$ \hat{\Psi} = \iint_S \hat{\rho} \kappa d\hat{z} d\hat{a}, \quad (59) $$

where

$$ \begin{cases} \hat{a} \equiv \frac{a}{L}, \\ \hat{z} \equiv \frac{z}{L}, \end{cases} \quad (60) $$

$$ \kappa = -2\sqrt{\frac{\hat{a}}{\hat{R}}} [k_{+} \mathbf{K}(k_{+}) + k_{-} \mathbf{K}(k_{-})], \quad (61) $$

$$ \mathbf{K}(k) = \int_{0}^{\pi/2} \frac{d\phi}{\sqrt{1 - k^2 \sin^2 \phi}} \quad (62) $$

is the complete elliptic integral of the first kind,

$$ k_{\pm} = \frac{2\sqrt{\hat{a}\hat{R}}}{\sqrt{(\hat{a} + \hat{R})^2 + \hat{\zeta}_{\pm}^2}} \quad (63) $$

is the modulus, $\hat{\zeta}_{+} = \hat{Z} - \hat{z}$ and $\hat{\zeta}_{-} = \hat{Z} + \hat{z}$. When $k \to 1$ (which occurs when $\hat{R} \to \hat{a}$ and $\hat{\zeta} \to 0$), $\kappa$ diverges logarithmically, and the determination of $\Psi$ can become inaccurate. This never happens here since the fluid and the grid boundary are never in contact as soon as $N_0 > 0$ (see Section 3.3).

### 3.6 Note on quadratures

We need a quadrature scheme which preserves the second order of the finite different scheme. Not only BVs are concerned, but also most output quantities (mass, volume, angular momentum, virial parameter, etc.). The trapezoidal rule is the natural choice. For any

function $\hat{f}$ which vanishes outside the fluid (if it depends on $\rho$), we can consider the full grid, i.e.

$$ \frac{1}{\hat{h}^2} \iint_{\text{fluid}} \hat{f} \, d\hat{a} \, d\hat{z} \approx \sum_{i \in [0, N], j \in [0, N]} \hat{f}_{i,j} \\ -\frac{1}{2} \sum_{i \in [0, N]} (\hat{f}_{i,0} + \hat{f}_{i,N}) \\ -\frac{1}{2} \sum_{j \in [0, N]} (\hat{f}_{0,j} + \hat{f}_{N,j}) \\ +\frac{1}{4} (\hat{f}_{0,0} + \hat{f}_{N,0} + \hat{f}_{0,N} + \hat{f}_{N,N}). \quad (64) $$

For geometrical quantities like the fluid volume, one must multiply $\hat{f}$ by a weighting function $w$, namely

$$ w = \begin{cases} 1 & \text{if } \hat{H} \ge 0 \quad (\text{inside the fluid}), \\ 0 & \text{if } \hat{H} < 0 \quad (\text{outside}). \end{cases} \quad (65) $$

There are two pending problems which can damage the second order. First, $\Gamma$ is rounded and does not follow the grid lines (see Fig. 1). Standard quadrature schemes are generally not well suited to irregular or curved integration domains. This is a natural motivation for space mapping (Grandclément et al. 2001; Rieutord, Espinosa Lara & Putigny 2016). Secondly, for hard EOS ($n < 1$), $\hat{p}$ does not reach $\Gamma$ through a finite gradient. Again, this is a situation where many quadrature schemes fail to be accurate. The question of determining the fluid boundary $\Gamma$ is then naturally addressed.

It is instructive to estimate the relative contribution of the numerical cells located just below the fluid surface. Let $\hat{f}_{\Gamma}$ denote the typical value of $\hat{f}$ on to or close to $\Gamma$, $\langle \hat{f} \rangle$ be the mean value of $\hat{f}$ inside the fluid, $\hat{p}$ the perimeter and $\hat{S}$ the fluid section. This contribution is given by

$$ \frac{\hat{f}_{\Gamma} \hat{p} \hat{h}}{\iint_{\text{fluid}} \hat{f} d\hat{a} d\hat{z}} \sim \frac{\hat{f}_{\Gamma} \hat{p}}{\langle \hat{f} \rangle \hat{S}} \hat{h}. \quad (66) $$

It is first order in $\hat{h}$ if $\hat{f}_{\Gamma}/\langle \hat{f} \rangle \sim 1$, but second order if $\hat{f}_{\Gamma}/\langle \hat{f} \rangle \sim \hat{h}$. Further, if $\hat{f}$ depends linearly on $\hat{p}$ (see Section 2.5), then this contribution is $\sim \hat{h}^{1+n}$. It immediately follows that all geometrical quantities like the fluid volume are only first-order accurate, and all quantities which are sensitive to the mass density (BVs, mass, angular momentum, virial kernels, etc.) are second-order accurate as soon as $n \ge 1$ (i.e. for soft EOS). For $n < 1$, the situation is worse (see Section 2.4): not only the fluid boundary is first order, but the gradient of mass density is infinite at $\Gamma$. Increasing the order of the quadrature scheme is not the right way to fix the problem. Instead, a specific treatment must be considered. This is discussed in Section 5.

## 4 PERFORMANCES AT SECOND ORDER

We have built a FORTRAN90 code named DROP for ‘Differentially Rotating Polytropes’ for solving toroidal figures of equilibrium by using methods described above. The main input parameters are (see Section 3.3)

(i) $\ell$ (or $N$) for the multigrid method and grid,

(ii) $n$, $\hat{R}_A$, $\hat{R}_B$ for the polytropic torus,

(iii) $N_0$ for the covering factor.
The starting guess $\hat{H}^{(0)}$ is obviously a critical point. We use a paraboloid with circular section (radius $\hat{R}_e$), which is, in general, a good choice. Regarding the convergence criterion, we use the Frobenius norm $\delta\hat{H}_F$, namely

$$ \delta\hat{H}_F^2 = \frac{1}{(N-1)^2} \sum_{i,j=1,N-1} \delta\hat{H}_{i,j}^2, \qquad (67) $$

which is a good compromise, and so the solution is found when

$$ \delta\hat{H}_F \le \epsilon, \qquad (68) $$

where $\epsilon$ is a small number. It is tempting to link this threshold to the scheme errors. In some cases however, the process seems to converge but finally oscillates without getting stable, which could lead to false solutions. This is observed for instance close to the mass shedding limit. In practice, we therefore let the field $\delta\hat{H}$ fall down at its minimum, which is generally close to the computer precision. This does not appear detrimental.

Once the solution (i.e. $\hat{\rho}$ and the three constants $C_1, C_2$ and $C_3$) is known, several output quantities can be inferred, in particular (in dimensionless form)

(i) the area of the meridional section $\hat{S} = \iint_S d\hat{a}d\hat{z}$,

(ii) the fluid volume $\hat{V} = 2\pi \iint_S d\hat{a}d\hat{z}$,

(iii) the mass $\hat{M} = 2\pi \iint_S \hat{\rho}d\hat{a}d\hat{z}$,

(iv) the mean mass density $\langle \hat{\rho} \rangle = \frac{\hat{M}}{\hat{V}}$,

(v) the angular momentum $\hat{J} = 2\pi \iint_S \hat{\rho}\hat{\Omega}\hat{a}^3d\hat{a}d\hat{z}$,

(vi) $\hat{W}, \hat{\Pi}, \hat{T}$ and the virial parameter VP.

## 4.1 An example with $n \ge 1$

As a first example, we consider a rigidly rotating torus with parameters $n = 1.5$, $\hat{R}_A = 0.5$, $\hat{R}_B = 1$. This case is for instance reported in Hachisu (1986) and others. For $\Omega = \Omega_0$, the centrifugal potential is quadratic with the radius, namely

$$ \begin{cases} \hat{\Omega} = 1, \\ \hat{\Phi} = -\frac{1}{2}\hat{R}^2. \end{cases} \qquad (69) $$

Regarding the numerical grid, we take $\ell = 7$ (i.e. $N = 128$) and $N_0 = 2$, which implies for the radial sampling $\hat{R}_2 = \hat{R}_A$ and $\hat{R}_{N-2} = \hat{R}_B$. According to Section 3.3, we have

$$ \left\{ \begin{aligned} \hat{h} &= 0.004\,032\,26\dots \\ \hat{R}_0 &= 0.491\,935\dots \\ \hat{R}_N &= 1.008\,065\dots \\ \hat{Z}_N &= \frac{1}{2}(\hat{R}_N - \hat{R}_0) = +0.258\,065\dots \end{aligned} \right. \qquad (70) $$

For these parameters, the absolute error in potential values is expected to be of the order of $\sim 2 \times 10^{-5}$ (see Section 3.4), which is transmitted to $\hat{H}$, $\hat{\rho}$, constants $C_1, C_3$, VP, etc. This means four correct digits typically, which is typical at such a spatial resolution (Eriguchi 1978; Hachisu 1986; Yoshida & Eriguchi 2006). The equilibrium solution is found after 65 iterations. The remarkable capability of convergence of the SCF loop is shown in Fig. 4, which displays $\delta\hat{H}$ a few steps after. Fig. 5 displays $\hat{\Psi}$, $\hat{H}$ and $\hat{\rho}$ at convergence. We notice the slight asymmetry of the enthalpy and density maps (point M is closer to edge A than edge B), due to curvature effects. Centrifugation shifts the gravitational potential well towards the inner edge. Equatorial forces are plotted in Fig. 6. We see that the zero-gravity and zero-pressure points do not coincide.

<figure><img src="image_6.png" /><figcaption>Figure 4. Deviation $\delta\hat{H}$ at convergence for the torus considered in Section 4.1, i.e. $n = 1.5$, $\hat{R}_A = 0.5$, $\hat{R}_B = 1$, $\ell = 7$ and $N_0 = 2$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5. Same legend as for Fig. 4 but for the potential $\hat{\Psi}$, the enthalpy $\hat{H}$ and the mass density $\hat{\rho}$ at equilibrium (from top to bottom).</figcaption></figure>

The virial kernel and the quantity $\frac{\kappa U}{2\kappa_T + \kappa_U}$ are displayed in Fig. 7. The gravitational energy density is more concentrated at the inner edge of the torus, while the internal and centrifugal energies dominate in the outer region. The kinetic term always exceeds the pressure term, and this is especially pronounced close to the fluid boundary. Table 2 shows that output quantities are in agreement with Hachisu (1986) who works with Simpson's fourth-order quadrature rule in spherical coordinates and rather low covering factor. Successful comparisons
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Forces at the equatorial plane versus the radius for the torus considered in Fig. 5.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Same legend as for Fig. 5 but for the virial kernel $\kappa_{Virial}$ (top) and the ratio $\frac{\kappa_{f}}{2\kappa_{T}+\kappa_{U}}$ (bottom).</figcaption></figure>

are observed for many other configurations. Deviations appear close to critical rotation and in the thin ring limit as well.

## 4.2 VP versus covering factor

By varying $N_0$, we modify the covering factor (fluid section-to-grid area ratio) defined by

$$ \frac{\hat{S}}{(\hat{R}_N - \hat{R}_0)(\hat{Z}_N - \hat{Z}_0)} \equiv \Lambda. \qquad (71) $$

In the above example, $\Lambda \approx 0.61$. When increasing $N_0$ (from 0 to $N'/3$), the resolution of the fluid is decreased, which decreases $\Lambda$. Fig. 8 displays VP versus $\Lambda$ in log scale. While the sensitivity remains weak, there is a net linear correlation, and we find $\log \text{VP} \approx -4.58 - 0.99\log \Lambda$. The more resolved the fluid, the

<table><caption>**Table 2.** Comparison with Hachisu (1986). At the actual resolution, only four figures are guaranteed ($\ast$estimated).</caption><thead><tr><th>Quantity</th><th>Hachisu (1986)</th><th>This work</th></tr></thead><tbody><tr><td>Covering factor Λ</td><td>0.091*</td><td>0.608</td></tr><tr><td>C<sub>1</sub></td><td>0.0842</td><td>0.08422</td></tr><tr><td>C<sub>2</sub> (i.e. Ω<sub>0</sub><sup>2</sup>)</td><td>0.207</td><td>0.2068</td></tr><tr><td>-C<sub>3</sub></td><td></td><td>0.3691</td></tr><tr><td>R&#x0302;<sub>e</sub></td><td>0.5</td><td>0.5</td></tr><tr><td>R&#x0302;<sub>p</sub></td><td></td><td>0.2071</td></tr><tr><td>S&#x0302;</td><td></td><td>0.1619</td></tr><tr><td>V&#x0302;</td><td>0.762</td><td>0.7615</td></tr><tr><td>M&#x0302;</td><td>0.219</td><td>0.2188</td></tr><tr><td>⟨ρ⟩</td><td>0.287*</td><td>0.2874</td></tr><tr><td>Max. enthalpy</td><td>0.0842</td><td>0.08422</td></tr><tr><td>Max. pressure</td><td>0.0337</td><td>0.03369</td></tr><tr><td>Max. density</td><td>0.0245*</td><td>0.02444</td></tr><tr><td>√C<sub>2</sub>J&#x0302;</td><td>0.0562</td><td>0.05613</td></tr><tr><td>C<sub>2</sub>T&#x0302;</td><td>0.0128</td><td>0.01276</td></tr><tr><td>-W&#x0302;</td><td>0.0401</td><td>0.04005</td></tr><tr><td>C<sub>1</sub>/n<sub>1</sub>U&#x0302;</td><td>0.0145</td><td>0.01453</td></tr><tr><td>log (VP)</td><td>(?)</td><td>-4.42</td></tr><tr><td>Iterations</td><td>?</td><td>65</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Variation of the virial parameter with the covering factor $\Lambda$. Values of $N_0$ are reported (see Fig. 5 for $N_0 = 2$ and Fig. 9 for $N_0 = 42$).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Same legend as for Fig. 5(b) but with $N_0 = 42$.</figcaption></figure>

better the virial parameter. Fig. 9 shows the enthalpy at equilibrium for the largest value $N_0 = 42$.

## 4.3 VP as a reliable indicator for precision

We have considered various numerical resolutions $\hat{h}$ for the finest grid, which is obtained by varying $\ell$, according to equation (46). We
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Variation of the virial parameter VP with the resolution number $N'$ for $\ell \in [3, 11]$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Precision observed on a few output quantities versus the resolution number $N'$. The volume is not precisely determined unless the fluid boundary is accounted for (see Section 5).</figcaption></figure>

have explored the range $\ell \in [3, 11]$. The convergence of the SCF loop is almost unchanged, with 65 iterations, whatever $\ell$. The virial parameter is plotted versus $N'$ in Fig. 10 for both single and double precision. The correlation is remarkable. A linear regression yields $\log \text{VP} \approx -0.16915 - 2.02367 \log N'$, or

$$ \text{VP} \approx 0.339 \times \hat{h}^{2.02}, \qquad (72) $$

which means that the virial parameter is second-order accurate. For single precision, accuracy is maximum for $\ell \sim 7$ and deteriorates for finer grids, in agreement with the error analysis (see Section 3.4).

We can easily check that all output quantities are improved when the resolution is increased by comparing the results at $\ell - 1$ and $\ell$, namely

$$ \left| 1 - \frac{\hat{f}(\ell - 1)}{\hat{f}(\ell)} \right| \sim \left( \frac{d \ln |\hat{f}|}{d \ell} \right)_{\ell-1}. \qquad (73) $$

Fig. 11 shows this quantity for the mass, the two constants $C_2$ and $C_3$, the angular momentum $\hat{J}$ and the volume $\hat{V}$. We see that the virial parameter is a good tracer of precision for most output quantities, within a factor of less than 2 typically. The exception

<figure><img src="image_14.png" /><figcaption>**Figure 12.** The $\omega^2-j^2$ diagram for $n = 1.5$ at two different scales. Near critical rotation (bottom panel), the sequence is not well defined unless the fluid boundary is detected and accounted for. The value given in Hachisu (1986) stands outside the plot.</figcaption></figure>

is for the volume, and for the meridional section area as well (not shown). This is in agreement with the discussion in Section 3.6.

## 4.4 Uncertainty of the $\omega^2-j^2$ diagram near critical rotation

A classical manner to visualize and compare equilibria obtained for different triplets ($C_1, C_2, C_3$) is the $\omega^2-j^2$ diagram (Chandrasekhar 1969; Hachisu 1986), where

$$ \left\{ \begin{aligned} j^2 &= \frac{1}{4\pi G(\rho) M^2 V^{4/3}}, \\ \omega^2 &= \frac{1}{4\pi G(\rho)} \Omega_0^2, \end{aligned} \right. \qquad (74) $$

and $J = \rho_0 \Omega_0 L^5 \hat{J}$. Using dimensionless quantities, we have

$$ (j^2, \omega^2) = \frac{C_2}{4\pi \hat{M}} \times \left( \frac{\hat{J}^2}{\hat{M}^2 \hat{V}^{1/3}}, \hat{V} \right). \qquad (75) $$

For one-ring sequences obtained by varying the torus axis ratio (while $n$ is held fixed), $\omega$ typically decreases when $j$ increases (Eriguchi & Sugimoto 1981; Hachisu 1986; Ansorg et al. 2003). Fig. 12 displays the diagram obtained by decreasing $\hat{R}_A/\hat{R}_B$. We notice that the ‘trajectory’ is hazardous close to the sequence end. This effect is a direct consequence of the uncertainty in the volume. Critical rotation, where the pressure gradient vanishes at point B (i.e. gravity and centrifugation are just balanced), is found for $\hat{R}_A/\hat{R}_B \sim 0.3238$, against 0.325 in Hachisu (1986).
<table><caption>**Table 3.** Results obtained for $n = 1$ and $\hat{R}_A/\hat{R}_B = 0.9$ compared with the spectral approach by Petroff & Horatschek (2008). Same grid parameters as for Table 2.</caption><thead><tr><th>Quantity</th><th>Petroff & Horatschek (2008)</th><th>This work</th></tr></thead><tbody><tr><td>Covering factor</td><td>1</td><td>0.731</td></tr><tr><td>C<sub>1</sub></td><td></td><td>0.005 363</td></tr><tr><td>C<sub>2</sub> (i.e. Ω<sup>2</sup><sub>0</sub>)</td><td>0.014 99</td><td>0.014 99</td></tr><tr><td>-C<sub>3</sub></td><td></td><td>0.040 55</td></tr><tr><td>i<sub>M</sub></td><td></td><td>63</td></tr><tr><td>ħ<sub>M</sub></td><td></td><td>1</td></tr><tr><td>ħ</td><td></td><td>0.007 786</td></tr><tr><td>ħ̂</td><td></td><td>0.046 47</td></tr><tr><td>ħ̃</td><td></td><td>0.002 004</td></tr><tr><td>⟨ρ⟩</td><td></td><td>0.4312</td></tr><tr><td>Max. enthalpy</td><td></td><td>0.005 363</td></tr><tr><td>Max. pressure</td><td></td><td>0.002 682</td></tr><tr><td>Max. density</td><td></td><td>0.005 363</td></tr><tr><td>√C<sub>2</sub>ħ</td><td></td><td>0.002 216</td></tr><tr><td>C<sub>2</sub>ħ</td><td></td><td>0.000 1357</td></tr><tr><td>-ħ</td><td></td><td>0.000 3720</td></tr><tr><td>C<sub>1</sub>/n+1 ħ̃</td><td></td><td>0.000 1006</td></tr><tr><td>log (VP)</td><td></td><td>-4.75</td></tr><tr><td>Iterations</td><td></td><td>28</td></tr><tr><td>ħ̃ (n+1/C<sub>1</sub>)<sup>3/2</sup></td><td>144.3</td><td>144.3</td></tr><tr><td>ħ̃ (n+1/C<sub>1</sub>)<sup>5/2</sup></td><td>5953</td><td>5952</td></tr><tr><td>ħ̃ (n+1/C<sub>1</sub>)<sup>5/2</sup></td><td>364.5</td><td>364.4</td></tr><tr><td>-ħ̃ (n+1/C<sub>1</sub>)<sup>5/2</sup></td><td>999.2</td><td>999.0</td></tr><tr><td>ħ̃ (n+1/C<sub>1</sub>)<sup>5/2</sup></td><td>90.07</td><td>90.06</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Enthalpy gradient at point B near critical rotation for $n = 1$ and two numerical resolutions.</figcaption></figure>

### 4.5 A case with $n = 1$

As a second example, we consider a rigidly rotating torus with $n = 1$ and $\hat{R}_A/\hat{R}_B = 0.9$. The solution is found after 28 SCF iterations. The internal structure resembles very much the one presented above. The torus boundary is much closer to a circle, with a polar radius $\hat{R}_p = 0.0495$. Output quantities are given in Table 3. The agreement with the results obtained from a spectral method by Petroff & Horatschek (2008) is remarkable. From this equilibrium, we can reach the critical rotation by decreasing the radius of the inner edge (the number of SCF iterations generally increases). Fig. 13

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Same legend as for Figs 5 and 7 but for $n = 1$ at critical rotation where $\hat{R}_A/\hat{R}_B \approx 0.253 221$.</figcaption></figure>

shows $\nabla_R \hat{H}$ at the outer edge versus $\hat{R}_A$. There is obviously a slight resolution effect, but the gradient vanishes for $\hat{R}_A/\hat{R}_B \sim 0.2532$, which is also the solution of Petroff & Horatschek (2008) obtained from a spectral approach. Fig. 14 displays the mass density, virial kernel and fraction $\frac{\kappa_{\text{J}}}{2\kappa_{\text{T}}+\kappa_{\text{J}}}$. The fluid core is dominated by the gravitational energy, and it remains more marked close to the inner edge. We see the elongated structure at the outer edge where matter is in pure Keplerian rotation.

### 4.6 Case with a hard EOS

We consider the same parameters as in Section 4.1 but with $n = 0.5$. This is a typical case where $\rho$ has an infinite gradient on to $\Gamma$. The equilibrium solution is now found after 34 iterations. In general, convergence is faster with low indices $n$. The mass density, virial kernel and ratio $\frac{\kappa_{\text{J}}}{2\kappa_{\text{T}}+\kappa_{\text{J}}}$ are displayed in Fig. 15. Again, the comparison with Hachisu (1986) does not reveal any deviations for the first three digits. The virial test indicates $\log(\text{VP}) \approx -4.2$. This is astonishing since, as argued in Section 3.6, quadratures are expected to be inaccurate for such polytropic index. This illustrates the fact that the virial test does not prove that the actual solution is the good one, but that it is numerically clean and self-consistent. We have run the code for various resolution numbers where $\ell \in [3, 11]$ and computed the variation of some
<figure><img src="image_17.png" /><figcaption>Figure 15. Same legend as for Fig. 5 but for $n = 0.5$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 16. Same legend as for Fig. 11 but for $n = 0.5$.</figcaption></figure>

quantities, as done above. The results are displayed in Fig. 16 and are to be compared with Fig. 11. While the global trend is still there, the order is oscillating. All quantities are concerned, not only geometrical ones.

<figure><img src="image_19.png" /><figcaption>Figure 17. Double integrals over the fluid section are improved by dividing the fluid into three subdomains. Close to the inner and outer edges, the radial quadrature is performed first. In the c-subdomain, the vertical quadrature is performed first.</figcaption></figure>

# 5 SUBGRID APPROACH: METHOD AND EXAMPLES

## 5.1 Detection of the fluid boundary $\Gamma$. Case of a soft EOS

Detecting the fluid boundary is basically achieved by scanning the whole grid radially or/and vertically and searching for nodes where the sign of enthalpy changes from one node to the other. There are degeneracies if $\Gamma$ is not convex, which case is met for instance when matter is pinched at one or more places around the equatorial plane (Ansorg et al. 2003). We then use a more straightforward method based on the Freeman chain code (e.g. Freeman 1961): $\Gamma$ is gradually revealed by visiting neighbouring cells following the same searching sequence. It is in fact much faster than by direct scan because the grid is not explored as a whole. The algorithm starts at point A and fastly reaches point B, which generates $N_\Gamma$ points with coordinates ($\hat{R}_\Gamma$, $\hat{Z}_\Gamma$) where $\hat{H} = 0$ (the root finding method is second-order accurate here). If $\Gamma$ has length $\hat{p}_e$, then $N_\Gamma \sim \hat{p}_e/\hat{h}$ typically.

As a matter of fact, we do not need all the $\Gamma$-nodes so detected to improve the quality of double integrals. Basically, the inner and outer edges where $d\hat{Z}_\Gamma/d\hat{R}$ is large must be oversampled. A good compromise consists in dividing the fluid section into three subdomains:

(i) an 'l-subdomain' for the inner edge,

(ii) an 'r-subdomain' for the outer edge,

(iii) a 'c-subdomain' located in between,

which requires two intermediate radii $\hat{R}_l$ and $\hat{R}_r$, as depicted in Fig. 17. For the l- and r-subdomains, we integrate in the radial direction first by considering the boundary nodes, then in the vertical direction. For the c-subdomain, we proceed in the reverse way. The double integral is therefore written in the form

$$ \frac{1}{2} \iint_{\text{fluid}} \hat{f} \, d\hat{a} \, d\hat{z} = \underbrace{\int_{0}^{\hat{Z}_{\Gamma}} d\hat{z} \int_{\hat{R}_{\Gamma}(\hat{z})}^{\hat{R}_{l}(\hat{z})} \hat{f} \, d\hat{a}}_{\text{l-subdomain}} + \underbrace{\int_{0}^{\hat{R}_{l}} d\hat{a} \int_{\hat{R}_{l}(\hat{z})}^{\hat{Z}_{\Gamma}(\hat{a})} \hat{f} \, d\hat{z}}_{\text{c-subdomain}} + \underbrace{\int_{0}^{\hat{Z}_{\Gamma}} d\hat{z} \int_{\hat{R}_{r}}^{\hat{R}_{\Gamma}(\hat{z})} \hat{f} \, d\hat{a}}_{\text{r-subdomain}} \quad (76) $$

and executed as is (the factor 2 accounts for matter below the equatorial plane). The quadrature scheme is almost the same as in equation (64), except at the fluid boundary where two nodes have a spacing smaller than $\hat{h}$. This technique is quite easy to implement. The two
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Fluid boundary for the torus with soft EOS considered in Section 4.1. By performing a vertical scan, the precise curvature of $\Gamma$ close to the inner and outer edges is not accessible. Similarly, the radial scan does not give much detail of the boundary at the top/bottom of the fluid.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Same legend as for Fig. 11 but when the fluid volume is computed at second order by accounting for the boundary (red).</figcaption></figure>

radii $\hat{R}_l$ and $\hat{R}_r$ could change inside the SCF iterations (depending on $d\hat{Z}_\Gamma/d\hat{R}$). For the present study, however, we simply take

$$ \begin{cases} \hat{R}_l = \hat{R}_A + \frac{1}{2}\hat{R}_e, \\ \hat{R}_r = \hat{R}_B - \frac{1}{2}\hat{R}_e. \end{cases} \qquad (77) $$

## 5.2 Results for $n = 1.5$

We go back to the example considered in Section 4.1 and recompute the solution including the detection of $\Gamma$, for different numerical resolutions. BVs, which feed the Poisson solver, are then estimated from equation (76) at each time step of the SCF loop. At convergence, output quantities (viral parameter, fluid volume, etc.) are computed in a similar manner. Fig. 18 shows the result obtained for the example given in Section 4.1. We find $N_\Gamma = 227$ for $\ell = 7$. Fig. 19 shows the new volume, denoted by $\hat{V}^*$ versus $N'$. We clearly see that error is now fully second-order accurate for the volume (and other geometrical quantities). As a consequence, the equilibrium sequence near critical rotation is now much better described as Fig. 20 shows.

<figure><img src="image_22.png" /><figcaption>**Figure 20.** Same conditions as for Fig. 12 but when the fluid boundary is detected and accounted for.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 21.** Poisson kernel $\kappa\rho$ limited to the half-plane $Z \ge 0$ versus $a$ and $z$ for the BV at $(\hat{R}_N, 0)$ for the torus with hard EOS considered in Section 4.6.</figcaption></figure>

## 5.3 Detection of the fluid boundary combined with kernel splitting. Case of a hard EOS

If we apply the preceding recipe to configurations with $n \le 1$, the second order is still not recovered for the reason evoked in Section 3.6: the presence of infinite gradients all along the fluid boundary, as a direct consequence of the polytopic assumption, makes the numerical integration too uncertain. This is illustrated in Fig. 21, which displays the full integrand $\kappa\hat{\rho}$ in equation (59) obtained when computing one of the BVs for $n = 0.5$ at convergence. Infinite gradients are clearly visible all along the fluid boundary. The problem is not fixed by increasing the order of the quadrature scheme. A convenient approach consists in treating the underlying power law by analytical means. Actually, by using the generic $q$-coordinate (which is $a$ or $z$ depending on the subdomain; see below), we have

$$ \hat{f}(q)\hat{\rho}(q) = \hat{f}_0\hat{\rho}_0 \left( \frac{q - q_\Gamma}{q_0 - q_\Gamma} \right)^n + \delta(\hat{f}\hat{\rho}), \qquad (78) $$

where $q_0 > q_\Gamma$ is a pivot, $\hat{f}_0\hat{\rho}_0 = \hat{f}(q_0)\hat{\rho}(q_0)$, and $\delta(\hat{f}\hat{\rho})$ is the residue, i.e. the deviation to the power law. By integrating over $q$,
<figure><img src="image_24.png" /><figcaption>**Figure 22.** Same legend as for Fig. 21 but for the residue $\delta(\hat{f}\hat{\rho})$ (i.e. after power-law splitting of the kernel).</figcaption></figure>

we get

$$ 
\begin{aligned}
\int_{q_\Gamma}^q \hat{f}(q')\hat{\rho}(q')dq' &= \int_{q_\Gamma}^q \hat{f}_0\hat{\rho}_0 \left(\frac{q-q_\Gamma}{q_0-q_\Gamma}\right)^n dq' + \int_{q_\Gamma}^q \delta(\hat{f}\hat{\rho})dq' \\
&= \hat{f}_0\hat{\rho}_0 \frac{(q-q_\Gamma)^{n+1}}{(1+n)(q_0-q_\Gamma)^n} + \int_{q_\Gamma}^q \delta(\hat{f}\hat{\rho})dq' \\
&= \hat{f}_0 \frac{\hat{\rho}(q)(q-q_\Gamma)}{(1+n)} + \int_{q_\Gamma}^q \delta(\hat{f}\hat{\rho})dq',
\end{aligned}
\quad (79) $$

where the first analytical term is precisely the contribution of the infinite gradient. The second term has finite gradient on to $\Gamma$ by construction and so, it can be estimated by numerical means without problems at second order from the trapezoidal rule. The kernel splitting method is general and applies for computing both BVs and output values. It has no special efficiency for $n \ge 1$. It requires the knowledge of the fluid boundary.

For the l- and r-subdomains where the radial integration is performed first, we have $q = a$, while this is $q = z$ in the c-subdomain. Note that the pivot point serves only to determine the amplitude of the power-law term through $f_0$ and must not be used in the quadrature scheme because $\delta(\hat{f}\hat{\rho}) = 0$ for both $q = q_\Gamma$ and $q = q_0$. An example of residue $\delta(\hat{f}\hat{\rho})$ is shown in Fig. 22. Infinite derivatives on $\Gamma$ have disappeared.

## 5.4 Results for $n = 0.5$

We have run the code for various resolutions and the results are given in Fig. 23, which is to be compared to Fig. 16. We see that the second-order accuracy is now well reproduced by using kernel splitting, for all quantities, in particular for the virial parameter and fluid volume as well. Table 4 lists output quantities and compares the different approaches. While the boundary $\Gamma$ becomes explicitly an unknown of the problem, which directly impacts on BVs inside the SCF iterations, the convergence process is not fragilized. We can still maintain a very low level of convergence close to the computer precision (see below) and still use inequality (68) as the pertinent convergence criterion. The computing time is substantially increased with kernel splitting, as expected but this is not critical.

<figure><img src="image_25.png" /><figcaption>**Figure 23.** Same legend as for Fig. 16 but with power-law splitting of the kernel.</figcaption></figure>

<table><caption>**Table 4.** Comparison with Hachisu (1986) for $n = 0.5$ (*estimated). Results obtained by accounting for the fluid boundary $\Gamma$ are given in the last two columns.</caption><thead><tr><th rowspan="2">Quantity</th><th rowspan="2">Hachisu (1986)</th><th colspan="2">This work</th></tr><tr><th>Subgrid approach<br>with $\Gamma$</th><th>splitting</th></tr></thead><tbody><tr><td>$\Lambda$</td><td>0.088*</td><td>0.594</td><td>0.594</td></tr><tr><td>$C_1$</td><td>?</td><td>0.1200</td><td>0.1201</td></tr><tr><td>$C_2$ (i.e. $\Omega_0^2$)</td><td>0.383</td><td>0.3825</td><td>0.3823</td></tr><tr><td>$-C_3$</td><td></td><td>0.7729</td><td>0.7725</td></tr><tr><td>$i_M$</td><td></td><td>61</td><td>61</td></tr><tr><td>$\hat{H}_M$</td><td></td><td>1</td><td>1</td></tr><tr><td>$\hat{\delta}$</td><td></td><td>0.1581</td><td>0.1583</td></tr><tr><td>$\hat{V}$</td><td>0.746</td><td>0.7451</td><td>0.7459</td></tr><tr><td>$\hat{M}$</td><td>0.471</td><td>0.4710</td><td>0.4707</td></tr><tr><td>$(\hat{\rho})$</td><td>0.631*</td><td>0.6321</td><td>0.6311</td></tr><tr><td>$\sqrt{C_2}\hat{J}$</td><td>0.173</td><td>0.1732</td><td>0.1730</td></tr><tr><td>$C_2\hat{T}$</td><td>0.536</td><td>0.5354</td><td>0.5349</td></tr><tr><td>$-\hat{W}$</td><td>0.171</td><td>0.1713</td><td>0.1712</td></tr><tr><td>$\frac{C_1}{n+1}\hat{U}$</td><td>0.0642</td><td>0.6423</td><td>0.6429</td></tr><tr><td>log (VP)</td><td>?</td><td>-4.23</td><td>-4.23</td></tr><tr><td>Iterations</td><td>?</td><td>34</td><td>33</td></tr></tbody></table>

# 6 ACCELERATION OF THE SCF ITERATIONS FROM INTERNAL PRE-CONDITIONING

## 6.1 A recipe

The evolution of $\delta\hat{H}_F$ during the SCF iterations is shown in Fig. 24 for the three configurations studied before (see Section 3.1). We see the remarkable stabilization of the SCF iteration at a very high level of precision (see also Fig. 4). As mentioned, this holds whatever the polytropic index, with or without subgrid approaches. The striking feature is the exponential decay which appears after a very few iterations, once the starting guess is almost forgotten. There is obviously a sensitivity to the polytropic index and axial ratio $\hat{R}_A/\hat{R}_B$ of the torus, but this property appears relatively universal. The decay is not only observed for $\delta\hat{H}_F$, but for the enthalpy field at almost all of the grid nodes. It has an origin in the SCF method. We can therefore try to reduce the number of iterations,
<figure><img src="image_26.png" /><figcaption>Figure 24. Variation of $\delta \hat{H}_F$ during the SCF iterations (black) down to convergence for the three configurations discussed in Section 4. Convergence is not impacted by kernel splitting (red dots).</figcaption></figure>

at least empirically, from this decay. Let us assume that, at time $t$, we have

$$ \ln |\delta \hat{H}^{(t)}| \sim a - bt, \qquad (80) $$

where $a > 0$ and $b > 0$. Since this expression also writes

$$ \delta \hat{H}(t + \tau) = \delta \hat{H}^{(t+\tau-1)} e^{-b}, \qquad (81) $$

we can estimate field values at the end of the next iteration, i.e. at time $t+2$. With $\tau=2$, we have

$$ \hat{H}^{(t+1)} + \delta \hat{H}^{(t)} e^{-b} \equiv \hat{H}^{(t+2), \text{estimate}}. \qquad (82) $$

Like $\hat{H}$, $b$ is a function of space (each node has its own path to convergence), namely

$$ b \sim - \left( \frac{d \ln |\delta \hat{H}|}{dt} \right)_{t+1}, \qquad (83) $$

which can be determined from finite differences. At second order, we have

$$ \begin{aligned} -b &\approx \frac{\delta \hat{H}^{(t+\frac{3}{2})} - \delta \hat{H}^{(t+\frac{1}{2})}}{\delta \hat{H}^{(t+1)}} \\ &\approx \frac{\hat{H}^{(t+1)} - 2\hat{H}^{(t)} + \hat{H}^{(t-1)}}{\frac{1}{2} [\hat{H}^{(t+1)} + \hat{H}^{(t-1)}]} \end{aligned} \qquad (84) $$

The mass density field to be considered in BVs when entering step $t+1$ of the SCF method is therefore given by

$$ \hat{H}_{i,j}^{(t+1)} = \sup [\hat{H}^{(t+2), \text{estimate}}, 0]^n, \qquad (85) $$

which requires the storage of the enthalpy over a three successive iterations, at $t-1$, $t$ and $t+1$.

The next question concerns the application of such a correction: where and when can we anticipate and provide an improved guess for the mass density at the beginning of iteration $t+1$? We must first recall that we are mainly interested in BVs, and so only nodes of positive enthalpy are concerned. Secondly, the correction must not perturb the SCF iterations too much, and cancel convergence. There is no risk at nodes where $b > 0$, but this may be critical with $b < 0$. A possible criterion for pre-conditioning is

$$ \left\{ \begin{array}{l} \hat{H}_{i,j} > 0 \\ \delta \hat{H}_{i,j}^{(t+1)} \times \delta \hat{H}_{i,j}^{(t)} > 0, \\ b_{i,j} > b_0 \gtrsim 0, \end{array} \right. \qquad (86) $$

<figure><img src="image_27.png" /><figcaption>Figure 25. Function $b$ computed after five iterations of the classical SCF loop for the torus considered in Section 4.1. The criterion for pre-conditioning is not satisfied everywhere (grey zones). At nodes where $b > 0$ (red zones), the enthalpy field has already started to converge and would be helped by pre-conditioning.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 26. Same legend as for Fig. 24 but when pre-conditioning the mass density inside the cycle (red). Convergence is then attained faster, by a factor of 2 typically. This holds with a subgrid approach (filled circles).</figcaption></figure>

where the second inequality means that the enthalpy field must vary monotonically at the actual node. An accelerated version of the SCF iterations can be obtained by modifying the algorithm given in Section 3.2 as follows.

(viii) Pre-conditioning (optional): if inequalities (86) are fulfilled, then $\rho$ is given by equation (85) when entering step $t+1$, i.e. $\hat{H}^{(t+1)}$ is replaced by $\hat{H}^{(t+2), \text{estimate}}$.

## 6.2 Results

Fig. 25 displays $b$ computed from equation (84) for the torus with $n=1.5$ after five steps of the classical SCF iterations. With a mean value $\langle b \rangle \sim 0.4$, the enthalpy is on the way to convergence as the deviation from one step to the other already diminishes, and pre-conditioning helps to accelerate the process. We have checked this recipe to many configurations and observed mostly no failure indeed. It works well with or without including boundary effects and kernel splitting for $n < 1$. The results are given in Fig. 26 for the three examples discussed in Section 4 where we have used $b_0 = -1$. By using the same convergence criterion, the SCF iterations require about half the normal time.
## 7 CONCLUSION

In this paper, we have reinvestigated toroidal figures of equilibrium with special focus on the convergence process of the SCF iterations and accuracy of solutions. In this context, we have carried out the full problem at second order and discussed the impact of discretization and quadrature schemes on output quantities. As argued, the polytropic hypothesis sets severe constraints on the determination of solution for $n < 1$ due to the non-derivability of the mass density at the fluid boundary. Though a few selected examples, we have shown how the gravitational, centrifugal and internal energy densities are distributed in space. We have demonstrated that the virial parameter is a good tracer of the precision of output quantities, except for geometrical quantities like the volume which fluctuates with the resolution. Such an uncertainty directly impacts on the location of critical rotation in the $\omega^2-j^2$ equilibrium diagrams. Rather than increasing the order of the quadrature scheme, we have shown that the second order can finally be recovered if the fluid boundary is detected and accounted for correctly. For polytropic indices lower than unity, quadrature troubles due to infinite gradients of the mass density can be circumvented through kernel splitting.

Finally, we have shown that it is possible to speed up the convergence of the SCF method by internal pre-conditioning of the enthalpy field. Actually, the temporal evolution of the enthalpy at the grid nodes shows an almost universal exponential decay, which can be used back to improve the new guess inside the SCF loop. This improvement is particularly interesting when building high-resolution equilibrium sequences. It is a step towards the understanding of the convergence of the SCF algorithm, which remains poorly understood and documented (Odrzywołek 2003). The recipe is not limited to toroidal structures and can be probably applied to spheroidal ones, which would be interesting to check.

## ACKNOWLEDGEMENTS

It is a pleasure to thank B. Boutin-Basillais. We would like to dedicate this paper to our colleague and friend Dr J. P. Zahn.

## REFERENCES

* Abramowicz M. A., Karas V., Lanza A., 1998, A&A, 331, 1143

* Amendt P., Lanza A., Abramowicz M. A., 1989, ApJ, 343, 437

* Ansorg M., Kleinwächter A., Meinel R., 2003, MNRAS, 339, 515

* Balmforth N. J., Meacham S. P., Spiegel E. A., Young W. R., 1992, Ann. New York Acad. Sci., 675, 53

* Barker K., 2001, Astron. Geophys., 42, 21

* Beech M., 1988, Ap&SS, 146, 299

* Binney J., Tremaine S., 1987, Galactic Dynamics. Princeton Univ. Press, Princeton, NJ

* Bodenheimer P., Ostriker J. P., 1973, ApJ, 180, 159

* Brandt A., Lanza A., 1988, Class. Quantum Grav., 5, 713

* Briggs W. L., Henson V. E., McCormick S. F., 2000, A Multigrid Tutorial, 2nd edn. SIAM, Philadelphia, PA

* Cai M. J., Shu F. H., 2003, ApJ, 583, 391

* Caimmi R., 1986, A&A, 159, 147

* Chandrasekhar S., 1969, The Silliman Foundation Lectures. Yale Univ. Press, New Haven

* Chandrasekhar S., Lebovitz N. R., 1962, ApJ, 135, 238

* Clement M. J., 1967, ApJ, 148, 159

* Cox J. P., 1968, Principles of Stellar Structure – Vol. 1: Physical Principles; Vol. 2: Applications to Stars

* Curry C. L., McKee C. F., 2000, ApJ, 528, 734

* Durand E., 1953, Electrostatique. Vol. 1: Les Distributions. Masson, Paris

* Dutrey A. et al., 2014, Nature, 514, 600

* Dyson F. W., 1893, Phil. Trans. R. Soc. A, 184, 1041

* Elvis M., 2012, J. Phys.: Conf. Ser., 372, 012032

* Eriguchi Y., 1978, PASJ, 30, 507

* Eriguchi Y., Mueller E., 1985, A&A, 147, 161

* Eriguchi Y., Mueller E., 1993, ApJ, 416, 666

* Eriguchi Y., Sugimoto D., 1981, Prog. Theor. Phys., 65, 1870

* Fang L.-Z., Xiang S.-P., Wang R.-C., 1983, Acta Astrophys. Sin., 3, 55

* Font J. A., Montero P. J., Shibata M., 2010, J. Phys.: Conf. Ser., 229, 012007

* Freeman H., 1961, IRE Trans. Electron. Comput., EC-10, 260

* Gheller C., Sartoretto F., Guidolin M., 2004, Mem. Soc. Astron. Ital. Suppl., 4, 55

* Goertz C. K., Ip W.-H., 1982, Planet. Space Sci., 30, 855

* Grandclément P., Bonazzola S., Gourgoulhon E., Marck J.-A., 2001, J. Comput. Phys., 170, 231

* Guillet T., Teissier R., 2011, J. Comput. Phys., 230, 4756

* Hachisu I., 1986, ApJS, 61, 479

* Hachisu I., Eriguchi Y., Nomaoto K., 1986, ApJ, 308, 161

* Hamerský J., Karas V., 2013, A&A, 555, A32

* Hashimoto M., Eriguchi Y., Arai K., Mueller E., 1993, A&A, 268, 131

* Horedt G. P., ed., 2004, Astrophysics and Space Science Library, Vol. 306, Polytropes – Applications in Astrophysics and Related Fields. Kluwer, Dordrecht

* Huré J.-M., 2005, A&A, 434, 1

* Ivory J., 1831, Phil. Trans. R. Soc. I, 121, 109

* James R., Kopal Z., 1963, Icarus, 1, 442

* Kellogg O. D., 1929, Foundations of Potential Theory. Springer, Berlin

* Kimura H., 1981, PASJ, 33, 273

* Kishimoto M., Hönig S. F., Antonucci R., Millour F., Tristram K. R. W., Weigelt G., 2011, A&A, 536, A78

* Kley W., 1996, MNRAS, 282, 234

* Kowalewsky S., 1885, Astron. Nachr., 111, 37

* Kuwahara F., 1988, Prog. Theor. Phys., 80, 449

* Lai D., Rasio F. A., Shapiro S. L., 1993, ApJS, 88, 205

* Lanza A., 1987, in Bruzzo U., Cianci R., Massa E., eds, General Relativity and Gravitational Physics. p. 67

* Lanza A., 1992, Class. Quantum Grav., 9, 677

* Mathews J. H., Fink K. D., 2004, Numerical Methods Using MATLAB. Prentice-Hall, Englewood Cliffs, NJ

* Murphy J. O., Fiedler R., 1985, Proc. Astron. Soc. Aust., 6, 219

* Nishida S., Eriguchi Y., 1994, ApJ, 427, 429

* Nishida S., Eriguchi Y., Lanza A., 1992, ApJ, 401, 618

* Odrzywołek A., 2003, MNRAS, 345, 497

* Ostriker J. P., Mark J. W.-K., 1968, ApJ, 151, 1075

* Petroff D., Horatschek S., 2008, MNRAS, 389, 156

* Reese D., 2006, Thesis, Université Paul Sabatier – Toulouse III

* Rieutord M., Espinosa Lara F., Putigny B., 2016, J. Comput. Phys., 318, 277

* Rouan D., Rigaut F., Alloin D., Doyon R., Lai O., Crampton D., Gendron E., Arsenault R., 1998, A&A, 339, 687

* Rucinski S. M., 1988, AJ, 95, 1895

* Schmitz F., 1990, A&A, 230, 509

* Singh M., 1970, Indian J. Pure Appl. Math., 1, 180

* Slaný P., Kovář J., Stuchlík Z., Karas V., 2013, ApJS, 205, 3

* Stuchlík Z., Slaný P., Kovář J., 2009, Class. Quantum Grav., 26, 215013

* Taniguchi K., Nakamura T., 2000, Phys. Rev. D, 62, 044040

* Tohline J. E., Hachisu I., 1990, ApJ, 361, 394

* Umemura M., Ikeuchi S., 1986, A&A, 165, 1

* Washimi H., Shibata S., Mori M., 1996, PASJ, 48, 23

* Wong C.-Y., 1974, ApJ, 190, 675

* Yoshida S., Eriguchi Y., 2006, ApJS, 164, 156

This paper has been typeset from a *TEX/LATEX* file prepared by the author.