# On the maximum mass of magnetised white dwarfs

D. Chatterjee, A.F. Fantina, N. Chamel, J. Novak, M. Oertel

► To cite this version:

D. Chatterjee, A.F. Fantina, N. Chamel, J. Novak, M. Oertel. On the maximum mass of magnetised white dwarfs. Monthly Notices of the Royal Astronomical Society, 2017, 469, pp.95-109. 10.1093/mn-ras/stx781. in2p3-01390118

HAL Id: in2p3-01390118

https://in2p3.hal.science/in2p3-01390118v1

Submitted on 27 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>

# On the maximum mass of magnetized white dwarfs

D. Chatterjee,<sup>1,2*</sup> A. F. Fantina,<sup>3,4*</sup> N. Chamel,<sup>3*</sup> J. Novak<sup>1*</sup> and M. Oertel<sup>1*</sup>

Accepted 2017 March 28. Received 2017 March 14; in original form 2016 February 24

## ABSTRACT

We develop a detailed and self-consistent numerical model, assuming a poloidal magnetic field, for extremely magnetized white dwarfs, which have been proposed as progenitors of overluminous Type Ia supernovae. This model can describe fully consistent equilibria of magnetic stars in axial symmetry, with rotation, general relativistic effects and realistic equations of state (including electrostatic interactions and taking into account Landau quantization of electrons due to the magnetic field). We study the influence of each of these ingredients on to the white dwarf structure and, in particular, on their maximum mass. We perform an extensive stability analysis of such objects, with their highest surface magnetic fields reaching ~10<sup>13</sup> G (at which point the star adopts a torus-like shape). We confirm previous speculations that, although very massive strongly magnetized white dwarfs, could potentially exist, the onset of electron captures and pycnonuclear reactions may severely limit their stability. Finally, the emission of gravitational waves by these objects is addressed – our results suggest that detection by the currently planned space-based detector *Laser Interferometer Space Antenna* is unlikely.

**Key words:** equation of state – magnetic fields – methods: numerical – white dwarfs.

## 1 INTRODUCTION

White dwarfs (WDs) are the stellar remnants of low- and intermediate-mass stars, that is, stars with masses < 10 M<sub>☉</sub> (M<sub>☉</sub> being the mass of our Sun; Shapiro & Teukolsky 1983). The interest in WD properties and in particular their mass-radius relation has been renewed by the recent discovery of overluminous Type Ia supernova (SNIa; Howell et al. 2006; Scalzo et al. 2010; Maoz, Mannucci & Nelemans 2014). The progenitors of such events are thought to be ‘super-Chandrasekhar’ WDs with a mass > 2 M<sub>☉</sub> (see e.g. Hillebrandt et al. 2013). These SNIa may result from the merger of two massive WDs or from the explosion of rapidly differentially rotating WDs (see e.g. Howell et al. 2006). Alternatively, it has been argued that WDs endowed with a strong magnetic field could be potential candidates for the progenitors of these peculiar SNIa (Kundu & Mukhopadhyay 2012; Das & Mukhopadhyay 2012a,b, 2013a,b; Das, Mukhopadhyay & Rao 2013). As a matter of fact, very massive strongly magnetized WDs were proposed much earlier by Shul’man (1976); this work, however, does

not seem to have attracted much attention. The existence of magnetized supermassive WD is of interest not only for astrophysics, but also for cosmology since SNIa have been used as ‘standard candles’ to measure cosmological distances assuming a unique astrophysical scenario for these events.

The mass-radius relation for non-magnetic WDs was established way back in 1935 by the pioneering work of Chandrasekhar, assuming a simple degenerate electron Fermi gas equation of state (EoS) at zero temperature. Subsequently, other works attempted to construct more realistic models of WDs. Apart from the magnetic field, different effects could alter the Chandrasekhar limit for the maximum mass of a WD.

(i) Inclusion of finite temperature (Marshak 1940) and slow rotation (Ostriker & Hartwick 1968) did not result in significant deviations from the Chandrasekhar WD mass-radius relation.

(ii) The electrostatic interaction between electrons and ions introduced by Hamada & Salpeter (1961) was found to lower the electron pressure resulting in a softer EoS and hence leads to slightly less massive configurations.

(iii) The effect of general relativity on the structure of non-magnetic WDs was found to be non-negligible, reducing by a few per cent the WD maximum mass (see e.g. Ibáñez 1983, 1984; Shapiro & Teukolsky 1983; Rotondo et al. 2011; De-Hua, He-Lei &

<footnote n="1"><sup>1</sup>LUTH, Observatoire de Paris, PSL Research University, CNRS, Université Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, F-92195 Meudon, France</footnote>

<footnote n="2"><sup>2</sup>Laboratoire de Physique Corpusculaire, ENSICAEN, 6 Boulevard Maréchal Juin, F-14050 Caen Cédex, France</footnote>

<footnote n="3"><sup>3</sup>Institut d'Astronomie et d'Astrophysique, CP-226, Université Libre de Bruxelles (ULB), B-1050 Brussels, Belgium</footnote>

<footnote n="4"><sup>4</sup>Grand Accélérateur National d'Ions Lourds (GANIL), CEA/DRF - CNRS/IN2P3, Bvd Henri Becquerel, F-14076 Caen, France</footnote>

<footnote>* E-mail: dchatterjee@lpccaen.in2p3.fr (DC); anthea.fantina@ganil.fr (AFF); nchamel@ulb.ac.be (NC); jerome.novak@obspm.fr (JN); micaela.oertel@obspm.fr (MO)</footnote>
Xiang-Dong 2014; Boshkayev et al. 2015; Chamel & Fantina 2015; Bera & Bhattacharya 2016).

The existence of magnetized WDs, with magnetic fields $B \gtrsim 10^6$ G (Shapiro & Teukolsky 1983), first predicted in the 1940s by Blackett (1947), was confirmed in 1970 by Kemp (see e.g. Jordan 2009, for a review). At the time of this writing, hundreds of magnetized WDs have already been discovered (Kepler et al. 2013). While surface magnetic fields up to about $10^9$ G can be inferred from Zeeman spectroscopy and polarimetry, or cyclotron spectroscopy (see e.g. Wickramasinghe & Ferrario 2000, for a review), internal magnetic fields are not directly accessible by observations. Therefore, very strong magnetic fields in the core of WDs, as high as $B \sim 10^{13}$ G (Fujisawa, Yoshida & Eriguchi 2012), cannot be ruled out from a simple estimate of the energetics (i.e. scalar virial theorem) or the conservation of magnetic flux of the progenitor star. The question therefore arises of how such a strong magnetic field affects the structure of WDs.

The study of the mass-radius relation of a magnetized WD has a long history and it was recognized early on that the impact of the magnetic field on both the stellar radius and mass could be large. However, simplifying assumptions have been made for convenience. For instance, the pioneering work by Ostriker & Hartwick (1968) considered a vanishing magnetic field at the surface of the star and neglected any magnetic field effect on the EoS as well as general relativistic effects and electrostatic interactions. The works of Adam (1986), Das & Mukhopadhyay (2012a) and Kundu & Mukhopadhyay (2012) focused on the effect of the magnetic field on the EoS, taking into account the Landau quantization of the electron gas, but were based on a Newtonian description of the star's structure in spherical symmetry, that is, neglecting the deformation of the star by the magnetic field. A similar approach was followed by Suh & Mathews (2000), where, however, the (spherical) stellar structure was calculated from the general relativistic Tolman-Oppenheimer-Volkoff equations. The deformations of the star under the influence of a strong magnetic field have been recently studied by Bera & Bhattacharya (2014) in the Newtonian framework. Das & Mukhopadhyay (2015) and Bera & Bhattacharya (2016) have computed the mass-radius relation of magnetic WDs for different magnetic field geometries within a general relativistic framework using the publicly available XNS code (Bucciantini, Pili & Del Zanna 2014).<fnref n="1" /> On the other hand, Landau quantization effects were neglected.

According to all these previous studies, the maximum possible mass of strongly magnetized WDs lies substantially above the Chandrasekhar limit, and more importantly above the mass inferred from observations of overluminous SNIa; see also Cheoun et al. (2013), Herrera & Barreto (2013), Herrera et al. (2014), Federbush, Luo & Smoller (2015) and Belyaev et al. (2015). Leaving aside the origin of such strong magnetic fields – more than three orders of magnitude higher than currently observed fields – the question as to whether such stars can possibly exist still needs to be further examined. Coelho et al. (2014), Chamel, Fantina & Davis (2013), Chamel et al. (2014) and Chamel & Fantina (2015) discussed various microscopic and macroscopic instabilities. In particular, Chamel et al. (2013, 2014) pointed out that the stability of very massive magnetized WDs will be limited not only by general relativity, but also by electron captures by nuclei and pycnonuclear reactions in the stellar core. These microscopic processes, whose onset depends on

the magnetic field strength (Chamel & Fantina 2015), should thus be carefully considered in modelling massive magnetized WDs.

In this study, we address these stability issues. To this end, we shall construct consistent global general relativistic models of strongly magnetized WDs, taking into account the magnetic field effects on both the EoS and the WD equilibrium structure equations. This work follows directly from our previous study of the structure of strongly magnetized neutron stars (Chatterjee et al. 2015). We shall account for the onset of electron capture and pycnonuclear reactions, including the shifts on the threshold densities and pressures due to electrostatic interactions and to the strong magnetic field. Recently, during the refereeing process of our paper, a very similar independent study by Otoniel et al. (2016) appeared in the literature; this work will be briefly discussed.

The paper is organized as follows. In Section 2, we describe the realistic dense matter EoS we considered, taking into account both the effects of electrostatic interactions and the presence of the magnetic field. In Section 3, we summarize the modified Maxwell equations and the hydrodynamic equations in presence of a magnetic field in general relativity and Newtonian theory. In Section 4, we describe the results concerning the influence of an extremely strong magnetic field on the WD structure. Section 5 describes the effects of rotation when combined with the magnetic field. Section 6 gives details about the influence of different microphysical ingredients on the WD model. Section 7 discusses physical instabilities that may occur within our equilibrium models. Section 8 makes a short assessment on the detectability of gravitational waves (GWs) emitted by strongly magnetized WDs. Finally, in Section 9, we give our conclusions. Further details on our model can be found in the Appendix.

## 2 EQUATION OF STATE OF STRONGLY MAGNETIZED WHITE DWARFS

The model adopted here was initially developed for the description of the outer crust of strongly magnetized neutron stars (Lai & Shapiro 1991; Chamel et al. 2012). The interior of magnetic WDs is assumed to be composed of fully ionized atoms. Moreover, we assume that the internal temperature has dropped below the crystallization temperature $T_m$ so that ions are arranged on a regular crystal lattice. The crystallization temperature is defined as (see e.g. Haensel, Potekhin & Yakovlev 2007)

$$T_m = \frac{e^2}{a_e k_B \Gamma_m} Z^{5/3}, \quad (1)$$

where $e$ is the proton electric charge, $a_e = (3/(4\pi n_e))^{1/3}$ is the electron-sphere radius, $n_e$ is the electron number density, $k_B$ is the Boltzmann's constant and $\Gamma_m$ is the Coulomb coupling parameter at melting. For simplicity, we shall consider crystalline structures made of only one type of ions $A_Z X$, with mass number A and atomic number Z (typically carbon or oxygen). In unmagnetized matter, $\Gamma_m \approx 175$ (Haensel et al. 2007). In the presence of a strong magnetic field, the crystalline structure is more stable since $\Gamma_m \lesssim 175$ (Potekhin et al. 2013). Because the electron Fermi temperature is typically much smaller than $T_m$, electrons are highly degenerate. In the following, we shall therefore neglect thermal effects. The matter energy density is given by

