© ESO 2020

1. Introduction

The smoothed particle hydrodynamics (SPH) technique has been widely used in astrophysics to study highly dynamical, geometrically distorted, and often catastrophic events, such as star formation (Springel & Hernquist 2003) and stellar encounters either from direct collisions (Lombardi et al. 1995) or with finite impact parameters (Davies et al. 1992). SPH has been used to model the gravitational radiation emission from the merging of white dwarfs (Lorén-Aguilar et al. 2005) and neutron stars (Rasio & Shapiro 1992). It has also been applied to simulate type Ia (García-Senz & Bravo 2005; Pakmor et al. 2012) and type II supernova explosions (Fryer & Warren 2002; Cabezón et al. 2018). It is also an important tool in the modeling of large-scale structure in the Universe (Springel 2005; Guedes et al. 2011).

Interestingly, the initial models in several of the scenarios above are simple barotropic configurations. In the particular case of the merging of two white dwarf (WD) stars, the final outcome could be a type Ia supernova (SNe Ia) explosion, being this double-degenerate scenario (DD) one of the favored production channels for these explosions (Hillebrandt et al. 2013). White dwarfs belonging to a compact binary system may have their rotational velocity substantially increased once the merger sets in, owing to the transfer of angular momentum from the accretion disk. Although it is thought that the accretion of matter from the companion star would lead to differential rotation (Yoon & Langer 2005), the assumption of rigid rotation is easy to implement and provides insight into many physical problems. Moreover, the rigid rotation hypothesis is adequate in those cases where the transport of angular momentum from the surface to the center of the star is very efficient (Piro 2008) as could be the case for degenerate objects like WDs. In particular, the fingerprint of the rotation in the thermonuclear explosion of a WD has been studied by Pfannes et al. (2010a,b) and García-Senz et al. (2018). In the latter, the SPH code SPHYNX (Cabezón et al. 2017) was used to simulate the explosion of a WD with mass ≃1 M_⊙ in fast rigid rotation. In that calculation, the initial model was built using the relaxation scheme proposed in the present work. The models above belong to another SNe Ia explosion channel, which involves a single WD accreting mass from a companion, non-degenerate star through the Roche-Lobe overflow. This second possibility is called the single-degenerate scenario (SD) channel for SNe Ia. Both scenarios (DD and SD) involve large amounts of angular momentum, so the question arises of how to adequately model self-gravitating fast-spinning rotators with SPH.

Another topic where having good initial rotating models is crucial is the study of the interaction between the accretion disk around millisecond pulsars with their magnetospheres. Numerical simulations of these scenarios require the correct modeling of steady accretion disks, characterized by the total mass, angular momentum, and particular rotation law. Although there have been a number of works dealing with the interaction between accretion disk and magnetosphere (Parfrey & Tchekhovskoy 2017, and references therein) none of them was carried out with particle-based codes. Clearly, it would be very useful to have steady accretion disks with the adequate resolution to perform simulations of these scenarios with SPH.

Unfortunately, there is no general procedure to build initial conditions for such self-gravitational, rotating equilibrium structures with SPH. A mapping procedure from an axisymmetric grid of points to a 3D distribution of equal-mass SPH particles was briefly described in Durisen et al. (1986) in connection with the study of rapidly rotating n = 3/2 polytropes. Smith et al. (1996) used a similar mapping procedure to simulate the development of a dynamical bar instability in a spinning polytrope with n = 3/2. However, the direct mapping from grid-points to particles usually leads to incomplete equilibrium because the models presented in those works display excessive numerical noise¹. A numerical scheme to handle disks in the presence of pressure gradients was discussed by Raskin & Owen (2016) and there are several recipes to approach the initial conditions prior to the dynamic merging of two WDs with SPH codes (Dan et al. 2011). There is, however, a rich literature concerning grid-based calculations of the equilibrium properties of barotropic, self-gravitating gases in rotation. An iterative algebraic method to obtain axisymmetric equilibrium structures is the self-consistent field (SCF) method developed by Ostriker & Mark (1968). Later, Hachisu (1986) successfully applied the SCF method to build zero-temperature spinning white dwarfs. Another approach was proposed by Eriguchi & Sugimoto (1981) in which the Poisson equation and the hydrostatic equilibrium equation (which includes the centrifugal force) are both used in integral form and iterated until convergence in a 2D grid of points. Nevertheless, these methods are not directly applicable to building initial models to carry out SPH calculations because the balance of gravity, pressure, and centrifugal forces is lost during the mapping procedure from the 2D ordered grid of points to a 3D distribution of particles.

As far as we know there is no public, well-documented, general procedure for obtaining self-gravitating structures in steady rotation with particle-based hydrodynamic codes. In this manuscript, we develop and test an easy yet practical novel relaxation scheme to build barotropic, P(ρ), 3D rotating structures in equilibrium under the SPH paradigm. In our proposal it is not necessary for the initial distribution of particles to match any prescribed SCF solution prior to relaxation. As an example, these relaxed structures can be used as suitable initial conditions to study the explosion of a rotating WD in the SD scenario, and as the outcome of the WDs merging in the DD scenario. In this work the relaxation scheme is applied to zero-temperature white dwarfs, to high-density polytropes, and to pseudo-Keplerian disks. To validate the method we provide quantitative comparisons with reference solutions, like those obtained by Hachisu (1986) and Eriguchi & Mueller (1985) using iterative time-independent methods such as the SCF scheme.

In Sect. 2 we describe the physical foundations of our proposal. Section 3 presents the application of the method to build 3D zero-temperature white dwarfs in rigid rotation. The scheme developed in Sect. 3 is also applied to the initial setting of two interacting WDs in a compact binary system in Sect. 4. The extension of the scheme to handle differential rotation in white dwarfs and polytropes is explained in Sect. 5. We apply the method to build pseudo-Keplerian accretion disks in equilibrium in Sect. 6. Finally, we present a summary of our findings and the prospects for future work in the conclusions in Sect. 7.