$$\mathcal{E} = n \frac{M'(A, Z)c^2}{A} + \mathcal{E}_e + \mathcal{E}_L - n_e m_e c^2, \quad (2)$$

where $M'(A, Z)$ is the mass of the nucleus $A_Z X$ (including the rest mass of nucleons and Z electrons), c is the speed of light, n is the

<footnote n="1">¹ http://www.arcetri.astro.it/science/ahead/XNS/index.html.</footnote>
baryon number density, $\mathcal{E}_e$ the energy density of electrons, $m_e$ the electron mass and $\mathcal{E}_L$ the lattice energy density. The mass density $\rho$ is defined as $\rho = mn$, where $m$ denotes the average mass per nucleon (approximated here by the unified atomic mass unit). The last term in equation (2), the electron rest-mass energy, is included to avoid double counting. For magnetic fields below ~$10^{17}$ G, the nuclear masses remain essentially unchanged (Peña et al. 2011; Stein et al. 2016). As in Lai & Shapiro (1991) and Chamel et al. (2012), we therefore assume that nuclear masses are the same as in the absence of magnetic fields. The nuclear mass $M'(A, Z)$ can be obtained from the tabulated atomic mass $M(A, Z)$ from the 2012 atomic mass evaluation (Audi et al. 2012) after subtracting out the binding energy of the atomic electrons (see equation A4 in Lunney, Pearson & Thibault 2003):

$$M'(A, Z)c^2 = M(A, Z)c^2 + 1.443\,81 \times 10^{-5} Z^{2.39} \\ +1.554\,68 \times 10^{-12} Z^{5.35}, \qquad (3)$$

where both masses are expressed in units of MeV c⁻².

In the presence of a strong magnetic field, the electron motion perpendicular to the field is quantized into Landau levels (see e.g. Chapter 4 in Haensel et al. 2007). Ignoring the small electron anomalous magnetic moment (see e.g. section 4.1.1 in Haensel et al. 2007, and references therein), and treating electrons as a relativistic Fermi gas, the energies of Landau levels are given by

$$\epsilon_v = \sqrt{c^2 p_z^2 + m_e^2 c^4 (1 + 2v_L b_*)} \qquad (4)$$

$$v_L = n_L + \frac{1}{2} + \sigma, \qquad (5)$$

where $n_L$ is any non-negative integer, $\sigma = \pm 1/2$ is the spin, $p_z$ is the component of the momentum along the magnetic field and $b_*$ = $b/b_{crit}$ represents the magnetic field strength $b$ in units of the critical magnetic field $b_{crit}$ defined by

$$b_{crit} = \frac{m_e^2 c^3}{e\hbar} \approx 4.4 \times 10^{13} \text{ G}. \qquad (6)$$

For a given magnetic field strength $b_*$, the number of occupied Landau levels is determined by the electron number density $n_e$

$$n_e = \frac{2b_*}{(2\pi)^2 \lambda_e^3} \sum_{v_L=0}^{v_L^{\max}} g_v x_e(v_L), \qquad (7)$$

$$x_e(v_L) = \sqrt{\gamma_e^2 - 1 - 2v_L b_*}, \qquad (8)$$

where $\lambda_e = \hbar/m_e c$ is the electron Compton wavelength and $\gamma_e$ is the electron Fermi energy in units of the electron rest-mass energy,

$$\gamma_e = \frac{\mu_e}{m_e c^2}, \qquad (9)$$

with $\mu_e = d\mathcal{E}_e/dn_e$, while the degeneracy $g_v$ of a Landau level is $g_v = 1$ for $v_L = 0$ and $g_v = 2$ for $v_L \ge 1$.

The electron energy density $\mathcal{E}_e$ and corresponding electron pressure $P_e$ are given by

$$\mathcal{E}_e = \frac{b_* m_e c^2}{(2\pi)^2 \lambda_e^3} \sum_{v_L=0}^{v_L^{\max}} g_v (1 + 2v_L b_*) \psi_+ \left[ \frac{x_e(v_L)}{\sqrt{1 + 2v_L b_*}} \right], \qquad (10)$$

and

$$P_e = \frac{b_* m_e c^2}{(2\pi)^2 \lambda_e^3} \sum_{v_L=0}^{v_L^{\max}} g_v (1 + 2v_L b_*) \psi_- \left[ \frac{x_e(v_L)}{\sqrt{1 + 2v_L b_*}} \right], \qquad (11)$$

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Matter contribution to the EoS (pressure $P$ versus mass density $\rho$) for a $^{12}$C WD, for different magnetic field strengths $b_*$, where $b_* = b/b_{crit}$ and $b_{crit}$ is defined in equation (6).</figcaption></figure>

respectively, where

$$\psi_{\pm}(x) = x\sqrt{1+x^2} \pm \ln(x + \sqrt{1+x^2}). \qquad (12)$$

A magnetic field will be referred to as strongly quantizing if only the lowest level $v_L = 0$ is filled, or equivalently whenever $n_e < n_{eb}$, where

$$n_{eb} = \frac{b_*^{3/2}}{\sqrt{2\pi^2 \lambda_e^3}}, \qquad (13)$$

which corresponds to the mass density

$$\rho_b = \frac{m}{y_e} \frac{b_*^{3/2}}{\sqrt{2\pi^2 \lambda_e^3}}, \qquad (14)$$

where $y_e = n_e/n$. Conversely, for a given mass density $\rho$ the strongly quantizing regime corresponds to magnetic field strengths

$$b_* > \left( \frac{\rho y_e \lambda_e^3 \sqrt{2\pi^2}}{m} \right)^{2/3} \approx 180 (2y_e \rho_{10})^{2/3}, \qquad (15)$$

where $\rho_{10} = \rho/(10^{10} \text{ g cm}^{-3})$.

According to the Bohr-van Leeuwen theorem (Van Vleck 1932), the lattice energy density is independent of the magnetic field, neglecting the small contribution due to the quantum zero-point motion of ions (Baiko 2009). Considering point-like ions arranged in a body-centred cubic (bcc) structure (Kozhberov 2016), the lattice energy density is given by

$$\mathcal{E}_L = C e^2 n_e^{4/3} Z^{2/3}, \qquad (16)$$

where $C \approx -1.44$ (Baiko, Potekhin & Yakovlev 2001). The corresponding lattice contribution to the pressure is given by

$$P_L = \frac{\mathcal{E}_L}{3}. \qquad (17)$$

In the recent work of Otoniel et al. (2016), nuclei are supposed to be arranged in a simple cubic lattice. However, such a structure is known to be unstable (Born 1940). The matter contribution to the pressure is therefore (the contribution from nuclei is negligible at low temperatures)

$$P = P_e + P_L. \qquad (18)$$

The matter pressure $P$ is plotted in Fig. 1 as a function of the mass density $\rho$ for different core magnetic field strengths $b_*$ for a WD composed of $^{12}$C. The kinks correspond to the complete filling of Landau levels. From this figure, it can already be anticipated that
the effect of the magnetic field on the EoS will have only a very small influence on the WD structure, except for magnetic fields $b_* \gg 1$.

The magnetization $M$ is defined as

$$M = \left. \frac{\partial P}{\partial b} \right|_{\mu}. \qquad (19)$$

The baryon chemical potential $\mu$, which coincides with Gibbs free energy per nucleon,

$$g = \frac{\mathcal{E} + P}{n}, \qquad (20)$$

is given by

$$\mu = \frac{M'(A, Z)c^2}{A} + \frac{Z}{A}(\mu_e - m_e c^2 + \mu_L), \qquad (21)$$

where

$$\mu_L = \frac{4 \mathcal{E}_L}{3 n_e}. \qquad (22)$$

Complete general expressions of the magnetization as well as of its derivatives can be found in Appendix A.

In the following, we shall consider stellar cores made of either $^{12}$C or $^{16}$O. Indeed, these elements are the most likely fusion products of the helium burning phase.

## 3 STELLAR STRUCTURE EQUATIONS

In this section, we outline the formalism to construct equilibrium WD configurations starting from the EoS described in Section 2. As general relativity has been shown to have a non-negligible effect on the maximum WD mass, we compute the structure of WDs within both Newtonian theory of gravity and the general relativistic framework, for comparison (see also Section 7.1). The matter properties enter the energy-momentum tensor that acts as a source of the Einstein equations. The equilibrium is determined by the conservation of energy and momentum, derived from the condition of vanishing divergence of the energy-momentum tensor. We then compute equilibrium numerical models of WDs to obtain the global structure properties such as mass or radius.

The Einstein-Maxwell and equilibrium equations described hereafter are solved using spectral methods within the numerical library LORENE<fnref n="2" /> (Gourgoulhon et al. 2016). We apply the numerical scheme originally designed for strongly magnetized neutron stars (see Chatterjee et al. 2015, for further details) to study strongly magnetized WDs. As elaborated in Bocquet et al. (1995), the electromagnetic field tensor can be chosen to be either purely poloidal or purely toroidal. Within the numerical scheme of LORENE, the magnetic field configuration is purely poloidal by construction. This is convenient for comparison with the observed surface poloidal fields as estimated from the spin-down measurements of pulsars. However, for neutron stars, it has been suggested that amplification of the seed poloidal magnetic field could result in the generation of strong toroidal fields. Following the same arguments, one may argue that the same holds for strongly magnetized WDs. Further, both purely poloidal and toroidal magnetic field configurations are unstable, resulting in the rearrangement of the configuration to a mixed one. That means, a purely poloidal magnetic field configuration is not necessarily the most general one. Studies of magnetized WDs in a general magnetic field configuration have been recently performed

employing the xns code (Bera & Bhattacharya 2014, 2016; Das & Mukhopadhyay 2015). These studies showed that large WD masses can be supported by a purely toroidal field, but these configurations are also unstable. For the twisted torus mixed configuration, the star deformations are almost completely due to the poloidal part of the magnetic field and the toroidal component is always energetically subdominant (at most 10 per cent of the total energy; Bucciantini, Pili & Del Zanna 2016). Therefore, our simplifying assumption of purely poloidal magnetic field should be relevant in the WD case too.

### 3.1 Energy-momentum tensor in presence of a magnetic field

In Chatterjee et al. (2015), some of us explicitly derived the generalized expression for the thermodynamic average of the microscopic energy-momentum tensor, which is required for the study of the macroscopic structure of a compact object, as recapitulated below:

$$\begin{aligned} \langle T^{\mu\nu} \rangle = (\mathcal{E} + P) & u^{\mu} u^{\nu} + P g^{\mu\nu} \\ & + \frac{1}{2} (F_{\tau}^{\nu} M^{\tau\mu} + F_{\tau}^{\mu} M^{\tau\nu}) \\ & - \frac{1}{\mu_0} \left( F^{\mu\alpha} F_{\alpha}^{\nu} + \frac{g^{\mu\nu}}{4} F_{\alpha\beta} F^{\alpha\beta} \right). \end{aligned} \qquad (23)$$

where $\mathcal{E}$ is the matter energy density, $P$ is the matter pressure (see Section 2), $u^\mu$ is the fluid four-velocity and $\mu_0$ is the vacuum magnetic permeability in SI units. The electromagnetic field strength tensor is derived from the electromagnetic potential 1-form $A_\mu$ through

$$F_{\mu\nu} = \frac{\partial A_{\nu}}{\partial x^{\mu}} - \frac{\partial A_{\mu}}{\partial x^{\nu}}, \qquad (24)$$

and $M^{\mu\nu}$ is the magnetization tensor (see Chatterjee et al. 2015, for complete derivation).

The first two terms on the right-hand side of equation (23) can be identified as the pure (perfect fluid) matter contribution, followed by the magnetization term and finally the usual electromagnetic field contributions to the energy-momentum tensor. For isotropic media, the magnetization is aligned with the magnetic field. Thus one may write $M_{\alpha\beta} = \frac{\chi}{\mu_0} F_{\alpha\beta}$, in terms of the dimensionless scalar

$$\chi = \frac{\mu_0 M}{b}. \qquad (25)$$

Finally, the electric and magnetic fields as measured by the Eulerian observer have only two non-vanishing components each (see Bonazzola et al. 1993):

$$E_r = \frac{1}{N} \left( \frac{\partial A_r}{\partial r} + N^\psi \frac{\partial A_\psi}{\partial r} \right), \qquad (26a)$$

$$E_{\theta} = \frac{1}{N} \left( \frac{\partial A_{\theta}}{\partial \theta} + N^{\psi} \frac{\partial A_{\psi}}{\partial \theta} \right), \qquad (26b)$$

$$B_r = \frac{1}{Cr^2 \sin \theta} \frac{\partial A_{\psi}}{\partial \theta}, \qquad (26c)$$

$$B_{\theta} = - \frac{1}{C \sin \theta} \frac{\partial A_{\psi}}{\partial r}, \qquad (26d)$$

where $C$, $N$ and $N^\psi$ are gravitational potentials defined by the metric in equation (27). Note that we denote the magnetic field norm in the matter comoving frame by $b$, to distinguish it from that in the Eulerian frame, denoted by $B$ with components $B_r$ and $B_\theta$.

<footnote n="2">² http://www.lorene.obspm.fr.</footnote>
## 3.2 Einstein–Maxwell and equilibrium equations

To construct numerical models of magnetic WDs, we adapt the general relativistic scheme following previous works by Bocquet et al. (1995) and Bonazzola et al. (1993). Under the assumptions of stationarity, axisymmetry and circularity of the space–time, we employ the maximally sliced quasi-isotropic coordinates in which the line element can be written as

$$ds^2 = -N^2 dt^2 + C^2 r^2 \sin^2 \theta (d\phi - N^\psi dt)^2 + D^2 (dr^2 + r^2 d\theta^2), \quad (27)$$

where $N, N^\psi, C$ and $D$ are functions of coordinates $(r, \theta)$. With these properties and coordinate choice, the Einstein equations reduce to a set of four elliptic partial differential equations for the four gravitational potentials, in which part of the source terms is derived from the energy-momentum tensor. These equations are then solved for a given matter content. Note that, in this section, Latin letters i, j, ... are used for spatial indices only, whereas Greek ones $\alpha, \mu, ...$ denote the space–time indices.

Whereas the homogeneous Maxwell equation is automatically verified by the expression (24) for the electromagnetic field tensor, the inhomogeneous one must be modified to include the contribution from the magnetization (Chatterjee et al. 2015):

$$\frac{1}{\mu_0} \nabla_{\mu} F^{\nu\mu} = j_{\text{free}}^{\nu} + \nabla_{\mu} \mathcal{M}^{\nu\mu}, \quad (28)$$

where $j_{\text{free}}^{\nu}$ is the free current that generates the electromagnetic field, as opposed to the bound current contribution associated with the magnetization (last term in the above equation). For a discussion about these free currents, which are not driven by the bulk motion of matter, see Bonazzola et al. (1993).

To obtain equilibrium configurations, one must solve the coupled Einstein–Maxwell equations, together with the condition of conservation of energy and momentum $\nabla_{\mu} T^{\mu\nu} = 0$, giving rise to the equilibrium condition in our stationary case. It was already illustrated in Chatterjee et al. (2015) that on inclusion of the magnetization current density in the equation of equilibrium, the Lorentz force associated with the magnetization current cancels with the magnetization contribution arising from the pressure gradient (key equations are recapitulated in Appendix B). Therefore, there is no explicit magnetization appearing in the equilibrium condition, and the first integral of fluid stationary motion remains unchanged retaining the same form as in the case without magnetization (see details in Appendix B):

$$\ln h(r, \theta) + \nu(r, \theta) - \ln \Gamma(r, \theta) + \Phi(r, \theta) = \text{const}, \quad (29)$$

where $h$ is the enthalpy (as we assume zero temperature, $h = g$ the Gibbs free energy, introduced in equation 20), $\nu = \log N$ a gravitational potential function, $\Gamma$ is the Lorentz factor and $\Phi$ is the electromagnetic term associated with the Lorentz force (see Bonazzola et al. 1993, for a discussion):

$$\frac{\partial \Phi}{\partial x^i} = - \frac{F_{i\mu} j_{\text{free}}^{\mu}}{\mathcal{E} + P}. \quad (30)$$

## 3.3 Newtonian limit

As WDs are not strongly relativistic objects, their structure is often calculated in a non-relativistic Newtonian framework. We summarize here the basic stellar structure equations that the general relativistic equations reduce to in the Newtonian limit.

Whereas for Maxwell equations one must consider the special-relativistic form, the Newtonian limit of the Einstein equations is the Poisson equation for the gravitational potential $\phi$:

$$\nabla^2 \phi = 4\pi G\rho, \quad (31)$$

where $G$ is the gravitational constant and $\rho$ is the mass density. Thus in the non-relativistic limit, the factor $\nu$ tends to the Newtonian gravitational potential $\phi$ (Gourgoulhon 2011).

The first integral of motion (equation 29), described in the previous section, reduces to a simple form in the Newtonian limit in the absence of an electromagnetic field

$$\hat{h} + \phi - \frac{1}{2} U^2 = \text{const}, \quad (32)$$

where $\hat{h} = (\hat{\mathcal{E}} + P) / \rho$ is the non-relativistic enthalpy (i.e. $\mathcal{E} = \rho c^2 + \hat{\mathcal{E}}, \hat{\mathcal{E}}$ being the internal energy, not taking into account rest-mass energy density) and $U$ the fluid velocity in the $\varphi$-direction. Note here that, as shown, e.g. by Gourgoulhon (2011) this is not the classical Bernoulli equation, but a first integral of motion, which is valid only if the fluid is in pure circular motion, while the Bernoulli theorem is valid for any stationary flow. Moreover, the expression for the Bernoulli theorem provides a constant along each fluid line, which may vary from one fluid line to another, whereas the first integral, equation (32), gives a constant valid in the entire star.

In the presence of free currents $j_{\text{free}}^{\mu}$, see equation (28), and of the electromagnetic field they generate, a force term of the form $f_i = F_{i\sigma} j_{\text{free}}^{\sigma}$ is added to the equation of stationary motion (Bonazzola et al. 1993). Now, using the definition of $F_{\mu\nu}$, equation (24), together with the Newtonian limit (i.e. the metric potentials $N, C, D \to 1$ and $N^\psi \to 0$) in equations (26), and writing the 4-currents $j^\mu$ as a charge density $j_{\text{free}}^0$ and a current vector $j_{\text{free}}$ separately, one gets the standard form for the force:

$$f = j_{\text{free}}^{0} \mathbf{E} + \mathbf{J}_{\text{free}} \times \mathbf{B}. \quad (33)$$

As we assume the matter inside the star to be a perfect conductor, we have the relation:

$$\mathbf{E} = -\mathbf{v} \times \mathbf{B}, \quad (34)$$

with $\mathbf{v} = \mathbf{\Omega} \mathbf{e}^\psi$ the fluid velocity$^3$ (thus, $\mathbf{U} = \mathbf{\Omega} \text{rsin} \theta$). The Newtonian interpretation of our requirement of having a first integral of motion, with the introduction of the potential $\Phi$, equation (30), is that the electromagnetic force be potential-like:

$$f = (\mathbf{J} - \mathbf{\Omega} \mathbf{e}^\psi j^0) \times \mathbf{B} = \rho \nabla \Phi. \quad (35)$$

The equilibrium equation in the Newtonian case then reads

$$\hat{h} + \phi - \frac{1}{2} U^2 + \Phi = \text{const}. \quad (36)$$

As demonstrated before in the relativistic case, note again that the above equations for equilibrium do not contain any contribution from the magnetization. They thus differ from those given in Bera & Bhattacharya (2014), where the magnetization has been artificially included.

## 3.4 Numerical methods

From a numerical point of view, the model detailed above reduces to (i) an equilibrium equation in the form of a first integral of motion

<footnote>$^3$ $\mathbf{\Omega}$ is the angular velocity and $\mathbf{e}^\psi$ is the local triad vector associated with the $\varphi$-coordinate.</footnote>
as equation (36) in the Newtonian case, or equation (29) in general relativity; (ii) elliptic-type partial differential equation(s) for the gravitational field, as the Poisson equation (31) in the Newtonian case, or a set of four non-linear Poisson-like equations in the general relativistic case and (iii) the inhomogeneous Maxwell equation (28), which are of Poisson type too. Equations for the Newtonian model have been shown in Section 3.3; in the general relativistic case, these equations have been derived and discussed in detail in Bonazzola et al. (1993) and Chatterjee et al. (2015). In both cases, one has to solve a set of elliptic partial differential equations, which are eventually non-linear and with non-compact support. To do so, we use spectral methods (see e.g. Grandclément & Novak 2009) based on the decomposition of fields on to a basis of known functions; here we use multidomain Chebyshev expansion for the radial coordinate $r$ and Legendre polynomials of the type $P_{\ell}(\cos \theta)$ ($m = 0$ spherical harmonics) for the angular coordinate.

The numerical technique to solve a scalar or vector Poisson equation has been described in Grandclément et al. (2001), and is a standard part of the LORENE library (Gourgoulhon et al. 2016). The overall strategy to solve for a given WD numerical model is the same as that used for neutron stars and outlined in Bonazzola et al. (1993); it is based on the fixed point algorithm to deal with non-linear terms in the equations. This strategy has been improved by Bocquet et al. (1995), where several accuracy tests are presented. The numerical neutron star models were later incorporated into the LORENE library by Novak & Marcq (2003) and updated to be used together with a magnetic-field-dependent EoS by Chatterjee et al. (2015). It is the version of LORENE library used and described by Chatterjee et al. (2015) that we use in this study too.

For a given numerical model, one must provide (in addition to an EoS, eventually depending on the magnetic field strength): a central value for the enthalpy $\hat{h}_c = \hat{h}(r=0)$, which shall set the WD mass; a rotation frequency $\Omega$, which determines the fluid velocity; and a current function relating the free currents $j_{\text{free}}^{\mu}$ and the electromagnetic potential $A_{\mu}$, whose amplitude determines the magnetic field strength and dipole moment. For more details, see Bonazzola et al. (1993) and Chatterjee et al. (2015). Note that, in addition to the standard convergence tests, the accuracy of these models can be assessed through the so-called *relativistic virial theorem* (Bonazzola & Gourgoulhon 1994; Gourgoulhon & Bonazzola 1994), which is an independent probe of the global accuracy of the numerical solution. Results given in this work typically achieved a global accuracy, thus measured by the virial indicators, lower than 1 per cent and usually around $10^{-5}$.

## 4 STELLAR STRUCTURE AND MAGNETIC FIELD

With the numerical setup described in Section 3, we construct models of magnetic $^{12}$C and $^{16}$O WDs in equilibrium and determine their structure for various magnetic field strengths. After comparing our results with those obtained in previous studies, we show that the relevant quantity to consider for assessing the global stability of strongly magnetized WDs is the magnetic dipole moment rather than the magnetic field strength. We compute the maximum mass and discuss the limitations of our model. Since the magnetic field dependence of the EoS is very weak (see the discussion in Section 6.1), we will neglect magnetization and magnetic field dependence of the EoS for simplicity within the calculations (unless otherwise stated).

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Masses *M* in units of solar mass (M⊙) versus equatorial radii *R*<sub>eq</sub> for ¹²C WDs in Newtonian gravity along five sequences of fixed central magnetic field *B*<sub>*</sub>, ranging from 0 to 10. *B*<sub>*</sub> = *B*/*b*<sub>crit</sub> and *b*<sub>crit</sub> is defined in equation (6).</figcaption></figure>

### 4.1 Role of a strong poloidal magnetic field on the white dwarf structure

Based on a non-rotating model in Newtonian gravity, we first vary central enthalpy along sequences of fixed central magnetic fields $B_*$, where $B_* = B/b_{crit}$, $b_{crit}$ being defined in equation (6) and *B* is the magnetic field measured by the Eulerian observer. The value remaining fixed along a given sequence corresponds to the value of magnetic field at the centre of the star. Indeed, for each choice of free current distribution, equation (28), one obtains a given magnetic field distribution inside the star and thus one can determine the value of central magnetic field strength.

Fig. 2 displays the mass *M* (in units of solar mass M⊙) as a function of equatorial radius $R_{\text{eq}}$ (in km) for such sequences in Newtonian gravity, for a $^{12}$C WD, as an illustrative example. For low magnetic fields, the sequences follow the mass–radius relation as for non-magnetic WDs, ultimately reaching the Chandrasekhar limit. As the field strength is increased, the deviation from the non-magnetic curve increases, resulting in configurations with higher mass at lower densities. This is clearly evident from Fig. 3, where for the same sequences we show the masses for varying central mass density $\rho_c$ (in g cm⁻³). These results in Newtonian gravity are well known and can be compared with those of Suh & Mathews (2000) and Bera & Bhattacharya (2014) and are found to be in accordance. For field strengths of about $10 b_{crit}$, masses as large as ~1.8 M⊙ are found to be supported by the magnetic field for $^{12}$C WDs (similar results are found for $^{16}$O WDs).

These figures clearly demonstrate that WD masses larger than the standard Chandrasekhar limit for the non-magnetic case can be supported by strong magnetic fields. However, the strongly magnetized sequences ($B_* = 1, 10$) seem to be gravitationally unstable, if one uses the standard criterion for non-magnetic stars $\partial R_{\text{eq}} / \partial M > 0$ or $\partial M / \partial \rho_c < 0$ (e.g. Shapiro & Teukolsky 1983). The important point here is that, in order to look for such instabilities along families of two-parameter equilibria, one must be careful in choosing the quantity that should be kept constant when varying the mass (Sorkin 1982). Indeed, just as in the case of a rotating body, it is the angular momentum and not the angular velocity that is the constant
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Masses *M* in units of solar mass (M⊙) as a function of central density ρ<sub>c</sub> for <sup>12</sup>C WDs in Newtonian gravity along five sequences of fixed central magnetic fields B<sub>⋆</sub>, ranging from 0 to 10. B<sub>⋆</sub> = *B*/b<sub>crit</sub> and b<sub>crit</sub> is defined in equation (6).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** <sup>12</sup>C WD masses *M* versus radii *R*<sub>eq</sub> along sequences of fixed magnetic moments *D* in Newtonian gravity.</figcaption></figure>

of motion, it can be shown that the conservation of magnetic flux implies that magnetic moment is a constant of motion and not the magnetic field. Therefore, sequences of fixed dipole magnetic moments *D* should be considered, and not of fixed central magnetic field value as in Figs 2 and 3. This dipole moment is computed from the asymptotic behaviour of the magnetic field, obtained from our numerical model (see Bocquet et al. 1995, for a definition).

In Fig. 4, we therefore construct mass–radius relations for increas-
ing central enthalpy along sequences of fixed dipole magnetic moments *D*. As seen from this figure, even very highly magnetized sequences, which can reach two solar masses configurations, are gravitationally stable in Newtonian gravity. Effects of relativistic gravity (general relativity) and the justification for choosing Newtonian gravity in the rest of this work shall be discussed later in Section 7.1.

Finally, as it is not the magnetic moment that is an astrophys-
ically observable quantity, but rather the polar magnetic field that

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Polar magnetic field *B*<sub>P</sub> versus magnetic dipole moment *D* for two numerically computed sequences of <sup>12</sup>C Newtonian WDs of baryon masses 0.7 and 1.3 M<sub>☉</sub>.</figcaption></figure>

is derived from the measurement of the Zeeman effect on stellar
spectra, we here give numerically computed relations between both
quantities, so that one can get an idea of typical magnetic fields aris-
ing for the magnetic dipolar moments quoted in this study. For two
exemplary typical values of WD masses (0.7 and 1.3 solar masses),
we display in Fig. 5 the values of polar magnetic fields B<sub>P</sub> (in G) for
a <sup>12</sup>C WD corresponding to the range of magnetic moments consid-
ered in this work. These relations are obtained from sequences of
our numerical models, while there is in general no analytic relation
between both quantities. Indeed, as in these models the polar mag-
netic field is given by a superposition of all multipolar momenta,
the relation with the dipole moment (which is the leading moment,
asymptotically) can only be computed numerically. From the ob-
servational point of view, one must here keep in mind that surface
magnetic fields of isolated magnetic WDs are observed to lie within
the range 10<sup>3</sup>–10<sup>9</sup> G (Ferrario, de Martino & Gänsicke 2015).

**4.2 Maximally distorted stellar configurations**

With stronger magnetic fields, the anisotropic Lorentz force acting
on the free currents (see equation 33) causes the stellar structure
to increasingly deviate from spherical symmetry, finally reaching a
toroidal shape. At this point, the code reaches the limit of its nu-
merical capability for describing the stellar surface, and hence more
massive stable WD configurations cannot be achieved within this
numerical framework. As an example, in Fig. 6, we illustrate the
maximally distorted shape of the stellar surface of a <sup>12</sup>C WD for a
magnetic dipole moment of 3 × 10<sup>34</sup> A m<sup>2</sup>. The induced magnetic
field at the pole rises up to 3.1 × 10<sup>13</sup> G (equivalently B<sub>*</sub> ~ 1),
whereas the magnetic field at the centre reaches B<sub>*</sub> ~ 10. The ratio
of magnetic to fluid pressure computed at the centre of the WD is
found to be ~0.79 for <sup>12</sup>C WDs and 0.83 for <sup>16</sup>O WDs. The end-
points of the curves for two-solar masses magnetized WDs in Fig. 4
correspond to these maximally distorted configurations, for which
the numerical approach starts being inaccurate. Strictly speaking,
this does not represent an instability, and more massive stars with
torus-like shape are not a priori excluded (see also the discussion in
the case of neutron stars by Cardall, Prakash & Lattimer 2001). Nev-
ertheless, it should be kept in mind that for these highly magnetized
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Enthalpy isocontours of a static stellar configuration in the (x, z) plane of a $^{12}$C WD for a magnetic dipole moment of $3 \times 10^{34}$ A m$^2$. The mass of the star is 1.99 $M_{\odot}$ and the polar magnetic field $B_{\perp} \sim 3 \times 10^{13}$ G.</figcaption></figure>

stars, the polar field is four orders of magnitude above the highest currently observed values and that, in particular, it is not clear in which way such a huge magnetic field could be formed within a WD. The stability of such toroidal stars has not been checked, either.

Another point that we considered here is a purely poloidal magnetic field configuration. Purely toroidal configurations might lead to much higher masses (see e.g. Das & Mukhopadhyay 2015). However, it remains to be shown that these configurations are stable. In addition, under the fossil field hypothesis, numerical simulations performed by Braithwaite & Spruit (2004) have shown that the WD magnetic field is probably dominated by the poloidal field. Finally, the results for poloidal-dominated mixed configurations are fairly similar to the results for purely poloidal fields discussed here (see e.g. Bera & Bhattacharya 2016). In Section 7, we will discuss other instabilities that may limit the mass of a magnetized WD.

## 5 STELLAR STRUCTURE AND ROTATION

It is already a well-known result that rotation provides centrifugal force that can support larger masses in compact stars. We estimated the relative increase in mass of a rotating magnetic WD compared to a non-rotating one. The observed rotation periods of isolated magnetic WDs span the range from 725 s to decades or centuries (Ferrario et al. 2015). Magnetic cataclysmic variables exhibit shorter periods, with AE Aqr having a period of 33 s. For such values of the rotation period, we find that the relative difference between the masses of non-rotating and rotating WDs at a given density and magnetic dipole moment is negligible (of the order of $10^{-4}$ or lower). On the contrary, the magnetic field plays a major role in supporting massive WDs. To demonstrate this fact, we display in Fig. 7 the maximum mass of a WD rotating with a period of 725 s (depicted by a solid line; these results are indistinguishable from those obtained for a non-rotating WD and that rotating at 33 s),

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Masses for Newtonian magnetic WDs composed of $^{16}$O rotating with a period of 725 s (solid line) and at Kepler frequency (dashed line) as a function of fixed magnetic moments $D$.</figcaption></figure>

as well as a WD rotating at the Kepler frequency corresponding to the mass shedding limit (dashed line) for increasing values of the magnetic dipole moment. We have chosen to plot the results for $^{16}$O WDs, but similar results are obtained for $^{12}$C WDs. It is evident that even for the fastest rotating observed WDs, rigid rotation does not lead to any significant increase in their mass. More importantly, the effects of rotation are comparatively less important for strongly magnetized WDs, even at Kepler frequency, as the main force resisting gravity is no longer the centrifugal buoyancy, but the Lorentz force due to the magnetic field.

On the other hand, differentially rotating compact stars can be significantly more massive than their non-rotating or uniformly rotating counterparts (see e.g. Subramanian & Mukhopadhyay 2015). However, the support due to differential rotation would ultimately be cancelled by magnetic braking and/or viscosity. These processes drive the star into uniform rotation (Shapiro 2000). As a matter of fact, Bonazzola et al. (1993) showed that a magnetic star in a stationary state is necessarily rigidly rotating.

## 6 STELLAR STRUCTURE AND THE EQUATION OF STATE

In this section, we perform a systematic analysis of the role of the EoS on the WD structure, and more specifically the effects of the magnetic field as well as electron-electron and electron-ion interactions.

### 6.1 Magnetic field dependence of the equation of state

As discussed in Section 2, the magnetic field influences the properties of dense matter in magnetic WDs through (i) Landau quantization (i.e. magnetic-field-dependent EoS) and (ii) the magnetization $M$. In order to assess the relative importance of these effects, we compute the structure of magnetic WDs for a fixed central magnetic field $B_c = b_{crit}$ considering the three different cases.

(i) Case 1: EoS ignoring magnetic field effects.

(ii) Case 2: EoS including Landau quantization but neglecting magnetization.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Mass *M* versus radius *R*<sub>eq</sub> for ¹²C Newtonian WDs for fixed central magnetic field *B*<sub>c</sub> = *b*<sub>crit</sub> defined in equation (6), with and without magnetic field dependence in EoS and magnetization. See the text for details about cases 1, 2 and 3.</figcaption></figure>

(iii) Case 3: EoS including both Landau quantization and magnetization.

As shown in Fig. 8 for the illustrative case of a ¹²C WD, the corresponding mass–radius relationships are practically indistinguishable, in accordance with the results previously obtained by Bera & Bhattacharya (2014). Considering this result, along with Fig. 1, the discussion in Section 4.2, and the fact that the maximum values of the magnetic field inside the star remain low, only of the order of $B_* \lesssim 10$, even for the maximally distorted configurations, we will adopt the simplest case 1 for the following calculations.

## 6.2 Electron–electron and electron–ion interactions

Most previous investigations of magnetic WDs have employed the EoS of an electron Fermi gas that was originally considered in the work by Chandrasekhar (1931) [see e.g. Das & Mukhopadhyay (2012b), Kundu & Mukhopadhyay (2012), Bera & Bhattacharya (2014), Das & Mukhopadhyay (2015) and Franzon & Schramm (2015)]. Bera & Bhattacharya (2016) and Otoniel et al. (2016) have recently computed the structure of magnetic WDs including the lattice correction to the EoS. Bera & Bhattacharya (2016) employed the expression obtained by Salpeter (1961) in the Wigner–Seitz approximation, while Otoniel et al. (2016) used the expression for a simple-cubic lattice. However, as already pointed out in Section 2, this lattice type is unstable. In all our calculations presented so far (bcc), lattice corrections were taken into account. But, in order to assess the importance of these effects, we now set $C = 0$ in equation (16) (see Fig. 9), where the results without lattice effects are compared with the results taking them into account, for the illustrative case of ¹⁶O WDs. We find that electron–electron and electron–ion interactions lead to a slightly lower maximum mass for WDs, as previously discussed by Chamel & Fantina (2015) for non-magnetic WDs. The maximum mass of ¹²C and ¹⁶O non-magnetic WDs thus decreases from ~ 1.44 $M_⊙$ to ~ 1.41 $M_⊙$ and ~ 1.40 $M_⊙$, respectively (see also Tables 1 and 2). On the other hand, lattice effects can have a much larger impact on the

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Mass *M* versus radius *R*<sub>eq</sub> for Newtonian ¹⁶O WDs for different values of fixed magnetic moments with (solid lines) and without (dashed lines) lattice effects.</figcaption></figure>

stellar radius, e.g. for a non-magnetic star with a mass of 1.4 $M_⊙$, it changes from 1799 km without lattice effects to 639 km on inclusion of lattice interactions. These effects are less pronounced in strongly magnetized WDs.

## 7 STABILITY OF STRONGLY MAGNETIZED WHITE DWARFS

In this section, we discuss several physical mechanisms that could affect the WD stability, such as the use of general relativity for the description of gravity, electron-capture instability, and pycnonuclear instability. A summary of the different maximum mass values, considering the mechanisms discussed in this section, can be found in Table 1 for ¹⁶O WDs and in Table 2 for ¹²C WDs.

### 7.1 General relativistic instability

Although most mass–radius models for WDs have been constructed in Newtonian gravity, it is well known that computing equilibrium configurations within general relativity reduces the maximum mass for the case of non-magnetic WDs (e.g. Ibáñez 1983, 1984; Rotondo et al. 2011; Boshkayev et al. 2015; Chamel & Fantina 2015). We compute the structure of magnetized WDs both in Newtonian theory and in general relativity. To estimate how general relativistic effects may limit the maximum masses, we plot in Fig. 10 the mass–radius relationship for the non-magnetic as well as the magnetic ¹²C WDs, considering both Newtonian theory and general relativity. The mass used in general relativity is the so-called *gravitational mass*, which is the mass felt by a test-particle orbiting around the WD. For a complete definition, details of its computation and the definition of the circumferential equatorial radius used there too; see Bonazzola et al. (1993). It is clear from this figure that general relativity has a non-negligible effect in limiting the maximum mass of non-magnetic WD. We observed that on inclusion of general relativistic effects, the maximum mass of WD was reduced from 1.41 to 1.38 $M_⊙$.

For the case of strongly magnetized WDs, it is found that it is not general relativity that plays a crucial role in limiting the maximum masses but the magnetic field that distorts the WD beyond the equilibrium poloidal shape. The reason is that general relativity affects
<table><caption>**Table 1.** Maximum masses (in M⊙) of strongly magnetized $^{16}$O WDs for various magnetic moments (in A m$^2$), limited by different mechanisms: B denotes the onset of torus-like shape due to magnetic field, EC stands for electron-capture instability and pycno is short for pycnonuclear fusion reactions.</caption><thead><tr><th rowspan="2">Magnetic moment<br>(A m<sup>2</sup>)</th><th colspan="3"><sup>16</sup>O, without lattice</th><th colspan="3"><sup>16</sup>O, with lattice</th></tr><tr><th>B</th><th>B+EC</th><th>B+pycno</th><th>B</th><th>B+EC</th><th>B+pycno</th></tr></thead><tbody><tr><td>0</td><td>1.44</td><td>1.44</td><td>1.25</td><td>1.40</td><td>1.40</td><td>1.21</td></tr><tr><td>10<sup>33</sup></td><td>1.45</td><td>1.44</td><td>1.25</td><td>1.41</td><td>1.40</td><td>1.21</td></tr><tr><td>5 × 10<sup>33</sup></td><td>1.60</td><td>1.50</td><td>1.25</td><td>1.59</td><td>1.46</td><td>1.22</td></tr><tr><td>10<sup>34</sup></td><td>2.01</td><td>1.68</td><td>1.26</td><td>1.97</td><td>1.67</td><td>1.23</td></tr><tr><td>2 × 10<sup>34</sup></td><td>2.00</td><td>2.00</td><td>1.32</td><td>1.96</td><td>1.96</td><td>1.29</td></tr><tr><td>3 × 10<sup>34</sup></td><td>1.96</td><td>1.96</td><td>1.40</td><td>1.94</td><td>1.94</td><td>1.38</td></tr></tbody></table>

<table><caption>**Table 2.** Maximum masses (in M⊙) of strongly magnetized $^{12}$C WDs for various magnetic moments (in A m$^2$), limited by different mechanisms: B denotes the onset of torus-like shape due to magnetic field, EC stands for electron-capture instability and pycno is short for pycnonuclear fusion reactions.</caption><thead><tr><th rowspan="2">Magnetic moment<br>(A m<sup>2</sup>)</th><th colspan="3"><sup>12</sup>C, without lattice</th><th colspan="3"><sup>12</sup>C, with lattice</th></tr><tr><th>B</th><th>B+EC</th><th>B+pycno</th><th>B</th><th>B+EC</th><th>B+pycno</th></tr></thead><tbody><tr><td>0</td><td>1.44</td><td>1.44</td><td>1.41</td><td>1.41</td><td>1.41</td><td>1.38</td></tr><tr><td>10<sup>33</sup></td><td>1.45</td><td>1.44</td><td>1.41</td><td>1.41</td><td>1.41</td><td>1.38</td></tr><tr><td>5 × 10<sup>33</sup></td><td>1.60</td><td>1.53</td><td>1.43</td><td>1.60</td><td>1.50</td><td>1.40</td></tr><tr><td>10<sup>34</sup></td><td>2.01</td><td>1.85</td><td>1.48</td><td>2.00</td><td>1.86</td><td>1.47</td></tr><tr><td>2 × 10<sup>34</sup></td><td>2.01</td><td>2.00</td><td>1.74</td><td>1.99</td><td>1.99</td><td>1.75</td></tr><tr><td>3 × 10<sup>34</sup></td><td>2.03</td><td>2.03</td><td>2.03</td><td>1.96</td><td>1.96</td><td>1.96</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Gravitational mass $M_G$ versus circumferential radius $R_{circ}$ for magnetic $^{12}$C WDs along sequences of fixed magnetic dipole moments $D$ in Newtonian theory (bold lines) and general relativity (dashed lines).</figcaption></figure>

the high-density part of the EoS, which determines the maximum mass for the non-magnetic case, but for the magnetic case it is rather the low-density part of the EoS that determines the maximum mass. To understand this, we recall Fig. 3, where we plotted the masses as a function of density, from which it was evident that it is the low-density part of the curves that bends away from the non-magnetic mass–radius relation and are responsible for the large masses, and the high-density part remains unaffected. The fact that the effect of general relativity on the maximum masses is small (<2 per cent for poloidal fields) for the case of strongly magnetized WDs was demonstrated earlier by Das & Mukhopadhyay (2015) and Bera & Bhattacharya (2016) employing the $x$NS and LORENE codes.

## 7.2 Electron-capture instability

Chamel et al. (2013) pointed out that the onset of electron captures by nuclei (with the emission of a neutrino),

$$A_Z X + e^- \rightarrow A_{Z-1} Y + v_e, \quad (37)$$

may induce a local instability in magnetic WDs, thus limiting their maximum gravitational mass. In the absence of magnetic fields, the onset of electron captures occurs at mass density $\rho_\beta \simeq 4.16 \times 10^{10}$ g cm$^{-3}$ (pressure $P_\beta \simeq 6.99 \times 10^{28}$ dyn cm$^{-2}$) for $^{12}$C and $\rho_\beta \simeq 2.06 \times 10^{10}$ g cm$^{-3}$ ($P_\beta \simeq 2.73 \times 10^{28}$ dyn cm$^{-2}$) for $^{16}$O (Chamel & Fantina 2015). In the presence of a strong magnetic field, the threshold density and pressure are shifted to either higher or lower values, depending on the magnetic field strength (Chamel & Fantina 2015), an effect that has been neglected until recently. Otoniel et al. (2016) have considered two limiting cases: (i) the absence of magnetic field and (ii) the presence of a strongly quantizing magnetic field such that only the lowest Landau level is filled and $\rho_\beta = \rho_b$. In both cases, Otoniel et al. (2016) neglected the effects of electron–electron and electron–ion interactions on $\rho_\beta$, as in the study of Bera & Bhattacharya (2016). In the present work, we take into account the full dependence of the threshold density $\rho_\beta$ on the magnetic field strength including lattice corrections, as computed by Chamel & Fantina (2015). In this way, both the EoS and the electron-capture threshold are calculated fully consistently.

The resulting mass–radius relations and electron-capture thresholds for magnetic WDs along fixed magnetic moment sequences are shown in Fig. 11, considering stars made of either $^{12}$C or $^{16}$O. The filled and open dots along the sequences indicate the onset of electron-capture instability for $^{12}$C and $^{16}$O, respectively. The mass–radius relations for $^{16}$O differ only marginally from those obtained for $^{12}$C, and are not displayed in the figure for better readability. Our calculations confirm the suggestion of Chamel et al. (2013, 2014) that the electron-capture instability can be a limiting factor for determining the maximum WD mass and in turn the maximum magnetic
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Mass *M* versus radius *R*<sub>eq</sub> for magnetic WDs in Newtonian theory, composed of <sup>12</sup>C along sequences of fixed magnetic moments including lattice effects in the EoS. The filled dots mark the onset of electron capture on <sup>12</sup>C nuclei and the threshold electron-capture densities for daughter nucleus obtained from pycnonuclear reactions are marked with filled squares. The corresponding onset for <sup>16</sup>O WDs are shown by open symbols. The mass–radius relationship for <sup>16</sup>O WDs is not shown since the difference to <sup>12</sup>C is maximally a few per cent, hardly visible in the figure.</figcaption></figure>

field. For magnetic moments as large as $10^{34}$ A m$^2$, electron-capture instability (in the case with lattice) limits the maximum gravitational mass to 1.86 $M_{\odot}$ for $^{12}$C instead of 2.00 $M_{\odot}$. For $^{16}$O, the onset of the instability is reached at even lower values of the magnetic field, thus leading to a maximum gravitational mass of 1.67 $M_{\odot}$ instead of 1.97 $M_{\odot}$. This can be clearly understood from fig. 1 of Chamel & Fantina (2015), where the onset of electron-capture instability was shown to occur at lower densities for $^{16}$O than for $^{12}$C. For higher magnetic moments (see e.g. $D=2 \times 10^{34}$ and $3 \times 10^{34}$ A m$^2$ in Tables 1 and 2), the star becomes maximally distorted before the electron-capture threshold is reached. In this case, the maximum gravitational mass could thus be potentially higher than that shown in Fig. 11 (see Tables 1 and 2, and also Otoniel et al. 2016). Fig. 12 shows the masses of WDs, taking the example of $^{16}$O as an illustrative case, as a function of the core magnetic field strengths for different fixed magnetic dipole moments. Here, we show again the onset of electron-capture instabilities by open dots along the curves. It is evident from the figure that the higher the magnetic moment is, the larger is the increase of the mass with the core magnetic field strength. All in all, no stellar configurations with a poloidal magnetic field are found with a mass above 2 $M_{\odot}$, if the star is not allowed to take a toroidal shape.

## 7.3 Pycnonuclear instability

As first discussed by Chamel et al. (2013), pycnonuclear fusion reactions, whereby nuclei transform as

$$^{A}_{Z}X + ^{A}_{Z}X \rightarrow ^{2A}_{2Z}Y, \qquad (38)$$

could play an important role in limiting the masses of strongly magnetized WDs. However, the rates of these processes remain highly uncertain (Yakovlev et al. 2006). We can estimate the strongest impact on the maximum mass by considering that pycnonuclear reactions set in at the threshold density $\rho_{\beta}$ for the onset of electron capture by the daughter nuclei $^{2A}_{2Z}Y$. From table I of Chamel

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Central magnetic field $B_c$, in units of $b_{crit}$ defined in equation (6), versus WD mass *M* for different fixed magnetic moments *D* for Newtonian $^{16}$O WDs. The open dots show the onset of electron capture.</figcaption></figure>

& Fantina (2015), it is evident that the daughter nucleus $^{2A}_{2Z}Y$ is generally much more unstable than the parent nucleus $^{A}_{Z}X$ (e.g. the threshold density $\rho_{\beta}$ for $^{32}$S is $1.69 \times 10^8$ g cm$^{-3}$ while that of $^{16}$O is $2.06 \times 10^{10}$ g cm$^{-3}$, including lattice effects). In Fig. 11, we mark the threshold densities $\rho_{\beta}$ for the daughter nucleus $^{32}$S (from the fusion of $^{16}$O) and $^{24}$Mg (from the fusion of $^{12}$C), respectively, with open and filled squares along the curve. We find that for the non-magnetic case, the maximum mass decreases from 1.40 to 1.21 $M_{\odot}$ for $^{16}$O (Table 1). Excluding lattice interactions increases both masses by about 3 per cent. The effect is much less pronounced for $^{12}$C, where the threshold density is higher. Here, the maximum mass in the non-magnetized case becomes 1.38 $M_{\odot}$ (Table 2). This is in agreement with the results quoted in Chamel & Fantina (2015). For large values of magnetic moment, the reduction of the maximum mass of WDs is even more spectacular. For example, for a magnetic moment of $10^{34}$ A m$^2$, the maximum mass decreases from 1.97 to 1.23 $M_{\odot}$ for $^{16}$O WDs and from 2.00 to 1.47 $M_{\odot}$ for $^{12}$C WDs. Because pycnonuclear reactions may actually occur at densities above the electron-capture threshold density for the daughter nuclei, the values we found for the maximum mass represent lower bounds. One must note that, for magnetic moment $3 \times 10^{34}$ A m$^2$ for $^{12}$C WDs, the star becomes distorted to a toroidal shape before the pycnonuclear threshold can be reached, hence the reduction in maximum mass cannot be calculated within our numerical framework and is not reflected in the quoted values in Table 2.

## 8 GRAVITATIONAL WAVE EMISSION

It was conjectured by Heyl (2000) that rotating magnetic WDs could be important sources of GWs detectable by the proposed LISA mission (Laser Interferometer Space Antenna; Amaro-Seoane et al. 2017).<fnref n="4" /> Heyl estimated the GW emission for a number of observed rotating magnetic WDs using simple approximations for the reduced quadrupole moment. Using the analytic prescription of Bonazzola & Gourgoulhon (1996), where the distortions due to the magnetic field and due

<footnote n="4">⁴ https://elisascience.org.</footnote>
to rotation are assumed to be decoupled, we calculated the GW amplitudes for the potential sources listed in table 1 of Heyl (2000). We numerically estimated the quadrupole moment in a self-consistent way within our improved microscopic model, assuming a gravitational mass of $0.6 M_{\odot}$ and an internal magnetic field of $10^{11}$ G. We obtained a lower value of $0.58 \times 10^{46}$ g cm$^2$ for quadrupole moment, close to the value $10^{47}$ g cm$^2$ obtained in Heyl (2000). However, although the discussion in Heyl (2000) concluded that the amplitude was within the detectable range of the former *LISA* interferometer project, our calculations of GW amplitudes for the rotating magnetic WDs estimated here are not likely to be detectable by the current project *LISA*, as they are not within the estimated sensitivity range.

## 9 CONCLUSIONS

Using the formalism developed by Chatterjee et al. (2015), we have studied the equilibrium structure of WDs endowed with a strong poloidal magnetic field focusing on the determination of the maximum mass. In our approach, the coupled equilibrium equations for magnetic and gravitational fields are solved taking consistently into account stellar deformations due to rotations and/or anisotropies introduced by the magnetic field.

For our investigation, we have employed both general relativity and Newtonian theory for the description of gravity, together with the EoS of a degenerate electron Fermi gas interacting with a pure ionic crystal lattice made of $^{12}$C or $^{16}$O. The magnetic field dependence of the EoS induced by Landau quantization of electron motion has been fully included, as well as the resulting magnetization of stellar matter. However, our numerical calculations have demonstrated that neither of these effects does significantly alter the structure of magnetized WDs, and they can thus be neglected. Still, the magnetic field can have a large impact on the stellar mass and radius. In particular, stellar configurations more massive than the standard Chandrasekhar limit can be obtained for magnetic field strengths higher than the critical value $b_{crit} \sim 4.4 \times 10^{13}$ G, as found in previous works. We have also examined more closely the role of electron-electron and electron-ion interactions, which have been often neglected in previous studies of strongly magnetized WDs. Because these interactions are attractive, the maximum mass is reduced by a few per cent. On the other hand, the stellar radius can be increased by up to ~50 per cent for non-magnetic WDs with a mass of $1.4 M_{\odot}$, and therefore these interactions should be included in any realistic model of WDs.

To explore whether strongly magnetized WDs are globally stable, we have performed calculations of sequences of equilibrium stellar configurations along fixed values of the magnetic dipole moment. We have found that the presence of a strong magnetic field can lead to large deformations of the star. For high enough values of the magnetic dipole moment, the stellar shape thus becomes toroidal. This extreme configuration is reached for gravitational masses below $2 M_{\odot}$ and a surface magnetic field of the order of $10^{13}-10^{14}$ G. Note that, with our models of poloidal magnetic field the core magnetic strength is about one order of magnitude higher than the surface one. More strongly magnetized stellar configurations cannot be treated within our numerical framework, as the star gets a torus-like shape (magnetic pressure becomes larger than fluid pressure in the centre). We have also investigated the role of rotation on the stellar structure, and we have found that the maximum mass of magnetic WDs is increased by a few per cent at most, the effects being the largest in the absence of magnetic fields.

Although torus-shaped strongly magnetized WDs with a mass above $2 M_{\odot}$ could potentially exist beyond the distorted configurations computed in this work, various instabilities can arise. First of

all, it is well known that general relativity can limit the maximum mass of non-magnetic WDs. We have thus computed the structure of strongly magnetized WDs in full general relativity and have found that for weakly magnetized WDs, general relativistic effects reduce the maximum mass, as already discussed by De-Hua et al (2014), Coelho et al. (2014) and Bera & Bhattacharya (2016). On the other hand, general relativity hardly plays any role in limiting the masses of strongly magnetized WDs, as anticipated by Kundu & Mukhopadhyay (2012). More importantly, the stability of strongly magnetized WDs is found to be mainly limited by the onset of electron capture and pycnonuclear reactions in the stellar core, as argued by Chamel et al. (2013). We have determined the threshold densities and pressures consistently with the EoS, including for the first time the effects due to both electrostatic interactions and the magnetic field. Given the high uncertainties on the pycnonuclear fusion reaction rates, we have estimated the maximum possible reduction of the WD mass by assuming that these processes set in at the same threshold density and pressure as electron captures by the daughter nuclei. These reactions lead to a drastic decrease in the maximum WD mass to values even below the Chandrasekhar limit and is more pronounced for $^{16}$O WDs than for $^{12}$C ones since the corresponding threshold densities are lower for $^{16}$O than for $^{12}$C. Our numerical results about mass limits are summarized in Tables 1 and 2. Additionally, we have estimated GW amplitudes emitted from rotating magnetized WDs, for which the magnetic and rotation axes are misaligned. We come to the conclusion that, from the currently observed WDs, none could be detected by the future space-based GW detector *LISA*.

Extremely magnetized WDs can thus in principle reach masses higher than $2 M_{\odot}$ and remain stable with respect to gravitational and electron-capture instabilities; the consequence is that they would have torus-like shapes. On the other hand, pycnonuclear instabilities could severely limit the existence of such stars, but one must keep in mind the high uncertainties associated with these reactions. In summary, the possibility of super-Chandrasekhar strongly magnetized WDs cannot be totally excluded from current theoretical considerations, but important issues still need to be addressed before any firm conclusions on their existence could be drawn. First, the dynamical stability of these extremely high magnetic fields in such cold dense crystallized stellar environment remains to be proved. Moreover, the magnetic field strengths expected at the surface of such strongly magnetized WDs appear to be four orders of magnitude higher than the upper limit of about $10^9$ G set by currently observed magnetic WDs. Finally, no realistic and quantitative astrophysical formation scenarios have been so far proposed to explain the origin of such strongly magnetized WDs.

## ACKNOWLEDGEMENTS

We warmly thank M. Mouchet and E. Gourgoulhon for fruitful discussions. DC would like to thank the Observatoire de Paris in Meudon, particularly the computing department, for their constant support. This work has been funded by the SN2NS project ANR-10-BLAN-0503, the 'Gravitation et physique fondamentale' action of the Observatoire de Paris, Fonds de la Recherche Scientifique - FNRS (Belgium) and the COST action MP1304 'New-CompStar'.

## REFERENCES

Adam D., 1986, A&A, 160, 95
Amaro-Seoane P. et al., 2017, preprint (arXiv:1702.00786)
Audi G., Wang M., Wapstra A. H., Kondev F. G., MacCormick M., Xu X., Pfeiffer B., 2012, Chin. Phys. C, 36, 1287

Baiko D. A., 2009, Phys. Rev. E, 80, 046405

Baiko D. A., Potekhin A. Y., Yakovlev D. G., 2001, Phys. Rev. E, 64, 057402

Belyaev V. B., Ricci P., Šimkovic F., Adam J., Truhlík E., 2015, Nucl. Phys. A, 937, 17

Bera P., Bhattacharya D., 2014, MNRAS, 445, 3951

Bera P., Bhattacharya D., 2016, MNRAS, 456, 3375

Blackett P. M. S., 1947, Nature, 159, 658

Bocquet M., Bonazzola S., Gourgoulhon E., Novak J., 1995, A&A, 301, 757

Bonazzola S., Gourgoulhon E., 1994, Class. Quantum Gravity, 11, 1775

Bonazzola S., Gourgoulhon E., 1996, A&A, 312, 675

Bonazzola S., Gourgoulhon E., Salgado M., Marck J. A., 1993, A&A, 278, 421

Born M., 1940, Proc. Cambridge Philos. Soc., 36, 160

Boshkayev K., Rueda J. A., Ruffini R., Siutsou I., 2015, Thirteenth Marcel Grossmann Meeting: On Recent Developments in Theoretical and Experimental General Relativity, Astrophysics and Relativistic Field Theories. World Scientific, Singapore, p. 2468

Braithwaite J., Spruit H. C., 2004, Nature, 431, 819

Bucciantini N., Pili A. G., Del Zanna L., 2014, ASP Conf. Ser., 488, 211

Bucciantini N., Pili A. G., Del Zanna L., 2016, J. Phys. Conf. Ser., 719, 012004

Cardall C. Y., Prakash M., Lattimer J. M., 2001, ApJ, 554, 322

Chamel N., Fantina A. F., 2015, Phys. Rev. D, 92, 023008

Chamel N. et al., 2012, Phys. Rev. C, 86, 055804

Chamel N., Fantina A. F., Davis P. J., 2013, Phys. Rev. D, 88, 081301

Chamel N., Molter E., Fantina A. F., Peña Arteaga D., 2014, Phys. Rev. D, 90, 043002

Chandrasekhar S., 1931, ApJ, 74, 81

Chandrasekhar S., 1935, MNRAS, 95, 207

Chatterjee D., Elghozi T., Novak J., Oertel M., 2015, MNRAS, 447, 3785

Cheoun M.-K., Deliduman C., Gündör C., Keleş V., Ryu C. Y., Kajino T., Mathews G. J., 2013, J. Cosmol. Astropart. Phys., 10, 021

Coelho J. G., Marinho R. M., Malheiro M., Negreiros R., Cáceres D. L., Rueda J. A., Ruffini R., 2014, ApJ, 794, 86

Das U., Mukhopadhyay B., 2012a, Int. J. Mod. Phys. D, 21, 1242001

Das U., Mukhopadhyay B., 2012b, Phys. Rev. D, 86, 042001

Das U., Mukhopadhyay B., 2013a, Int. J. Mod. Phys. D, 22, 1342004

Das U., Mukhopadhyay B., 2013b, Phys. Rev. Lett., 110, 071102

Das U., Mukhopadhyay B., 2015, J. Cosmol. Astropart. Phys., 5, 16

Das U., Mukhopadhyay B., Rao A. R., 2013, ApJ, 767, L14

De-Hua W., He-Lei L., Xiang-Dong Z., 2014, Chin. Phys. B, 23, 089501

Federbush P., Luo T., Smoller J., 2015, Arch. Ration. Mech. Anal., 215, 611

Ferrario L., de Martino D., Gänsicke B. T., 2015, Space Sci. Rev., 191, 111

Franzon B., Schramm S., 2015, Phys. Rev. D, 92, 083006

Fujisawa K., Yoshida S., Eriguchi Y., 2012, MNRAS, 422, 434

Gourgoulhon E., 2010, Lectures given at the conference CompStar 2010: School and Workshop on Computational Tools for Compact Star Astrophysics, preprint (arXiv:1003.5015)

Gourgoulhon E., Bonazzola S., 1994, Class. Quantum Gravity, 11, 443

Gourgoulhon E., Grandclément P., Marck J.-A., Novak J., Taniguchi K., 2016, Astrophysics Source Code Library, record ascl:1608.018

Grandclément P., Novak J., 2009, Living Rev. Relativ., 12, 1

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

Haensel P., Potekhin A. Y., Yakovlev D. G., 2007, Neutron Stars 1: Equation of state and structure, Springer, Berlin

Hamada T., Salpeter E. E., 1961, ApJ, 134, 683

Herrera L., Barreto W., 2013, Phys. Rev. D, 87, 087303

Herrera L., Di Prisco A., Barreto W., Ospino J., 2014, Gen. Relativ. Gravit., 46, 1827

Heyl J. S., 2000, MNRAS, 317, 310

Hillebrandt W., Kromer M., Röpke F. K., Ruiter A. J., 2013, Frontiers Phys., 8, 116

Howell D. A. et al., 2006, Nature, 443, 308

Ibáñez Cabanell J. M., 1983, Rev. Mex. Astron. Astrofís., 8, 7

Ibáñez Cabanell J. M., 1984, A&A, 135, 382

Jordan S., 2009, in Strassmeir K. G., Kosovichev A. G., Beckman J. E., eds, Proc. IAU Symp. 162, Cosmic Magnetic Fields: From Planets, to Stars and Galaxies. Cambridge Univ. Press, Cambridge, p. 369

Kemp J. C., 1970, ApJ, 162, 169

Kepler S. O. et al., 2013, MNRAS, 429, 2934

Kozhberov A. A., 2016, Ap&SS, 361, 256

Kundu A., Mukhopadhyay B., 2012, Mod. Phys. Lett. A, 27, 1250084

Lai D., Shapiro S. L., 1991, ApJ, 383, 745

Lunney D., Pearson J. M., Thibault C., 2003, Rev. Mod. Phys., 75, 1021

Maoz D., Mannucci F., Nelemans G., 2014, ARA&A, 52, 107

Marshak R. E., 1940, ApJ, 92, 321

Novak J., Marcq E., 2003, Class. Quantum Gravity, 20, 3051

Ostriker J. P., Hartwick F. D. A., 1968, ApJ, 153, 797

Otoniel E., Franzon B., Malheiro M., Schramm S., Weber F., 2016, preprint (arXiv:1609.05994)

Peña Arteaga D., Grasso M., Khan E., Ring P., 2011, Phys. Rev. C, 84, 045806

Potekhin A. Y., Chabrier G., 2013, A&A, 550, A43

Rotondo M., Rueda J. A., Ruffini R., Xue She-Sheng, 2011, Phys. Rev. D, 84, 084007

Salpeter E. E., 1961, ApJ, 134, 669

Scalzo R. A. et al., 2010, ApJ, 713, 1073

Shapiro S. L., 2000, ApJ, 544, 397

Shapiro S. L., Teukolsky S. A., 1983, Black Holes, White Dwarfs, and Neutron Stars. Wiley, New York

Shul'man G. A., 1976, SvA, 20, 689.

Sorkin R., 1982, ApJ, 257, 847

Stein M., Maruhn J., Sedrakian A., Reinhard P.-G., 2016, Phys. Rev. C, 94, 035802

Subramanian S., Mukhopadhyay B., 2015, MNRAS, 454, 752

Suh I.-S., Mathews. G. J., 2000, ApJ, 530, 949

Van Vleck J. H., 1932, The Theory of Electric and Magnetic Susceptibilities. Oxford Univ. Press, London

Wickramasinghe D. T., Ferrario L., 2000, PASP, 112, 873

Yakovlev D. G., Gasques L. R., Afanasjev A. V., Beard M., Wiescher M., 2006, Phys. Rev. C, 74, 035803

## APPENDIX A: GENERAL EXPRESSION FOR THE MAGNETIZATION

We report here for completeness and future reference the general expressions for the magnetization and its first derivative, which are needed with the EoS to compute the WD structure. Starting from the definition of $M \equiv \partial P/\partial b|_{\mu}$, equation (19), we first observe that constant chemical potential implies

$$d\mu = 0 \Rightarrow d\mu_e = -d\mu_L, \quad (A1)$$

where we have used equations (21) and (22). Writing explicitly the dependence of the pressure on the magnetic field $b$ and the electron chemical potential $\mu_e$ as $P = P(b, \mu_e)$, we obtain

$$M \equiv \left. \frac{\partial P}{\partial b} \right|_{\mu} = \left. \frac{\partial P}{\partial b} \right|_{\mu_e} - \left. \frac{\partial P}{\partial \mu_e} \right|_{b} \left. \frac{\partial \mu_L}{\partial b} \right|_{\mu_e} \left[ 1 + \left. \frac{\partial \mu_L}{\partial \mu_e} \right|_{b} \right]^{-1}, \quad (A2)$$

with

$$\frac{\partial P}{\partial b} = \frac{\partial P_L}{\partial b} + \frac{\partial P_e}{\partial b}, \quad (A3)$$

$$\frac{\partial P}{\partial \mu_e} = \frac{\partial P_L}{\partial \mu_e} + \frac{\partial P_e}{\partial \mu_e}. \quad (A4)$$
Using the definition of $b_{crit}$ and $\gamma_e$ from equations (6) and (9), we also note that

$$ \frac{\partial}{\partial b} = \frac{1}{b_{crit}} \frac{\partial}{\partial b_{*}}, \qquad (A5) $$

$$ \frac{\partial}{\partial \mu_e} = \frac{1}{m_e c^2} \frac{\partial}{\partial \gamma_e}. \qquad (A6) $$

We thus find for the first derivative of the pressure (lattice and electron contribution) with respect to $b$:

$$ \left. \frac{\partial P_L}{\partial b} \right|_{\mu_e} = \frac{4 C e^2}{9 b_{crit}} Z^{2/3} n_e^{1/3} \left. \frac{\partial n_e}{\partial b_*} \right|_{\gamma_e}, \qquad (A7) $$

with

$$ \left. \frac{\partial n_e}{\partial b_*} \right|_{\gamma_e} = \frac{n_e}{b_*} + \frac{2b_*}{(2\pi)^2 \lambda_e^3} \sum_{\nu_L=0}^{\nu_L^{\max}} g_\nu \frac{\partial x_e}{\partial b_*}, \qquad (A8) $$

$$ \frac{\partial x_e}{\partial b_*} = -\frac{\nu_L}{x_e}; \qquad (A9) $$

$$ \left. \frac{\partial P_e}{\partial b} \right|_{\mu_e} = \frac{P_e}{b} - \frac{b_* m_e c^2}{(2\pi)^2 \lambda_e^3 b_{crit}} \times \sum_{\nu_L=0}^{\nu_L^{\max}} 2g_\nu \nu_L \ln \left( \frac{\gamma_e + x_e}{\sqrt{1 + 2\nu_L b_*}} \right). \qquad (A10) $$

For the first derivative of the pressure (lattice and electron contribution) with respect to $\mu_e$, we find

$$ \left. \frac{\partial P_L}{\partial \mu_e} \right|_b = \frac{4 C e^2}{9 m_e c^2} Z^{2/3} n_e^{1/3} \left. \frac{\partial n_e}{\partial \gamma_e} \right|_{b_*}, \qquad (A11) $$

with

$$ \left. \frac{\partial n_e}{\partial \gamma_e} \right|_{b_*} = \frac{2b_*}{(2\pi)^2 \lambda_e^3} \sum_{\nu_L=0}^{\nu_L^{\max}} g_\nu \frac{\partial x_e}{\partial \gamma_e}, \qquad (A12) $$

$$ \frac{\partial x_e}{\partial \gamma_e} = \frac{\gamma_e}{x_e}; \qquad (A13) $$

$$ \left. \frac{\partial P_e}{\partial \mu_e} \right|_b = \frac{b_*}{(2\pi)^2 \lambda_e^3} \sum_{\nu_L=0}^{\nu_L^{\max}} 2g_\nu x_e = n_e. \qquad (A14) $$

The first derivatives of the chemical potential can be written as

$$ \left. \frac{\partial \mu_L}{\partial b} \right|_{\mu_e} = \frac{4 C e^2}{9 b_{crit}} Z^{2/3} n_e^{-2/3} \left. \frac{\partial n_e}{\partial b_*} \right|_{\gamma_e}, \qquad (A15) $$

using equation (A8), and

$$ \left. \frac{\partial \mu_L}{\partial \mu_e} \right|_b = \frac{4 C e^2}{9 m_e c^2} Z^{2/3} n_e^{-2/3} \left. \frac{\partial n_e}{\partial \gamma_e} \right|_{b_*}, \qquad (A16) $$

using equation (A12).

To calculate the derivative of $M$,

$$ \frac{\partial M}{\partial \mu} = \frac{\partial^2 P}{\partial \mu \partial b} = \frac{\partial M}{\partial \mu_e} \left. \right|_b \left( \frac{\partial \mu}{\partial \mu_e} \right) \left. \right|_b^{-1}, \qquad (A17) $$

we write explicitly the terms appearing in equation (A17) as

$$ \left. \frac{\partial M}{\partial \mu_e} \right|_b = \frac{1}{m_e c^2} \left\{ \left. \frac{\partial}{\partial \gamma_e} \frac{\partial P}{\partial b} \right|_{\mu_e} - \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial P}{\partial \mu_e} \right|_{b_*} \right) \frac{\partial \mu_L}{\partial b} \left|_{\mu_e} \left( 1 + \frac{\partial \mu_L}{\partial \mu_e} \right) \right|_{b_*}^{-1} - \frac{\partial P}{\partial \mu_e} \left|_{b_*} \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial \mu_L}{\partial b} \right|_{\mu_e} \right) \right) \left( 1 + \frac{\partial \mu_L}{\partial \mu_e} \right) \right|_{b_*}^{-1} + \frac{\partial P}{\partial \mu_e} \left|_{b_*} \frac{\partial \mu_L}{\partial b} \right|_{\mu_e} \left( 1 + \frac{\partial \mu_L}{\partial \mu_e} \right) \right|_{b_*}^{-2} \times \left. \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial \mu_L}{\partial \mu_e} \right|_{b_*} \right) \right\}, \qquad (A18) $$