2. Relaxing rotating white dwarfs with SPH

In this section we aim to describe the physical basis of our relaxation scheme to obtain stable rotating structures of barotropic bodies, with particular emphasis in degenerate structures such as cold white dwarfs. Non-rotating equilibrium configurations of WDs can be obtained by relaxing a sample of particles with initial spherical coordinates (r, ϕ, θ). These SPH particles are radially distributed according to the density profile, but randomly in angles ϕ, θ. Usually, a damping force is added to the momentum equation so that after a few sound-crossing times the sample of particles relaxes to a stable configuration. It is worth noting, however, that as the mass of the WD approaches the Chandrasekhar-mass limit this equilibrium is not perfect and the degenerate star undergoes small radial oscillations. Fortunately, when simulating SNe Ia explosions the thermonuclear destruction of a massive WD is so fast that it is enough to keep the equilibrium only during a few sound-crossing times, t_sc (typically t_sc ∼ 0.46 s at ρ ≃ 10⁹ g cm⁻³).

In the absence of rotation, the structure of the WD after the relaxation process is spherically symmetric and follows the well-known solution of the Lane-Emden (LE) equation. Thus, for an equation of state (EOS) dominated by a zero-temperature electron gas (Chandrasekhar 1939),

$$\begin{array}{cc}& P_{\mathrm{e}}=a[x(2x^2-3){(x^2+1)}^{\frac{1}{2}}+3{sinh}^{-1}x],\hfill \end{array}$$(1)

$$\begin{array}{cc}& \rho =b{x}^3,\hfill \end{array}$$(2)

where a = 6.00 × 10²² dynes cm⁻², b = 9.82 × 10⁵ μ_e g cm⁻³, and x is the Fermi momentum of electrons in relativistic units. The only parameters determining the density and pressure profiles are the mass of the WD and the electron molecular weight {M_WD, μ_e}. Unfortunately, there is not a simple description that is equivalent to the Lane-Emden equation but for rotating stars. Assuming that the rotating white dwarf has axisymmetric geometry, its structure is set by the triad {M_WD,  J_WD,  μ_e} where J_WD is its total angular momentum. For differential rotators it is necessary to additionally specify the rotational law followed by the angular velocity Ω(s), where s is the distance to the rotation axis.

Once these parameters are defined, the maximum density (ρ_max) and radius (R_WD) of the resulting configuration will be uniquely determined. It is important to note that the final values (ρ_max and R_WD) are different from those obtained when rotation is not present, even for the same combination of {M_WD,  μ_e}. Once the mass and its composition are fixed, a rotating WD will have a lower ρ_max and larger R_WD in the rotating plane than a non-rotating WD.

When studying the efficiency of nuclear burning in SNe Ia explosions, the dominant magnitude is the density. Therefore, we fix ρ_max between rotating and non-rotating models which, as a consequence, implies a change in the total mass of the WD.

In this work, we consider the following rotation law:

$$\begin{array}{c}\hfill \mathrm{\Omega }(s)=\frac{{\mathrm{\Omega }}_{\mathrm{c}}}{{(1+\frac{s^2}{R_{\mathrm{c}}^2})}^m}·\end{array}$$(3)

Here s is the distance to the rotation axis, R_c is a parameter that sets the size of the central core with nearly rigid rotation Ω_c, and m is a parameter linked to the type of rotation (see Sects. 3 and 5). The current choices of m are m = 0 (rigid rotation), m = 1/2, and m = 1 (shellular). Rigid rotation is also attained for R_c ≫ R_WD, independently of m.

Our relaxation method works as follows. First, we choose the values of ρ_max, M_WD, and J_WD from Hachisu (1986) data tables (their Tables 4 and 5), and we take μ_e = 2 in the electron zero-temperature EOS. We then build an initial model with spherical symmetry (i.e., without rotation) and maximum density, ρ_max. This model is obtained after integrating the Lane-Emden equation with inner boundary condition ρ(r = 0) = ρ_max. We note that the total mass of the WD obtained from the LE equation, $M_{\mathrm{WD}}^{\mathrm{LE}}$, is not the same as the M_WD of the rotating model in Hachisu’s tables for the same central density ρ_max. The density profile is then mapped to a 3D distribution of equal-mass SPH particles, which ensures that the number density of the particles reflects the current density structure at any point. The mass of the particles is then re-scaled by a factor $M_{\mathrm{WD}}/M_{\mathrm{WD}}^{\mathrm{LE}}$, so that the total mass becomes M_WD. This changes the central density, but only transiently because its value is rapidly restored to ρ_max during the relaxation process. Next, a velocity is given to each SPH particle so that the total angular momentum is J_WD,

$$\begin{array}{c}\hfill J_{\mathrm{WD}}={\displaystyle \sum _b}{m}_b{s}_b{v}_b=\left[{\displaystyle \sum _b}{m}_b{s}_b^2{(1+\frac{s_b^2}{R_{\mathrm{c}}^2})}^{-m}\right]{\mathrm{\Omega }}_{\mathrm{c}}(t),\end{array}$$(4)

from which the instantaneous value of Ω_c(t) is obtained:

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{c}}(t)=\frac{J_{\mathrm{WD}}}{\left[{\sum }_b{m}_b{s}_b^2{(1+\frac{s_b^2}{R_{\mathrm{c}}^2})}^{-m}\right]}·\end{array}$$(5)

It should be noted that when R_c → ∞ the angular velocity becomes Ω(t) = J_WD/I_WD(t), where I_WD(t) is the time-dependent momentum of inertia of particles around the rotation axis, which is computed at each integration step during the relaxation. For other choices of R_c and m, Eq. (5) gives the correct angular velocity as a function of the generalized momentum of inertia of the system, provided J_WD is known at t = 0. This original and simple recipe is physically sound and robust because it always preserves the total mass and total angular momentum of the system, even for highly deformed axisymmetric structures.