and

$$ \left. \frac{\partial \mu}{\partial \mu_e} \right|_b = y_e \left( 1 + \frac{\partial \mu_L}{\partial \mu_e} \right) \left|_b \right. , \qquad (A19) $$

where each contribution of the pressure derivative can be calculated according to equations (A3) and (A4). We thus obtain, for the second derivative of the pressure (electron and lattice contribution):

$$ \left. \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial P_L}{\partial b} \right|_{\mu_e} \right) \right|_{b_*} = \frac{4 C e^2}{9 b_{crit}} Z^{2/3} \left[ \frac{1}{3} n_e^{-2/3} \left. \frac{\partial n_e}{\partial \gamma_e} \right|_{b_*} \frac{\partial n_e}{\partial b_*} \right]_{\gamma_e} + n_e^{1/3} \left. \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial n_e}{\partial b_*} \right|_{b_*} \right) \right] \qquad (A20) $$

using equations (A8) and (A12), and

$$ \left. \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial n_e}{\partial b_*} \right|_{\mu_e} \right) \right|_{b_*} = \frac{2}{(2\pi)^2 \lambda_e^3} \sum_{\nu_L=0}^{\nu_L^{\max}} g_\nu \frac{\partial x_e}{\partial \gamma_e} + \frac{2b_*}{(2\pi)^2 \lambda_e^3} \sum_{\nu_L=0}^{\nu_L^{\max}} g_\nu \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial x_e}{\partial b_*} \right) \right|_{b_*}, \qquad (A21) $$

with $\partial x_e/\partial \gamma_e$ from equation (A13) and

$$ \frac{\partial}{\partial \gamma_e} \left( \frac{\partial x_e}{\partial b_*} \right) = \frac{\nu_L}{x_e^3}; \qquad (A22) $$

$$ \left. \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial P_e}{\partial b} \right|_{\mu_e} \right) \right|_{b_*} = \frac{n_e m_e c^2}{b} - \frac{b_* m_e c^2}{(2\pi)^2 \lambda_e^3 b_{crit}} \sum_{\nu_L=0}^{\nu_L^{\max}} 2 g_\nu \nu_L; \qquad (A23) $$

$$ \left. \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial P_L}{\partial \mu_e} \right|_{b_*} \right) \right|_{b_*} = \frac{4 C e^2}{9 m_e c^2} Z^{2/3} \left[ \frac{1}{3} n_e^{-2/3} \left( \left. \frac{\partial n_e}{\partial \gamma_e} \right|_{b_*} \right)^2 + n_e^{1/3} \left. \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial n_e}{\partial \gamma_e} \right|_{b_*} \right) \right], \qquad (A24) $$

using equation (A12) and

$$ \left. \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial n_e}{\partial \gamma_e} \right|_{b_*} \right) \right|_{b_*} = \frac{2b_*}{(2\pi)^2 \lambda_e^3} \sum_{\nu_L=0}^{\nu_L^{\max}} g_\nu \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial x_e}{\partial \gamma_e} \right) \right), \qquad (A25) $$