Once Ω_c(t) is known, Eq. (3) gives Ω(s, t), so that the centripetal acceleration of each particle Ω × (Ω × r) is obtained. The particle distribution is henceforth relaxed in a co-moving reference frame. We then let the system freely evolve with the SPH code until equilibrium.

To relax the system and dissipate the spurious numerical noise stored in the velocity field the velocities are periodically set to zero (in the co-moving frame of reference). For example, all models in Table 1 with N_p = 5 × 10⁵ particles were obtained by setting v = 0 every Δt ≃ t_sc/3, being t_sc the sound-crossing time (t_sc ≃ 0.46 s in models with central density ρ_c = 10⁹ g cm⁻³) until time t ≃ 5t_sc/3. Afterwards the velocities were set to zero every Δt ≃ 0.8t_sc s. Such simple recipe works well and can be adapted to handle different particle numbers and densities. After several sound crossing times, the final equilibrium configuration is attained and only a residual numerical noise remains. To decide when the particle distribution has converged to a stable, time-independent configuration, each calculated model has to fulfill two criteria: that its central density and equatorial radius become constant (except small oscillations with amplitude 1–2% around their average value), and that the oscillations around the average value of the central density, the polar and equatorial radius, and the total kinetic, internal and gravitational energies remain lower than ≃2% when the system is left free to evolve in the current inertial frame (i.e., without re-setting the velocities to zero).

3. Isolated WDs with rigid rotation

As a first application we focus on the evolution towards equilibrium of models calculated assuming rigid rotation for both single WDs and double-degenerate stars in compact binaries. Some insights into how to handle differential rotation while obeying the rotation law given by Eq. (3) with m ≠ 0 are provided in Sect. 5. The simulations presented here, were carried out with the hydrodynamic code SPHYNX, using the electron zero-temperature EOS given by Eq. (1) with μ_e = 2. The calculations were performed with the default values of several parameters, for example α = 4/3 in the artificial viscosity (AV), with the Balsara limiters (Balsara 1995) to the AV turned on. The gravitational force is approached by a multipole expansion up to quadrupole terms, with the tolerance parameter θ (Hernquist 1987), set to θ = 0.6. An exponent parameter, p (0 ≤ p ≤ 1), allows choosing among different volume elements (VEs; Cabezón et al. 2017). The present simulations were obtained with p = 0, which is equivalent to the common choice VE = m/ρ.

A cautionary remark concerning rigid fast rotation close to Keplerian values is necessary here. For these models we observed that the combination of solid rigid rotation and a high rotational speed led to a harmful feedback between the centrifugal, gravitational, and pressure forces at the equator surface. Let us assume that there is a particle with a slight excess of angular velocity which displaces it farther in radius, thus increasing the centrifugal force (but weakening the gravity). This particle will move farther out, which again increases the centrifugal force, and so on. Even if a similar feedback only affects a handful of surface particles located on the equator, it is enough to spoil the convergence. If we want to make quantitative comparisons between our results and those by the SCF method for perfect rigid rotators, a solution is to (artificially) reduce a bit the gradient of pressure at the surface so that the feedback above is broken, even in the presence of low amounts of numerical noise. A simple way to do this while keeping the relaxation stable is to reduce the pressure in the low-density regions of the WDs by multiplying the electron pressure, Eq. (1), by a density cutoff:

$$\begin{array}{c}\hfill P=\{\begin{array}{cc}{P}_{\mathrm{e}}\hfill & \rho >{\rho }_{\mathrm{crit}}\hfill \\ P_{\mathrm{e}}\times \frac{\rho }{{\rho }_{\mathrm{crit}}}\hfill & \rho \le {\rho }_{\mathrm{crit}}\hfill \end{array}.\end{array}$$(6)

Here P is the pressure used in the calculations and ρ_crit = αρ_max, being ρ_max the maximum density, which in the models shown in Table 1 is attained at the center of the configuration. We empirically found that α = 5 × 10⁻⁴ works well in almost all the rigidly rotating models shown in Table 1. Slightly lower values, α = 2.5 × 10⁻⁴ and α = 2 × 10⁻⁴, were used to relax the more massive models A₂₁ and A₂₂, to reduce the error in the equatorial radius. Other magnitudes, such as the central density and total kinetic, internal, and gravitational energies, are rather insensitive to the precise value of α. We note again that this density cutoff is merely a way to facilitate the convergence of the relaxation. It affects a negligible amount of mass, and makes it possible to perform a direct comparison with the results obtained with the SCF method for rigid rotators. In real calculations, however, it is preferable not to consider absolute rigid rotation and to allow for some amount of differential rotation at the external layers. In this case, introducing the above density cutoff becomes unnecessary.

Here we carry out the relaxation of several WDs with different {M_WD, J_WD}. An estimation of the accuracy of the resulting equilibrium configurations is done by comparing our results to those of Hachisu (1986), which were obtained using the SCF method.

We provide a representative example of the evolution towards equilibrium of one of our rotating models in the left panel of Fig. 1, which depicts the evolution of the central density ρ_c (ρ_c = ρ_max in these models), equatorial radius, R_eq, and angular velocity, Ω. As we can see, the central density and angular velocity evolve in a similar manner. They start from relatively high values, decrease fast over a few tenths of a second, and stabilize at t ≃ 1 s. The equatorial radius follows the opposite trend: it is low at the beginning, then it rises quickly to asymptotically stabilize at t ≃ 2.0 s. From Fig. 1 it is obvious that the characteristic relaxation times of ρ_c and R_eq are quite different. To check that the star was stable we stop resetting the velocities to zero at t = 2.5 s, so that the system evolves freely in the current inertial frame of reference. The only source of dissipation being the artificial viscosity (AV) term, as given by Monaghan (1997), including the Balsara limiters. As can be seen, the central density remained stable during at least Δt ≃ 3 s (≃7 times the sound-crossing time). We show a color map slice of the density of model A₂₁ at t = 2.5 s in Fig. 1 (right). The white dwarf is neatly oblated with a polar-to-equator radius $F=\frac{R_{\mathrm{p}}}{R_{\mathrm{eq}}}=0.716$ which is ≃7.3% larger than that given by the SCF method. This is, however, the largest discrepancy found across all the calculated models in Table 1. The differences with respect the results by Hachisu remain, for the most part, below 4%.

Fig. 1. Simulation results of models A21 and A22 in Table 1. Left: evolution of central density (in 109 g cm−3), radius (in 108 cm), and angular velocity Ω. Solid lines correspond to model A21 (5 × 105 particles) and dashed lines to model A22 (2 × 106 particles). Once the periodic resetting to zero of the velocities is removed at t = 2.5 s, the central density remains stable during several sound-crossing times (tsc ≃ 0.4 s). Right: density color map and isodensity contours of a 2D meridional slice of the rotating WD at t = 2.5 s. The slice has a thickness of four times the local smoothing length (4 h), which represents roughly 10% of the total number of particles. Making cuts with the local value of h ensures a similar number of particles at any region of the color map.

Table 1 summarizes the relevant information of the calculated models. Models A refer to the SPHYNX calculations and H refer to the SCF models by Hachisu. Models A₁ and H₁ are non-rotating, with a mass close to the Chandrasekhar-mass limit. As we can see, the fit is excellent. The larger discrepancy, ≃4%, is in the radius of the configuration. Actually, this is something we expected because a white dwarf nearing the Chandrasekhar-mass limit does not have a well-defined scale length. We thus expect that the larger differences with respect the SCF method affect the equatorial radius and the oblateness of the WDs at central densities ρ₉ ≥ 1. Model A₂₂ in Table 1 is the same as A₂₁, but calculated with four times more particles. It leads to a stable model with similar relative errors in the central density, equatorial-to-polar radius and total energies, with respect to those of the reference model H₂. Therefore, our relaxation method is able to match the results by Hachisu in a wide range of stellar masses, 0.67 M_⊙ ≤ M_WD ≤ 1.44 M_⊙. The lower mass is close to that of a standard WD and the higher mass is actually at the Chandrasekhar-mass limit of a non-rotating white dwarf. The case of a super-Chandrasekhar-mass white dwarf stabilized by rotation is discussed in Sect. 5.2.

4. WDs in a double-degenerate binary system

A straightforward and timely extension of the relaxation procedure described above can be used to generate suitable initial models to study the DD production channel for SNe Ia. In the DD channel two WDs, settled in a compact orbit, move closer because of gravitational-wave radiative losses. At some point the gravitational pulling from the more massive WD breaks the lighter compact star and an accretion disk around the surviving WD forms. The further accretion of the debris would eventually provoke the explosion of the initially more massive white dwarf (Hillebrandt et al. 2013).

A key technical point of the simulations of the DD scenario is how to set the initial conditions immediately prior to the merging. The currently accepted standard procedure involves a two-step process (Rosswog et al. 2004; Dan et al. 2011). First, both stars are relaxed in isolation. Then, both stars are placed in a wide enough binary orbit in order to prevent any immediate mass transfer episodes. Subsequently, the binary system is evolved in the co-rotating frame where an artificial acceleration term is introduced in order to continuously shrink the orbital distance between the two white dwarfs. The orbital distance will be decreased at a sufficiently slow rate so the stars remain relaxed at all times, and so they slowly and continuously deform without introducing any spurious oscillations. The relaxation process will be finished when the secondary star starts to overflow its Roche lobe. This method is, however, computationally expensive and somewhat artificial because both WDs need to be relaxed all the time along the path from the initially detached position until they reach the onset of the merger (usually achieved through the introduction of additional artificial dissipative terms).

A somewhat more elegant and fast procedure is to relax both stars taking into account their angular momentum once they are already settled in a close orbit just prior to the merging, thus avoiding the intermediate relaxation stages. To do this we first set the orbit parameters so that the gravitational pulling onto the surface of the less massive WD becomes a sizable fraction (β) of its own gravity. That constraint leads to the following expression for the distance:

$$\begin{array}{c}\hfill D_{1,2}=R_2\times (1+\sqrt{\frac{M_{\mathrm{WD}1}}{\beta M_{\mathrm{WD}2}}})·\end{array}$$(7)

Here M_WD1 and M_WD2 are the masses of the more massive and lighter components, respectively. Once D_1, 2 is known, and assuming a circular orbit, we calculate the velocity of the center of mass of each star with respect to an inertial reference frame located at rest at the center of mass of the binary system. The total angular momentum of the binary system J_sys is calculated afterwards so that we can benefit from the scheme developed above in Sects. 2 and 3.

As an example, we have calculated the stable initial configurations in two cases. We first consider a system of twin white dwarfs with canonical masses M_WD1 = M_WD2 = 0.606 M_⊙ and parameter β = 1/4. The second case is for M_WD1 = 0.796 M_⊙, M_WD2 = 0.606 M_⊙, and the same value of β. The center of the compact stars is supposed to move in circles around the center of mass of the binary system, with angular velocity Ω. The total angular momentum of the system is J_sys = J_orb + J_spin. In this work we focus on pairs of white dwarfs which are tidally locked so that Ω_spin = Ω_orb = Ω. Therefore, once Ω is deduced from

$$\begin{array}{c}\hfill \mathrm{\Omega }=\sqrt{G\frac{(M_{\mathrm{WD}1}+M_{\mathrm{WD}2})}{D_{1,2}^3}},\end{array}$$(8)

the total angular momentum J_sys = J_orb + J_spin is easily obtained². The post-relaxed binary configuration in rigid rotation will be determined by {M_WD1, M_WD2, J_sys}. The EOS was that of a zero-temperature electron gas given by Eq. (1). A summary of the parameters used in these simulations is given in Table 2.

The evolution of the central densities of both WDs in model DD₂ is shown in Fig. 2 (left). The initial values of central density for both stars are considerably higher at the beginning of the relaxation process, as expected from spherically symmetric initial models. As the simulation proceeds the system rapidly achieves a stable configuration with lower densities due to the influence of rotational effects and tidal forces and to the re-ordering of the particles. This stable configuration is reached very quickly (∼1 s for the more massive star and ∼5 s for the other). At t = 20 the periodic reset to zero of the velocities was turned off without any appreciable effect on the evolution of the two stars. We also provide the equilibrium configuration of the WDs in the same figure (right panel), obtained in the co-moving non-inertial frame of reference located at the center of mass. Each dot represents an SPH particle within a thin cut along the plane Z = 0. The color-coding shows the logarithm of density. Comparing the two snapshots, it is clear that the initially spherical random particle distribution reaches a more ordered distribution where the least massive star is asymmetrically elongated in the direction of the more massive, as expected. After the relaxation we checked that the final structure remains stable during at least one complete orbit in the inertial frame. Although the final binary configuration is stable, it is on the verge of Roche-lobe overflow. To trigger the catastrophic merging of the WDs it is enough to slightly shrink the distance D_1, 2.

Fig. 2. Initial setting of the DD scenario. Left: evolution of the central density of both WDs during the relaxation (t ≤ 20 s) and the free evolving (t ≥ 20 s) stages. The orbital period is P ≃ 70 s. Right: slice in the orbital plane depicting the density color map of each white dwarf at t = 0 s (initial spherically symmetric configuration) and t = 20 s (final relaxed model). The slice has a thickness of four times the local smoothing length value (4 h). The center of mass of both configurations is located at (0, 0) km, but the snapshots have been shifted 14 000 km to the left and to the right to avoid the superposition of the images.

5. Differential rotation

It is feasible to apply the proposed relaxation method to handle differential rotation. We analyze here the impact of taking m ≠ 0 in Eq. (3) with different values of the parameter R_c. We first discuss the case of rotating polytropes of index n = 3/2 with central densities and sizes close to those characterizing neutron stars (NSs). The case of a zero-temperature rotating white dwarf with differential rotation is discussed later. In all cases the results were checked with well-known existing solutions obtained via other independent methods.

5.1. Case 1: differentially rotating high-density polytropes

To check the ability, as well as the limits, of our relaxation scheme to handle non-rigid rotation we chose a polytropic relation

$$\begin{array}{c}\hfill P=K{\rho }^{\gamma }\end{array}$$(9)

with γ = 5/3. In a non-rotating model this value produces a very stable spherically symmetric object whose structure is obtained after solving the Lane-Endem equation with n = 1/(γ −1) = 3/2.

Our initial pre-relaxed model is a spherically symmetric object with mass M = 2 M_⊙ and central density ρ_c = 10¹⁴ g cm⁻³. This combination sets the polytropic constant in Eq. (9) to K = 1.72 × 10¹⁰. The integration of the LE equation leads to an object of radius R_NS ≃ 38 km, roughly compatible with a NS size. Some amount of total angular momentum J was then added to the polytrope and the structure relaxed with the method explained in Sect. 2. The angular momentum J was chosen so that the square of the dimensionless angular momentum j,

$$\begin{array}{c}\hfill j^2=\frac{J^2}{4\pi {\mathrm{GM}}^{\frac{10}{3}}{\rho }_{max}^{-\frac{1}{3}}},\end{array}$$(10)

matches the values in Tables 1 and 2 in Eriguchi & Mueller (1985, hereafter EM) so that meaningful comparisons with their results can be made. We note, however, that the value of ρ_max in Eq. (10) initially differs from that in EM because we take ρ_max = ρ_c of the spherically symmetric polytrope, instead of the maximum density of the rotating structure used by EM. This initial choice of ρ_max is motivated by the fact that we do not know the true value of the density of the rotating structure prior to relaxation. Nevertheless, this does not pose a problem because the value of j can be recalculated with the value of ρ_max obtained once the relaxation has ended. This new value of j (Col. 6 in Table 3) was the one used to check our results with those given in the EM tables.

We calculated three models using j-constant (m = 1) and two more with v-constant (m = 1/2) rotational laws. In each case we consider R_c = AR_NS with A = 2 and A = 0.2 in Eq. (3), respectively. The case A = 2 is actually close to rigid rotation, but the case A = 0.2 is representative of models with large differential rotation. Table 3 shows a summary of the results, where the meaning of the columns is as follows: J (Col. 4) is the total angular momentum; ρ_max is the maximum density of the rotating body (not necessarily achieved at the center); E_k, E_I, and E_G is the total kinetic, internal, and gravitational energy, respectively; and E₀ is a normalization energy defined by

$$\begin{array}{c}\hfill E_0=\frac{{(4\pi G)}^2{M}^5}{J^2}·\end{array}$$(11)

The symbol V.T refers to the virial theorem defined as V.T =|(2E_k + E_G + 3(γ −1)E_I)/E_G|, Ω_c is the angular velocity at the center of the polytrope and, finally, F = R_p/R_eq is the polar-to-equatorial radius.

As shown in Table 3, the global energies and the polar-to-equator ratio agree with the results by EM within a few percentage points. Nevertheless, models B_n do not fulfill the virial theorem as well as they do in the EM_n calculations. There is almost a factor of ten difference between the two estimations. This discrepancy arises from the very different approach to the structure equilibrium of the rotating polytropes. In the case of the SPH calculation the equilibrium is approached dynamically in full 3D, while in the EM calculation the equilibrium equations are solved by means of a time-independent iterative procedure in a fixed 2D grid of points.

Figure 3 depicts the evolution of ρ_max (left panel) and Ω_c (right panel) for several of the models in Table 3. The profiles are similar to those shown in Fig. 1 of an isolated rotating white dwarf, but with a temporal scale of milliseconds instead of seconds. The final equilibrium value of Ω_c oscillates slightly around a stable value that, depending on the model, is in the range 10³ s⁻¹ ≤ Ω_c ≤ 10⁴ s⁻¹. This corresponds to periods of 0.6 ms ≤P ≤ 6 ms, typical of millisecond pulsars. Figure 4 shows several density color maps along meridional slices of models B_n. The first column depicts the color maps of pseudo-rigid rotators (A = 2), being the case m = 1 (model B₁) less elongated than m = 1/2 (model B₄). According to EM, this last case is very close to the critical value, R_p/R_eq = 0.5988, at which the gravitational and centrifugal forces become equal at the surface of the configuration when the rotation is 100% rigid (Eriguchi & Mueller 1985).

Fig. 3. Approach to the equilibrium of models described in Table 3 as a function of time: density (left panel) and Ωc (right panel).

Fig. 4. Color map of density in a meridional cut (with thickness 4h and axis units in km) of several models, as described in Table 3.

Models B₃ and B₅ (A = 0.2) host a strong differential rotation and become highly deformed, with ρ_max achieved far from the geometrical center of the configuration (see Fig. 4, right column). In particular, model B₃ is close to the limiting configuration at which our relaxation method gives satisfactory results when the trial initial configuration (prior to relaxation) is a spherically symmetric LE model. This occurs because of the large contrast between the angular velocity at the center and at the surface equator, which in model B₃ is Ω_c/Ω_eq ≃ 40. Obtaining suitable models for much higher angular velocity ratios would require a more refined trial initial configuration, for instance those obtained with the SCF method, and increased resolution, which is left for a future work.

5.2. Case 2: isolated WD with differential rotation

We consider again a spinning zero-temperature white dwarf, but this time with differential rotation. In particular, we focus on initial models obeying the rotation law given by Eq. (3) and m = 1/2. This choice leads to rigid rotation at the center (s ≤ R_c) of the configuration, but becomes Keplerian at distances s ≫ R_c. As previously stated, the case m = 1/2 is usually referred to as v-constant in the literature. As can be seen in models A_6, 7 in Table 1, we find a good agreement between our models and those reported by Hachisu for the case m = 1/2 (models H_6, 7 in Table 1). For these models R_c = 0.1 R_eq, where the value of R_eq was taken from the work by Hachisu (models 3 and 4 with ρ_max = 10⁹ g cm⁻³ in their Table 4). Another example of differential rotation, with an exponent intermediate between m = 1/2 and m = 1, is discussed in the following section.

We show a summary of our results for the two cases with m = 1/2 in Table 1 and Fig. 5. The two models differ from the SCF calculations by less than 5%, being stable enough for further hydrodynamic calculations. Figure 5 (left), shows the evolution of the maximum density and central value of the angular velocity for model A₇. At t ≃ 1 s the magnitudes ρ_max and Ω_c both become stable. The profile of the angular velocity Ω has been superposed on the density color map (right). It is clear that the maximum density is not located at the center of mass of the configuration, which is a typical signature of models with high angular momentum. The value of the angular velocity Ω is maximum at the center of the configuration and is very high, Ω(s = 0) ≃ 18 s⁻¹.

Fig. 5. Time evolution of the maximum density and central angular velocity in model A7, with differential rotation (left). Superposed on the density color map, with thickness 4h, is the 1D profile of angular velocity (green) along a 1D cut on the equator plane (right).

After being relaxed, the rotating white dwarf is allowed to evolve freely during ≃10 complete orbits of the fastest particle. During that time the central density remains constant (left panel in Fig. 5), while the central angular velocity Ω_c slightly oscillates around 18.4 s⁻¹. What is shown in Fig. 5 at t >  1.2 s is the average value of Ω_c for all particles with s ≥ 3 × 10⁶ cm. An estimate of Ω_c(s_a) is obtained from the tangential velocity of particle a, v_t(s_a, t) through Eq. (3):

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{c}}(s_a,t)=\frac{v_t(s_a,t)}{s_a(t)}\times {(1+\frac{s_a^2(t)}{R_{\mathrm{c}}^2})}^m·\end{array}$$(12)

6. Pseudo-Keplerian disks orbiting a central mass-point

Steady rotating disks are a universal phenomena in astrophysics (e.g., Frank et al. 2002) and often appear in merging stellar binary systems and during stellar and planetary formation (Armitage 2011). Because of the thermal pressure contribution, these disks cannot be described in a purely Keplerian way. Pressure gradient effects are important and have to be taken into account to adequately model the structure of steady or pseudo-steady disks. Nevertheless, the simulation of pressure supported accretion disks with any hydrodynamic method has been proven difficult, especially when keeping track of the structure during many orbital periods. In the particular case of SPH codes, there is an extra difficulty in building a stable enough initial distribution of particles with pressure gradients in differential rotation (Owen 2004).

With the aim of checking the abilities of our relaxation scheme, we have implemented the generalized disk test problem discussed by Raskin & Owen (2016). An idealized disk was assumed to have cylindrical geometry and can be studied in 2D. As shown by these authors, the gravitational potential ϕ(r), pressure (P(r) = Kρ^γ, with γ = 3/2), density, and tangential velocity (v_θ) profiles follow precise analytic relationships, which we reproduce here for completeness:

$$\begin{array}{cc}& \varphi (r)=-\frac{\mathrm{GM}}{{(r^2+r_{\mathrm{s}}^2)}^{\frac{1}{2}}}\hfill \end{array}$$(13)