with:

$$ \frac{\partial}{\partial \gamma_e} \left( \frac{\partial x_e}{\partial \gamma_e} \right) = - \frac{1 + 2\nu_L b_*}{x_e^3}; \qquad (A26) $$
$$ \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial P_e}{\partial \mu_e} \right|_{b_*} \right)_{b_*} = \frac{2b_*}{(2\pi)^2 \lambda_e^3} \sum_{\nu_L=0}^{\nu_L^{\max}} \frac{g_\nu \gamma_e}{x_e}. \quad (A27) $$

For the second derivatives of the chemical potential, we have

$$ \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial \mu_L}{\partial b} \right|_{\mu_e} \right)_{b_*} = \frac{4 C e^2}{9 b_{crit}} Z^{2/3} \times \left[ -\frac{2}{3} n_e^{-5/3} \left( \left. \frac{\partial n_e}{\partial \gamma_e} \right|_{b_*} \right) \left( \left. \frac{\partial n_e}{\partial b_*} \right|_{\mu_e} \right) + n_e^{-2/3} \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial n_e}{\partial b_*} \right|_{\mu_e} \right)_{b_*} \right], \quad (A28) $$

where we have used equations (A12) and (A21), and

$$ \frac{\partial}{\partial \gamma_e} \left( \left. \left( \frac{\partial \mu_L}{\partial \mu_e} \right)_{b_*} \right|_{b_*} \right) = \frac{4 C e^2}{9 m_e c^2} Z^{2/3} \times \left[ -\frac{2}{3} n_e^{-5/3} \left( \left. \frac{\partial n_e}{\partial \gamma_e} \right|_{b_*} \right)^2 + n_e^{-2/3} \frac{\partial}{\partial \gamma_e} \left( \left. \frac{\partial n_e}{\partial \gamma_e} \right|_{b_*} \right)_{b_*} \right], \quad (A29) $$