$$\begin{array}{cc}& \rho (r)={\left[\frac{\mathrm{GM}(\gamma -1)}{K\gamma {(r^2+r_s^2)}^{\frac{1}{2}}}\right]}^{\frac{1}{(\gamma -1}}·\hfill \end{array}$$(14)

The expressions above correspond to a pressure-supported disk with zero angular momentum. To introduce rotation, a reduction factor f_p in the pressure is assumed so that the EOS is P(r) = Kρ^γ f_p. To keep the equilibrium, the loss in the pressure-gradient force is compensated by adding a centripetal force created by a tangential velocity field,

$$\begin{array}{c}\hfill v_{\theta }^2(r)=(1-f_{\mathrm{p}})\frac{GM{r}^2}{{(r^2+r_s^2)}^{\frac{3}{2}}}·\end{array}$$(15)

For this test we set GM = 1, r_s = 0.5, and f_p = 0.5. The value of K in the EOS is obtained from Eq. (14), assuming ρ(r = 0) = ρ₀ = 1, giving K = 2/3. The resulting disk is pseudo-Keplerian with an important contribution of pressure. Finally, we want to study both the resulting profiles of our relaxed initial models and the ability of SPHYNX to keep the disk in steady rotation during many orbits.

Because it is a 2D calculation, it is feasible to build an initial model by spreading the particles in an ordered array according to the density profile (Eq. (14)). As shown by Owen et al. (1998) and Raskin & Owen (2016), it is enough to distribute the particles in rings so that the radial separation between particles is adjusted in each annulus such that the angular separation between particles at a given radial coordinate is constant. In this work we propose a different procedure, which is capable of generating good initial models, but with glass-like structure.

First, we fit the angular velocity Ω(r) = v_θ(r)/r obtained from Eq. (15) with good accuracy, by the rotation law given by Eq. (3) with Ω_c(t = 0) = 2, m = 0.730, and r_c = 0.478. The total angular momentum of the disk, J_D, is estimated from the fitted Ω(r) and ρ(r), with 0 ≤ r ≤ R_D, where R_D = 10 is the adopted disk radius. The main features of our initial model are summarized in Table 4. A sample of N = 9 × 10⁴ particles³ was radially distributed in a plane according to the density profile, while their angular position 0 ≤ θ ≤ 2π is set at random. Afterwards, the particle sample is relaxed with the procedure described in Sect. 2.

Figure 6 shows the density, pressure, and tangential velocity profiles of the particles once the relaxation has ended (t = 0 cases). A density color map showing the central region of the disk with r ≤ 3 is depicted in the left semi-plane of the bottom right panel in the same figure. In spite of the disordered glass-like pattern, the distribution of the SPH-particles in the disk matches the analytic profile very well. Furthermore, after the relaxation, this pseudo-Keplerian disk should retain its structure during as many orbits as possible. It was shown by Raskin & Owen (2016) that many current SPH schemes fail to preserve the disk features after about a dozen periods of the particle at maximum velocity. Satisfactory results were obtained with the conservative reproducing kernel (CRKSPH) method by Frontiere et al. (2017) in combination with the AV switches by Cullen & Dehnen (2010).

Fig. 6. Profiles of density, tangential velocity, and pressure of the disk after the relaxation period (green lines) compared to the analytic values (in red) and after t = 200 s calculated with the Balsara limiter (blue) or the AV switches by Read & Hayfield (2012) (magenta). Bottom right panel: density color map in the central region of the disk. The left semi-plane (x <  0) is for the relaxed model at t = 0 s; the right semi-plane (x ≥ 0) shows the density color map after t = 200 s, with the Balsara limiter. Unlike in Raskin & Owen (2016), the SPH particles are still settled in a glass-like configuration.

The stability of the disk obtained with our relaxation scheme was further checked with the SPHYNX code. We want to know if the integral approach to the gradients in combination with either the Balsara limiter (Balsara 1995) or the AV switches (Cullen & Dehnen 2010; Read & Hayfield 2012) is capable of keeping the disk in steady state over many orbital periods. For this calculation the EOS was changed to P = (γ −1)ρu (γ = 5/3), where u is the specific internal energy of the gas. The evolution of the internal energy was obtained by evolving the corresponding energy equation.

We have tracked the evolution until t = 200 s, so that the results can be compared to those obtained by Raskin & Owen (2016) with the CRKSPH method. At t = 200 s the particle with the highest velocity has completed ≃40 orbits around the center of the disk. As shown in Fig. 6, the deviation of ρ(s),v_t(s), and P(s) from either the analytic or the t = 0 s profile is very small. Particularly good is the fit of the tangential velocity, which still matches the initial profile after ≃40 orbits. Nevertheless, a closer look into the density profile in the central region indicates that its value is slowly growing with time. We agree with Raskin & Owen (2016), that such small growth of the central density is an artifact of the particular implementation of the artificial viscosity. It is slightly less pronounced when the switches to the AV are in command of the dissipation. It is worth noting, however, that the pressure profile remains very close to the analytic value.

The x ≥ 0 semi-plane in the bottom right panel in Fig. 6 depicts the density color map at t = 200 s in the central part of the disk. As we can see, the differences with the color map at t = 0 s (x <  0 semi-plane) are small, being only relevant at the very center of the configuration. Additionally, the model still retains a good glass-like granulation, with no indications of pairing instability after ≃40 revolutions.

7. Conclusions

A common problem in simulating the evolution of compact objects with SPH is that there is not a general procedure to obtain stable initial models when these objects are spinning quickly. In this paper we propose and test an easy and versatile relaxation scheme for building stable rotating models of self-gravitating bodies whose EOS is of barotropic type. The hypothesis of barotropic EOS allows to handle many interesting objects such as isothermal white dwarfs, neutron stars, and polytropic structures. As detailed in Sect. 2, our method relies on the exceptional angular momentum conservation properties of the SPH technique.

We apply our method to get stable rotating configurations of zero-temperature WDs with different masses and total angular momentum. We were able to build stable models of rotating white dwarfs and polytropes spanning a wide mass range, 0.67 M_⊙ ≤ M_WD ≤ 2 M_⊙, with both rigid and differential rotation. The main magnitudes (central density; total kinetic, internal, and gravitational energies; equatorial and polar radius) match the semi-analytical results by Hachisu (1986) and Eriguchi & Mueller (1985), in general within a few percentage points. Given the current uncertainties in the particular rotation law followed by these compact objects, this precision is enough to explore many issues concerning their evolution, either for explosion or collapse scenarios. Additionally, we show that our method is able to produce stable configurations when it is applied to a pair of white dwarfs orbiting in a compact binary system. This last scenario is of capital importance to understanding the double-degenerate route to SNe Ia explosions.

Finally, we have applied the method to produce a steady accretion disk with cylindrical symmetry. Even though the assumption of cylindrical geometry is not realistic, this configuration has the advantage that it has an analytical solution for comparison, while still retaining many of the features of real disks. The ensuing relaxed disk is characterized by a glass-like distribution of particles in differential rotation, partially supported by pressure effects and by rotation. As shown in Sect. 6, this pseudo-Keplerian disk is in a good steady state, being stable over several dozen orbits of the particle at maximum velocity.

Presently, two cautionary remarks have to be taken into account before applying our procedure to astrophysical calculations. First, the centrifugal force has to remain weaker than the gravitational pull at the equator when approaching the equilibrium configuration. Otherwise a harmful tradeoff between gravity (∝1/s²) and the centrifugal force (∝s) may appear which finally leads to the ejection of particles located in the outer shells along the equatorial plane. Second, in differential rotation the ratio of the angular velocity at the center to that at the surface, prior to relaxation, should be not too extreme. In particular, the combination of a shellular (m ≃ 1) rotation with a high total angular momentum and low values of R_c in Eqs. (3) and (5) lead to toroidal structures where r_p/r_eq → 0 (Eriguchi & Mueller 1985). Difficulties may also appear when trying to relax bodies with extremely large density contrasts between the core and the surface layers, owing to their very different characteristic timescales. In these cases, the choice of spherically symmetric Lane-Endem models as initial trial configurations is not adequate because they are too far from the final equilibrium structure and the relaxation to a steady state could be hard or even impossible. This difficulty can probably be overcome using better initial trial models, for example those obtained with the SCF method or even envisaging the final toroidal structure as the result of the coalescence of two spherically symmetric objects (as was done in Sect. 4).

The practical cases studied in this work cover only a small subset of the potential applications. An immediate additional application of our method may consist in building 3D stable post-Newtonian models of fast-rotating neutron stars, either isolated or in binary systems. This task should not be difficult because the EOS of cold neutron stars can also be described assuming zero temperature.

Prospects for future extensions of this work may include the impact of finite temperature gradients in the relaxed rotating structures. Currently, the steady spinning polytropes obtained in Sect. 5.1 with a barotropic EOS, P = Kρ^γ, can be easily mapped into an axisymmetric distribution of temperature if we assume a realistic EOS, P(ρ, T), i.e., that of a non-degenerate ideal gas. This suggests that relaxing a trial initial model with a realistic EOS besides a superimposed axisymmetric temperature profile, T(s), from the onset is feasible.

Acknowledgments

This work has been supported by the MINECO Spanish project AYA2017-86274-P and the Generalitat of Catalonia SGR-661/2017 (DG), and by the Swiss Platform for Advanced Scientific Computing (PASC) project SPH-EXA: Optimizing Smooth Particle Hydrodynamics for Exascale Computing (RC and DG). JMBI acknowledges the support by the FPU fellowship and wants to thank the Ministerio de Educación, Cultura y Deporte from Spain. The authors acknowledge the support of sciCORE (http://scicore.unibas.ch/) scientific computing core facility at University of Basel, where part of these calculations were performed.

References

All Tables

All Figures

Fig. 1. Simulation results of models A21 and A22 in Table 1. Left: evolution of central density (in 109 g cm−3), radius (in 108 cm), and angular velocity Ω. Solid lines correspond to model A21 (5 × 105 particles) and dashed lines to model A22 (2 × 106 particles). Once the periodic resetting to zero of the velocities is removed at t = 2.5 s, the central density remains stable during several sound-crossing times (tsc ≃ 0.4 s). Right: density color map and isodensity contours of a 2D meridional slice of the rotating WD at t = 2.5 s. The slice has a thickness of four times the local smoothing length (4 h), which represents roughly 10% of the total number of particles. Making cuts with the local value of h ensures a similar number of particles at any region of the color map.

In the text

Fig. 2. Initial setting of the DD scenario. Left: evolution of the central density of both WDs during the relaxation (t ≤ 20 s) and the free evolving (t ≥ 20 s) stages. The orbital period is P ≃ 70 s. Right: slice in the orbital plane depicting the density color map of each white dwarf at t = 0 s (initial spherically symmetric configuration) and t = 20 s (final relaxed model). The slice has a thickness of four times the local smoothing length value (4 h). The center of mass of both configurations is located at (0, 0) km, but the snapshots have been shifted 14 000 km to the left and to the right to avoid the superposition of the images.

In the text

	Fig. 3. Approach to the equilibrium of models described in Table 3 as a function of time: density (left panel) and Ωc (right panel).
In the text

	Fig. 4. Color map of density in a meridional cut (with thickness 4h and axis units in km) of several models, as described in Table 3.
In the text

	Fig. 5. Time evolution of the maximum density and central angular velocity in model A7, with differential rotation (left). Superposed on the density color map, with thickness 4h, is the 1D profile of angular velocity (green) along a 1D cut on the equator plane (right).
In the text

Fig. 6. Profiles of density, tangential velocity, and pressure of the disk after the relaxation period (green lines) compared to the analytic values (in red) and after t = 200 s calculated with the Balsara limiter (blue) or the AV switches by Read & Hayfield (2012) (magenta). Bottom right panel: density color map in the central region of the disk. The left semi-plane (x <  0) is for the relaxed model at t = 0 s; the right semi-plane (x ≥ 0) shows the density color map after t = 200 s, with the Balsara limiter. Unlike in Raskin & Owen (2016), the SPH particles are still settled in a glass-like configuration.

In the text