where we made use of equations (A12) and (A25).

## APPENDIX B: DERIVATION OF THE FIRST INTEGRAL IN GENERAL RELATIVITY

The constraint for the conservation of energy and momentum in equilibrium can be expressed as

$$ \nabla_{\mu} T^{\mu\nu} = 0. \quad (B1) $$

Inserting the expression for the energy-momentum tensor as in equation (23) into the above equation yields (Chatterjee et al. 2015)

$$ (\mathcal{E} + P) \left( \frac{1}{\mathcal{E} + P} \frac{\partial P}{\partial x^i} + \frac{\partial \nu}{\partial x^i} - \frac{\partial \ln \Gamma}{\partial x^i} \right) - F_{i\rho} j_{\text{free}}^{\rho} - \frac{\chi}{2\mu_0} F_{\mu\nu} \nabla_i \Gamma^{\mu\nu} = 0, \quad (B2) $$

where one can identify the first part as the perfect-fluid contribution to the energy-momentum tensor and the Lorentz force terms arising

from free currents and the magnetization. Defining the enthalpy as a function of both baryon density and magnetic field norm $b = \sqrt{b^\mu b_\mu}$ in the comoving frame:

$$ h = h(n, b) = \frac{\mathcal{E} + P}{m_B n c^2}, \quad (B3) $$

where $m_B$ is the mean mass of a baryon. One can rewrite the first term in the parentheses of equation (B2) in terms of the enthalpy as

$$ \frac{\partial P}{\partial x^i} = (\mathcal{E} + P) \frac{\partial \ln h}{\partial x^i} + \mathcal{M} \frac{\partial b}{\partial x^i}. \quad (B4) $$

The magnetic field in the comoving frame $b^\mu$ is defined by

$$ F_{\mu\nu} = \epsilon_{\alpha\beta\mu\nu} u^\beta b^\alpha \quad (B5) $$

and

$$ \mathcal{M}_{\mu\nu} = \epsilon_{\alpha\beta\mu\nu} u^\beta \mathcal{M}^\alpha. \quad (B6) $$

Further, the last term in equation (B2) can be written in terms of the magnetic field $b^\mu$ in the comoving frame as

$$ \begin{aligned} \frac{\chi}{2\mu_0} F_{\mu\nu} \nabla_i F^{\mu\nu} &= \frac{\chi}{\mu_0} b_\mu \nabla_i b^\mu - b_\mu b^\mu u_\nu \nabla_i u^\nu \\ &= b \nabla_i b = \mathcal{M} \frac{\partial b}{\partial x^i}. \end{aligned} \quad (B7) $$

Thus, inserting equations (B4) and (B7) into equation (B2), one gets an expression without explicit dependence on magnetization

$$ \frac{\partial \ln h}{\partial x^i} + \frac{\partial \nu}{\partial x^i} - \frac{\partial \ln \Gamma}{\partial x^i} - \frac{F_{i\rho} j_{\text{free}}^{\rho}}{\mathcal{E} + P} = 0, \quad (B8) $$

which has the same form as in absence of magnetization (Bonazzola et al. 1993; Bocquet et al. 1995). Finally, as in Chatterjee et al. (2015), we introduce the current function, whose primitive is noted $\Phi(r, \theta)$ and defined by equation (30). This enables us to get a first integral of motion as given in equation (29).

This paper has been typeset from a TeX/Latex file prepared by the author.