Strongly Magnetized White Dwarfs and Their Instability Due to Nuclear Processes

It is generally accepted that stars with masses below 10 solar masses end their evolutions as white dwarfs (WDs) (Weber 1999; Shapiro & Teukolsky 1983; Glendenning 2012). With a typical composition mostly made of carbon, oxygen, or helium, WDs possess central densities up to ∼10¹¹ g cm⁻³. They can be fast-rotating (Kepler et al. 2007; Kleinman et al. 2013; Ferrario et al. 2015) and strongly magnetized (Kepler et al. 2013; García-Berro et al. 2016). The observed surface magnetic fields range from 10⁶ G to 10⁹; see, e.g., G Putney (1995), Schmidt & Smith (1995), Reimers et al. (1996), Kemp et al. (1970), Angel (1978), and Terada et al. (2008). The internal magnetic fields of white dwarfs are not known, however, they are expected to be larger than their surface magnetic fields. This is due to the fact that in ideal magnetohydrodynamics, the magnetic field B is "frozen-in" with the fluid and B ∝ ρ, with ρ being the local mass density (see, e.g., Mestel 2012; Landau et al. 1958). A simple estimate of the internal magnetic field strength follows from the virial theorem by equating the magnetic field energy with the gravitational binding energy, which leads to an upper limit for the surface magnetic fields inside WDs of about ∼10¹³ G. On the other hand, analytical and numerical calculations, in the framework of both Newtonian gravity and general relativity, show that WDs may have internal magnetic fields as large as 10^12 − 16 G (see, e.g., Angel 1978; Shapiro & Teukolsky 1983; Das & Mukhopadhyay 2014a; Franzon & Schramm 2015a, 2017; Das & Mukhopadhyay 2015; Bhattacharya et al. 2018; Chatterjee et al. 2017).

The relationship between the gravitational stellar mass, M, and the radius, R, of non-magnetized white dwarfs was first determined by Chandrasekhar (1939). Recently, mass–radius relationships of magnetic white dwarfs have been explored in different studies (see, e.g., Suh & Mathews 2000; Bera & Bhattacharya 2014; Franzon & Schramm 2015a; Chatterjee et al. 2017). These studies show that the masses of white dwarfs increase in the presence of strong magnetic fields. This is due to the Lorentz force, which acts against gravity, thereby supporting stars with higher masses.

Based on recent observations of several superluminous type Ia supernovae (SN 2006gz, SN 2007if, SN 2009dc, SN 2003fg) in Howell et al. (2006), Hicken et al. (2007), Kepler et al. (2007), Yamanaka et al. (2009), Scalzo et al. (2010), Silverman et al. (2011), and Taubenberger et al. (2011), it has been suggested that the progenitor masses of such supernovae significantly exceed the Chandrasekhar mass limit of M_Ch ∼ 1. 4 M_⊙ (Ilkov & Soker 2012). Super-heavy progenitors were studied as a result of mergers of two massive white dwarfs (Ji et al. 2013; Moll et al. 2014; van Rossum et al. 2016). Alternatively, the authors Liu et al. (2014) and Carvalho et al. (2018a) obtained super-Chandrasekhar white dwarfs for electrically charged stars and suggested that they may be considered possible progenitors of supernovae, although there is not yet strong evidence for the existence of high electric fields in white dwarfs. In addition, super-Chandrasekhar white dwarfs were investigated in the presence of strong magnetic fields by Das & Mukhopadhyay (2014b). Adam (1986) and Ostriker & Hartwick (1968) calculated WDs models with magnetic fields in the framework of Newtonian physics. A recent study by Subramanian & Mukhopadhyay (2015) of differentially rotating, magnetized white dwarfs has shown that differential rotation might increase the maximum mass of magnetized white dwarfs up to 3.1 M_⊙. Also, as shown by Bera & Bhattacharya (2016), purely toroidal magnetic field components can increase the masses of white dwarfs up to 5 M_⊙.

According to Das & Mukhopadhyay (2012), the effects of an extremely large and uniform magnetic field on the equation of state (EoS) of a white dwarf could increase its critical mass up to 2.58 M_⊙. This mass limit is reached for extremely large magnetic fields of ∼10¹⁸ G. Nevertheless, as already discussed in the literature (Chamel et al. 2013; Coelho et al. 2014; Liccardo et al. 2018), the breaking of spherical symmetry due to magnetic fields and microphysical effects, such as electron capture reactions and pycnonuclear reactions, can severely limit the magnetic fields inside white dwarfs.

Franzon & Schramm (2015a) computed the mass–radius relationship of highly magnetized white dwarfs using a pure degenerate electron Fermi gas. However, according to Salpeter (1961), many-body corrections modify the EoS and therefore the mass–radius relationship of white dwarfs. The purpose of our paper is twofold. First, we model the EoS of white dwarfs, taking into account not only the electron Fermi gas contribution, but also the contribution from electron-ion interactions (Chamel & Fantina 2015). Second, we perform a stability analysis of the matter in the cores of white dwarfs against electron capture and pycnonuclear fusion reactions. The Landau energy levels of electrons are modified by relativistic effects if the magnetic field strength is higher than the critical QED magnetic field strength of B_cr = 4.4 × 10¹³ G. However, as already shown by Bera & Bhattacharya (2014), the global properties of white dwarfs, such as masses and the radii, are nearly independent of Landau quantization. For this reason, we do not take into account magnetic field effects in the EoS when calculating the global properties of WDs. Aside from the importance of magnetic white dwarfs for superluminous type Ia supernovae, they can also constitute rotation-powered pulsars (Marsh et al. 2016; Mukhopadhyay et al. 2017).

Our paper is organized as follows. In Section 2, we discuss the stellar interior of white dwarfs and the details of the EoS used in our study. In Section 3, we briefly discuss the equations that are being solved numerically to obtain the structure of stationary magnetized white dwarfs. In Section 4 we introduce the Einstein–Maxwell tensor and the metric tensor used to solve Einstein's field equations of general relativity. Our results are presented in Section 5 and summarized in Section 6.

White dwarf matter is largely composed of atomic nuclei immersed in a fully degenerate electron gas (Hamada & Salpeter 1961; Salpeter 1961). In this work, the EoS of such matter is determined for more recent atomic mass data (Audi et al. 2012; Wang et al. 2012). Magnetic field effects on the EoS of white dwarf matter are ignored. The model adopted for the nuclear lattice in the outer crust of a neutron star is described, for instance, in Shapiro & Teukolsky (1983). The pressure in the cores of white dwarfs is produced by degenerate electrons and the ionic lattice, leading to the total pressure given by

$$\begin{eqnarray}&&P={P}_{e}+{P}_{L}(Z),\end{eqnarray} \tag{ 1 }$$

where P_e denotes the electron pressure (determined in Salpeter 1961), and P_L(Z) is the lattice pressure of ions, where Z is the proton number of the ions. The lattice pressure can be written in terms of the energy density of the ionic lattice (see Pearson et al. (2011)),

$$\begin{eqnarray}&&{P}_{L}(Z)=\displaystyle \frac{1}{3}{{ \mathcal E }}_{L}.\end{eqnarray} \tag{ 2 }$$

In this work we consider Z = 12, which implies that the white dwarf is composed of carbon ions.

The lattice pressure of ions arranged in a regular body-centered-cubic (bcc) crystal does not depend on the magnetic field, apart from a small contribution due to the quantum zero-point motion of ions. In this case, the lattice energy density reads as (Chamel et al. 2013)

$$\begin{eqnarray}&&{{ \mathcal E }}_{L}={{Ce}}^{2}{n}_{e}^{4/3}{Z}^{2/3},\end{eqnarray} \tag{ 3 }$$

where the lattice constant is C = −1.444, for a body-centered-cubic (bcc) structure, e is the electron charge, and n_e is the electron number density.

The total energy density ${ \mathcal E }$ of the system consists of the energy density of the ions, the degenerate electron energy density, and the energy density of the ionic lattice, ${{ \mathcal E }}_{L}$, according to

$$\begin{eqnarray}&&{ \mathcal E }={n}_{x}M(Z,A){c}^{2}-{n}_{e}{m}_{e}{c}^{2}+{{ \mathcal E }}_{L},\end{eqnarray} \tag{ 4 }$$

where n_x is the number density of atomic nuclei of mass M(Z, A) and m_e denotes the electron mass.

3.1. Inverse β Decay

In this section we describe the treatment of inverse β-decay processes in white dwarf matter. Since these processes depend much more strongly on the magnetic field than the bulk properties of WDs, the magnetic field is taken into account in the calculation of the inverse β-decay processes. As known from the work of Gamow (1939; see also Shapiro & Teukolsky 1983), the matter inside white dwarfs can experience instabilities due to inverse β-decay processes,

$$\begin{eqnarray*}&&A(N,Z)+{e}^{-}\to A(N+1,Z-1)+{\nu }_{e}.\end{eqnarray*}$$

Because of this reaction, atomic nuclei become more neutron-rich and as a consequence the electron energy density and pressure are reduced, leading to a softer EoS. Using the thermodynamic relation (at zero temperature) ${{ \mathcal E }}_{e}+{P}_{e}={n}_{e}{\mu }_{e}$, the Gibbs free energy per nucleon, g, is obtained as

$$\begin{eqnarray}&&g(A,Z)={{mc}}^{2}+\displaystyle \frac{M(Z,A){c}^{2}}{A}+{\gamma }_{e}\left[{\mu }_{e}-{m}_{e}{c}^{2}+\displaystyle \frac{4}{3}\displaystyle \frac{{{ \mathcal E }}_{L}}{{n}_{e}}\right],\end{eqnarray} \tag{ 5 }$$

where γ_e = Z/A is the ration between the proton number Z and atomic number A, with m being the neutron mass and μ_e being the electron chemical potential.

Inverse β-decay reactions are believed to occur in the cores of white dwarfs if the condition (Chamel et al. 2014)

$$\begin{eqnarray}&&g(A,Z)\geqslant g(A,Z-1)\end{eqnarray} \tag{ 6 }$$

is fulfilled, where g(A, Z) and g(A, Z − 1) follow from Equation (5).

From the inequality (6), we obtain the following relation (Chamel & Fantina 2015):

$$\begin{eqnarray}&&{\mu }_{e}+{{Ce}}^{2}{n}_{e}^{1/3}f(Z,Z-1)\geqslant {\mu }_{e}^{\beta },\end{eqnarray} \tag{ 7 }$$

with the electron number density n_e and mass density ρ of a magnetized electron gas given, respectively, by

$$\begin{eqnarray}&&{n}_{e}=\displaystyle \frac{2{B}_{\star }}{{\left(2\pi \right)}^{2}{\lambda }_{e}^{3}}\displaystyle \sum _{\nu }^{{\nu }_{\max }}{g}_{\nu 0}\sqrt{{x}_{F}^{2}-1-2\nu {B}_{\star }}.\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\rho =\displaystyle \frac{1}{{\gamma }_{e}}{{mn}}_{e},\end{eqnarray} \tag{ 9 }$$

where λ_e = ℏ/m_ec is the electron Compton wavelength. The sum over ν is over all occupied energy levels. Moreover, B_⋆ = B/B_c with B_c = 4.414 × 10¹³ G for the critical magnetic field. The quantities x_F in Equation (8) and ${\mu }_{e}^{\beta }$ in Equation (7) are defined as ${x}_{F}\equiv {p}_{F}/{m}_{e}c$, where p_F represents the Fermi momentum, and

$$\begin{eqnarray}&&{\mu }_{e}^{\beta }(A,Z)\equiv M(Z-1,A){c}^{2}-M(Z,A){c}^{2}+{m}_{e}{c}^{2}.\end{eqnarray} \tag{ 10 }$$

For an electron gas consisting of only one type of ion, we have

$$\begin{eqnarray}&&f(Z,Z-1)={Z}^{5/3}-{\left(Z-1\right)}^{5/3}+\displaystyle \frac{1}{3}{Z}^{2/3}.\end{eqnarray} \tag{ 11 }$$

In the limit where only the ground state (ν = 0) is fully occupied by electrons, one has

$$\begin{eqnarray*}&&{n}_{e}={n}_{{eB}}\propto {B}_{\star }^{2/3},\end{eqnarray*}$$

where n_eB is the number density when the ground state is fully occupied (see Haensel et al. (2007) for more details about n_eB). The chemical potential of the electrons in this case is given by

$$\begin{eqnarray}&&{\mu }_{e}\approx \displaystyle \frac{2{\pi }^{2}{m}_{e}{c}^{2}{\lambda }_{e}^{3}{n}_{{eB}}}{{B}_{\star }}.\end{eqnarray} \tag{ 12 }$$

Chamel et al. (2014) estimated that the maximum magnetic field inside a white dwarf, before the onset of β-inverse reactions, is given by

$$\begin{eqnarray}&&{B}_{\star }^{\beta }\approx \displaystyle \frac{1}{2}{\left(\displaystyle \frac{{\mu }_{e}^{\beta }(A,Z)}{{m}_{e}{c}^{2}}\right)}^{2}{\left[1+{\left(\displaystyle \frac{4}{\pi }\right)}^{2/3}\displaystyle \frac{C\alpha }{3}f(Z,Z-1)\right]}^{-2},\end{eqnarray} \tag{ 13 }$$

where α = e²/(ℏc). We note that because of the second term on the right side of Equation (13), which originates from lattice contributions, the maximum value of ${B}_{\star }^{\beta }$ increases if lattice contributions are taken into account.

In Figure 1 we show the numerical solution of (7) for white dwarf matter made of only carbon ions immersed in a magnetized electron gas. The oscillatory density behavior is caused by the Landau level contributions to the number density, given by Equation (8). For high values of B with only the ground state occupied, the dependence of density on B becomes linear, as can be seen in Figure 1.

Zoom In Zoom Out Reset image size

3.2. Pycnonuclear Reaction

In this section, we will focus on nuclear fusion reactions (pycnonuclear fusion reactions) among heavy atomic nuclei, schematically expressed as ${}_{Z}^{A}{ \mathcal Y }{+}_{Z}^{A}{ \mathcal Y }{\to }_{2Z}^{2A}{ \mathcal K }$. An example of such a reaction is carbon on carbon, ${}^{12}{\rm{C}}+{}^{12}{\rm{C}}$. Pycnonuclear reactions have been found to occur over a significant range of stellar densities (see Gasques et al. (2005)), including the density range found in the interiors of white dwarfs (Chamel et al. 2013, 2014). The nuclear fusion rates at which pycnonuclear reactions proceed, however, are highly uncertain because of some poorly constrained parameters (see Gasques et al. (2005 and Yakovlev et al. 2006). We note that magnetic fields affect the zero-point energy of ions and hence the pycnonuclear reaction rates. This effect, however, is neglected in this paper. The reaction rates have been calculated for different models. According to Gasques et al. (2005), the pycnonuclear reaction rates are defined as

$$\begin{eqnarray}&&{R}_{\mathrm{pyc}}=\displaystyle \frac{{n}_{x}}{2}S({E}_{\mathrm{pk}})\displaystyle \frac{{\hslash }}{{{mZ}}^{2}{e}^{2}}\,{P}_{\mathrm{pyc}}\,{F}_{\mathrm{pyc}},\end{eqnarray} \tag{ 14 }$$

where S(E_pk) is the astrophysical S-factor used by Gasques et al. (2005) for the NL2 nuclear model parameterization. Following Gasques et al. (2005), an analytic equation for the S-factor is given by

$$\begin{eqnarray}&&\begin{array}{l}S({E}_{\mathrm{pk}})\\ =\,5.15\times {10}^{16}\exp \left[-0.428{E}_{\mathrm{pk}}-\displaystyle \frac{3{E}_{\mathrm{pk}}^{0.308}}{1+{e}^{0.613(8-{E}_{\mathrm{pk}})}}\right],\end{array}\end{eqnarray} \tag{ 15 }$$

where S(E_pk) is in units of MeV barn. The factors P_pyc and F_pyc in Equation (14) are given by

$$\begin{eqnarray}&&{P}_{\mathrm{pyc}}=\exp \left(-{C}_{\exp }/\sqrt{{\rm{\Lambda }}}\right),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{pyc}}=8\,{C}_{\mathrm{pyc}}\,11.515\,/\,{{\rm{\Lambda }}}^{{C}_{\mathrm{pl}}},\end{eqnarray} \tag{ 17 }$$

with C_exp, C_pyc, and C_pl being dimensionless parameters for a regular bcc-type crystal lattice (at zero temperature). Their numerical values are listed in Table 1.

Table 1.  Coefficients C_exp, C_pyc, and C_pl Related to Pycnonuclear Reaction Rates at Zero Temperature, Computed for Nuclear Model NL2 (see Cândido Ribeiro et al. 1997; Chamon et al. 1997)

Model	Cexp	Cpyc	Cpl
bcc; static lattice	2.638	3.90	1.25

Download table as:  ASCIITypeset image

The inverse-length parameter Λ in Equations (16) and (17) has the form (Gasques et al. 2005; Yakovlev et al. 2006)

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Lambda }} & = & \displaystyle \frac{{{\hslash }}^{2}}{{{mZ}}^{2}{e}^{2}}{\left(\displaystyle \frac{{n}_{x}}{2}\right)}^{1/3}\\ & = & \displaystyle \frac{1}{{{AZ}}^{2}}{\left(\displaystyle \frac{1}{A}\displaystyle \frac{\rho {X}_{i}}{1.3574\times {10}^{11}{\rm{g}}{\mathrm{cm}}^{-1}}\right)}^{1/3}.\end{array}\end{eqnarray} \tag{ 18 }$$

For mass densities, ρ, less than the neutron drip density one has X_i = 1 (Gasques et al. 2005) and the pycnonuclear reaction rates are given by

$$\begin{eqnarray}\begin{array}{rcl}{R}_{\mathrm{pyc}} & = & \rho {{AZ}}^{4}S({E}_{\mathrm{pk}}){C}_{\mathrm{pyc}}{10}^{46}{{\rm{\Lambda }}}^{3-{C}_{\mathrm{pl}}}\\ & & \times \,\exp \left(-{C}_{\exp }/\sqrt{{\rm{\Lambda }}}\right),\end{array}\end{eqnarray} \tag{ 19 }$$

with R_pyc given in units of cm⁻³ s⁻¹. The zero-point oscillation energy, E_pk, of ¹²C nuclei at ρ = 10¹⁰ g cm⁻³ is given by Shapiro & Teukolsky (1983)

$$\begin{eqnarray}&&{E}_{\mathrm{pk}}={\hslash }\omega ={\hslash }{\left(\displaystyle \frac{4\pi {e}^{2}{Z}^{2}\rho }{{{AM}}^{2}}\right)}^{1/2}.\end{eqnarray} \tag{ 20 }$$

The time it takes for the complete fusion of atomic nuclei of mass Am is obtained from (Gasques et al. 2005; Boshkayev et al. 2013)

$$\begin{eqnarray}&&{\tau }_{\mathrm{pyc}}=\displaystyle \frac{{n}_{x}}{{R}_{\mathrm{pyc}}}=\displaystyle \frac{\rho }{\mathrm{Am}{R}_{\mathrm{pyc}}}.\end{eqnarray} \tag{ 21 }$$

As already mentioned above, the reaction rates are rather uncertain, and the analytic astrophysical S-factor has an uncertainty of ∼3.5, which considerably affects the density thresholds of pycnonuclear reactions and their reaction times. In Figure 2 we show pycnonuclear fusion reaction rates and pycnonuclear reaction timescales for carbon burning at zero temperature as a function of mass density. The bcc crystal lattice for nuclear model NL2 was employed to produce Figure 2.

Zoom In Zoom Out Reset image size

The numerical technique used in this work to study axisymmetric magnetic fields was first applied to neutron stars (Bonazzola et al. 1993; Bocquet et al. 1995), and more recently (Franzon & Schramm 2015b; Franzon et al. 2016a, 2016b; Chatterjee et al. 2017). The same formalism was used to study rotating and magnetized white dwarfs by Franzon & Schramm (2015a). Here we compute stellar equilibrium configurations by solving the Einstein–Maxwell field equations in a fully general relativistic treatment. It is worth noting that, although the maximum mass is not affected by general relativistic effects, the stellar radius can be very different; see, for instance, Carvalho et al. (2018b). For more details about the theoretical formalism and numerical procedure, see, for instance, Gourgoulhon (2012). Below we show the basic electromagnetic equations, combined with the gravitational equations, are solved numerically by means of a spectral method (LORENE⁹ ). The stress-energy tensor T_αβ is composed of the matter and the electromagnetic source term,

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\alpha \beta } & = & ({ \mathcal E }+P){u}_{\alpha }{u}_{\beta }+{{Pg}}_{\alpha \beta }\\ & & +\,\displaystyle \frac{1}{{\mu }_{0}}\left({F}_{\alpha \mu }{F}_{\beta }^{\mu }-\displaystyle \frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }{{\rm{g}}}_{\alpha \beta }\right).\end{array}\end{eqnarray} \tag{ 22 }$$

Here F_αμ is the antisymmetric Maxwell tensor defined as ${F}_{\alpha \mu }={\partial }_{\alpha }{A}_{\mu }-{\partial }_{\mu }{A}_{\alpha }$, with A_μ denoting the electromagnetic four-potential A_μ = (A_t, 0, 0, A_ϕ). The total energy density of the system is ${ \mathcal E }$, the pressure is denoted by P, u_α is the fluid 4-velocity, and the metric tensor is g_αβ. The first term in Equation (22) represents the isotropic (ideal) matter contribution to the energy momentum-tensor, while the second term is the anisotropic electromagnetic field contribution.

The metric tensor in axisymmetric spherical-like coordinates $(r,\theta ,\phi )$ can be read off from the line element

$$\begin{eqnarray}\begin{array}{rcl}{{ds}}^{2} & = & -{N}^{2}{{dt}}^{2}+{{\rm{\Psi }}}^{2}{r}^{2}{\sin }^{2}\theta {\left(d\phi -{N}^{\phi }{dt}\right)}^{2}\\ & & +\,{\lambda }^{2}({{dr}}^{2}+{r}^{2}d{\theta }^{2}),\end{array}\end{eqnarray} \tag{ 23 }$$

where N, N^ϕ, Ψ, and λ are functions of the coordinates (r, θ) (Bonazzola et al. 1993). The gravitational field is derived from the integration of a coupled system of four elliptic partial differential equations for the metric functions. The final system of gravitational equations can be put in the form

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{2}[(N{\rm{\Psi }}-1)r\sin \theta ]=8\pi N{\lambda }^{2}{\rm{\Psi }}r\sin \theta ,({S}_{r}^{r}-{S}_{\theta }^{\theta }),\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{2}[\mathrm{ln}\lambda +\nu ]=8\pi {\lambda }^{2}{S}_{\phi }^{\phi }+\displaystyle \frac{3{{\rm{\Psi }}}^{2}{r}^{2}{\sin }^{2}\theta }{4{N}^{2}}\partial \omega \partial \omega -\partial \nu \partial \nu ,\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{3}\nu =4\pi {\lambda }^{2}(E+S)+\displaystyle \frac{{{\rm{\Psi }}}^{2}{r}^{2}{\sin }^{2}\theta }{2{N}^{2}}\partial \omega -\partial \nu \partial \nu (\nu +\mathrm{ln}{\rm{\Psi }}),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\begin{array}{l}\left[{{\rm{\Delta }}}_{3}-\displaystyle \frac{1}{{r}^{2}{\sin }^{2}\theta }\right]=-16\pi \displaystyle \frac{N{\lambda }^{2}}{{{\rm{\Psi }}}^{2}}\displaystyle \frac{{J}_{\phi }}{r\sin \theta }r\\ \quad \times \,\sin \theta \partial \omega \partial (\nu -3\mathrm{ln}{\rm{\Psi }}),\end{array}\end{eqnarray} \tag{ 27 }$$

where we have introduced the definitions

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{2}=\displaystyle \frac{{\partial }^{2}}{\partial {r}^{2}}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {\theta }^{2}},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{3}=\displaystyle \frac{{\partial }^{2}}{\partial {r}^{2}}+\displaystyle \frac{2}{r}\displaystyle \frac{\partial }{\partial r}+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {\theta }^{2}}+\displaystyle \frac{1}{{r}^{2}\tan \theta }\displaystyle \frac{\partial }{\partial \theta },\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\nu =\mathrm{ln}N.\end{eqnarray} \tag{ 30 }$$

In addition, for the gravitational field equations (Equations (24)–(27)), ∂ω∂ω is shorthand for

$$\begin{eqnarray}&&\partial \omega \partial \omega =\displaystyle \frac{\partial \omega }{\partial r}\displaystyle \frac{\partial \omega }{\partial r}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial \omega }{\partial \theta }\displaystyle \frac{\partial \omega }{\partial \theta },\end{eqnarray} \tag{ 31 }$$

and the total energy density, momentum density, and stress tensors of the system are

$$\begin{eqnarray}&&E={{\rm{\Gamma }}}^{2}\left({ \mathcal E }+P\right)-P+{E}^{{EM}},\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{J}_{\phi }=(E+P){\lambda }^{2}{\rm{\Psi }}r\sin \theta U+{J}_{\phi }^{{EM}},\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{S}_{r}^{r}=P+{S}_{r}^{{EM}\,r},\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{S}_{\theta }^{\theta }=P+{S}_{\theta }^{{EM}\,\theta },\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{S}_{\phi }^{\phi }=P+(E+P){U}^{2}+{S}_{\phi }^{{EM}\,\phi },\end{eqnarray} \tag{ 36 }$$

where U represents the fluid velocity, and Γ is the Lorentz factor, which connects the Eulerian and fluid co-moving observers. Finally ${E}^{{EM}},\,{J}_{\phi }^{{EM}},\,{S}_{r}^{{EM}\,r},\,{S}_{\theta }^{{EM}\,\theta }$, and ${S}_{\phi }^{{EM}\,\phi }$ are the electromagnetic contributions to the energy, momentum, and the stress tensor of the system (for further details, see, e.g., Franzon & Schramm 2015b; Franzon et al. 2016a, 2016b).

As in Bonazzola et al. (1993), the equation of motion for a star endowed with a magnetic field, reads

$$\begin{eqnarray}&&H\left(r,\theta \right)+\nu \left(r,\theta \right)+M\left(r,\theta \right)=\mathrm{const},\end{eqnarray} \tag{ 37 }$$

where H(r, θ) is the heat function defined in terms of the baryon number density n,

$$\begin{eqnarray}&&H={\int }_{0}^{n}\displaystyle \frac{1}{e({n}_{1})+p({n}_{1})}\displaystyle \frac{{dP}}{{dn}}({n}_{1}){{dn}}_{1}.\end{eqnarray} \tag{ 38 }$$

The quantity ν(r, θ) in Equation (37) is defined as in Equation (26), and the magnetic potential M(r, θ) is given by

$$\begin{eqnarray}&&M\left(r,\theta \right)=M\left({A}_{\phi }\left(r,\theta \right)\right)\equiv -{\int }_{{A}_{\phi }\left(r,\theta \right)}^{0}f\left(x\right){dx},\end{eqnarray} \tag{ 39 }$$

where f(x) denotes the current function. Magnetic stellar models are obtained by assuming a constant value, f₀, for the latter (Franzon & Schramm 2015b). According to Bocquet et al. (1995), other choices for f(x) are possible, but the general conclusions presented in this work remain the same. The constant current function is a standard way to self-consistently generate a dipolar magnetic field throughout the star.

In this section, we discuss the effects of strong magnetic fields on the global properties of stationary white dwarfs, taking into account instabilities due to inverse β-decay and pycnonuclear fusion reactions in their cores. The instabilities related to the microphysics are fundamental, as they put strong constraints on the equilibrium configurations and also limit the maximum magnetic fields that these stars can have (Chamel et al. 2013). In addition to the magnetic profiles, which have already been computed by Franzon & Schramm (2015a), we also compute stellar models at constant magnetic dipole moments μ. Franzon & Schramm (2015a) used a simple Fermi gas to model the EoS of WDs, but the microphysical issues were not addressed. In our study, the maximum white dwarf mass for non-magnetized stars is smaller than that considered in Franzon & Schramm (2015a), as the lattice contribution softens the EoS. The density threshold at which the inverse beta-decay sets in is taken to be ρ ∼ 3 × 10¹⁰ g cm⁻³, which is the minimum value at which the inverse beta-decay process occurs (see Figure 1).

In Figure 3, we show the gravitational mass versus central density of white dwarf sequences computed for different fixed magnetic dipole moments, μ, and current functions, f₀. The magnetic dipole moment is defined as (see Bonazzola et al. (1993))

$$\begin{eqnarray}&&\displaystyle \frac{2\mu \cos \theta }{{r}^{3}}=B(r){| }_{r\to \infty },\end{eqnarray} \tag{ 40 }$$

which is the radial (orthonormal) component of the magnetic field of a magnetic dipole seen by an observer at infinity.

Zoom In Zoom Out Reset image size

As can be seen in Figure 3, the masses of magnetized white dwarfs, with fixed magnetic dipole moments, increase monotonically with central density. This behavior is very different if the value of the current function is kept constant, in which case non-monotonic mass-density relationships are obtained. The cross-hatched area in Figure 3 shows the density regime where pycnonuclear fusion reactions become possible. The threshold density (9.25 × 10⁹ g cm⁻³) at which this reaction occurs is marked with a solid black square in Figure 3. The pycnonuclear reaction timescale considered to obtain the threshold density was 10 Gyr. For a central white dwarf density of 1.59 × 10¹⁰ g cm⁻³ the fusion reaction timescale decreases to 0.1 Myr (see Figure 2), which we considered an unstable configuration, because at this density the pycnonuclear reactions may induce a supernova outburst. White dwarfs subject to inverse β-decay reactions in their cores are located in the yellow area (marked "β-inverse") of Figure 3. An overview of the density thresholds discussed above is provided in Table 2 for white dwarfs with different magnetic field values and magnetic dipole moments. To derive the maximum stable masses the stability condition for non-magnetized white dwarfs, dM/dρ_c > 0, is employed, because it is also valid for sequences with constant magnetic dipole moment (Sorkin 1982). The curves in Figure 3 with fixed μ respect the stability condition, but the last point (point of maximum mass) on each curve has dM/dρ_c < 0. The most massive stable white dwarf that is not subject to microscopic instability reactions and dynamical instability (end point of the curve with μ = 2 × 10³⁴ Am²) has a mass of ∼2.14 M_⊙ and an equatorial radius of ∼1096 km; see also Figure 4).

Zoom In Zoom Out Reset image size

The mass–radius relationship of magnetized white dwarfs, for different (fixed) magnetic dipole moments μ, is shown in Figure 4. One can see that increasing values of μ lead to white dwarfs with larger radii, because of the added magnetic field energy. The strength of the magnetic field can be inferred from Figure 5, which shows the gravitational mass as a function of surface (B_s) and central (B_c) magnetic fields, the circumferential equatorial radius (R_circ), and the baryon number density (n_b), for two sample of magnetic dipole moments of μ = 0.5 × 10³⁴ Am² (red line) and μ = 4.0 × 10³⁴ Am² (orange line).

Zoom In Zoom Out Reset image size

In Figure 5 (top panels), the curves with μ = 0.5 × 10³⁴ Am² and μ = 4.0 × 10³⁴ Am² cross each other. This is due to the fact that the magnetic field scales as ∼μ/r³, with r being the stellar radius (see Equation (40)). The locations of stars with fixed baryon masses of M_B = 1.00 M_⊙ and M_B = 1.80 M_⊙ are shown in Figure 5 by dashed horizontal lines. According to Equation (40), the magnetic field is determined by the size of the star along the curves with constant μ. However, along the lines with fixed baryon masses the strength of the magnetic field is a combination of the magnetic dipole moment μ and the stellar radius r.

Next, we discuss the behavior of the magnetic dipole moments of white dwarfs whose magnetic fields are weakening. From the M versus B_s and M versus B_c relationships shown in Figure 5 (top panels), two different scenarios are possible, depending on the mass and the magnetic field strengths of white dwarfs. If located above the crossing point of the $\mu =0.5\times {10}^{34}\,{\mathrm{Am}}^{2}$ (red line) and μ = 4.0 × 10³⁴ Am² (orange line) curves, white dwarfs with weakening magnetic fields evolve from right to left in the two upper panels of Figure 5, as shown (back arrow) for a white dwarf with a constant baryon mass of M_B = 1.80 M_⊙. The magnetic dipole moment of such white dwarfs increases from μ = 0.5 × 10³⁴ Am² to μ = 4.0 × 10³⁴ Am². This is accompanied by an increase of the stellar radius (see M versus R curve) and a decrease of the central baryon density (see M versus n_b curve). The situation is reversed for white dwarfs located below the crossing. For those white dwarfs, a reduction of the magnetic field is accompanied by a decrease of the magnetic dipole moment, as shown in Figure 5 for a sample white dwarf with a constant baryon mass of M_B = 1.00 M_⊙ (black arrows). In this case, white dwarfs become smaller and therefore more dense at the center (see M versus R and M versus n_b curves shown in Figure 5).

As discussed, in (Figure 5), the equatorial radii of white dwarfs located above the crossing point increase as their magnetic fields become smaller. This increase in radius (at a fixed baryon mass) is due to the Lorentz force. However, the stellar magnetic field scales as μ/r³, which means that for a star with a mass of M_B = 1.80 M_⊙ the increase in the magnetic dipole moment is canceled by the increase in the radius, thus reducing the magnetic field. This is the opposite of what is expected for stars with lower masses. For example, a star with M_B = 1.00 M_⊙ decreases its magnetic dipole moment and its radius. However, in this case, the decrease in radius is not enough to cancel the reduction in μ. The net result is a decrease of the magnetic field. This can be understood by looking at the variation in the circular equatorial radius of the stars with M_B = 1.80 M_⊙ and M_B = 1.00 M_⊙. For the latter, the change in radius is much smaller than the radial change of a star with M_B = 1.80 M_⊙ and a change in the magnetic dipole moment of $| {\rm{\Delta }}\mu | =3.5\times {10}^{34}$ Am².

In Figure 6, we compare the global properties of white dwarfs with fixed baryon masses of M_B = 1.00 M_⊙ and M_B = 1.80 M_⊙. The top panels show the central baryon density as a function of the central magnetic field (top left panel) and the circular equatorial radius (top right panel) for a white dwarf with M_B = 1.80 M_⊙. For such stars, as the magnetic field decreases, the central baryon density becomes smaller due to the fact that the radius is increasing. On the other hand, for lighter white dwarfs, with a mass of M_B = 1.00 M_⊙, the central baryon number density increases as the magnetic field decreases, because the stellar radius becomes smaller.

Zoom In Zoom Out Reset image size

In this work, we presented axisymmetric and stationary models of magnetized white dwarfs obtained by solving the Einstein–Maxwell equations self-consistently and taking into account stability considerations related to neutronization, due to electron capture reactions, as well as pycnonuclear fusion reactions among carbon nuclei in the cores of white dwarfs.

We also investigated the influence of magnetic fields on the structure of white dwarfs. This is an important problem, because super-massive magnetized white dwarfs, whose existence is partially supported by magnetic forces, could simplify the explanation of observed ultra-luminous explosions of type Ia supernovae. The Lorentz force induced by strong magnetic fields breaks the spherical symmetry of stars and increases their masses, as the force acts in the radial outward direction against the inwardly directed gravitational pull.

In this paper, we employed the EoS for a degenerate electron gas with electron-ion interactions (body-centered-cubic lattice structure) to describe the matter inside white dwarfs. We have shown that the EoS becomes softer if nuclear lattice contributions are included in addition to the electron pressure. This is due to the fact that the repulsive force between electrons is smaller in the presence of an ionic lattice, causing a softening of the EoS. We note that the density thresholds for pycnonuclear fusion reactions and inverse β-reactions are reduced when magnetic fields are present in the stellar interior, as can be seen in Table 2.

We have shown that the masses of white dwarfs can increase up to M = 2.14 M_⊙ (with a corresponding magnetic dipole moment of μ = 2.0 × 10³⁴ Am² (see, e.g., Figure 3) even when microphysical instabilities are considered. In particular, this white dwarf star has an equatorial radius of ∼1100 km with magnetic fields of B_c = 3.85 × 10¹⁴ G and B_s = 7.21 × 10¹³ G at the center and at the stellar surface, respectively. For this white dwarf, the ratio between the magnetic pressure and the matter pressure at the center is 0.789. Although the surface magnetic fields obtained here are higher than the observed ones for white dwarfs, these figures provide an idea of the maximum possible magnetic field strength that can be reached inside of these objects and may also be used to assess the effects of strong magnetic fields on both the microphysics and the global structure of magnetized white stars.

The maximum magnetic field found in this work is an order of magnitude smaller than that established by Franzon & Schramm (2015a). This is because we modeled the stellar interior with a more realistic EoS than just a simple electron gas. In addition, the magnetic field here is slightly larger than that found in Chatterjee et al. 2017, as the pycnonuclear reaction rates have been calculated using a different model. The pycnonuclear reaction rates are defined using the astrophysical S-factor for the NL2 nuclear model for the carbon element, while in Chatterjee et al. 2017 they considered the mass-density thresholds for the onset of electron capture by the daughter nuclei ${}_{2Z}^{2A}{ \mathcal K }$, a difference that allowed us to reach values of mass densities before instability larger than the values used in Chatterjee et al. 2017, hence also yielding larger masses. In addition, we considered the density threshold for pycnonuclear fusion reactions for a 10 Gyr fusion reaction timescale, which restricts the central density of white dwarfs to ∼9.25 × 10⁹ g cm⁻³ (see Table 2), thus limiting the stellar masses, and therefore the radii, which, for very massive and magnetized white dwarfs, cannot be smaller than R ∼ 1100 km. However, it is important to mention that the pycnonuclear reaction timescales are somewhat uncertain. In our case we have a factor of uncertainty of approximately 3.5 in the calculation of the astrophysical S-factor (see (Gasques et al. 2005; Yakovlev et al. 2006)).

Our results show that the surface magnetic field, B_s, is about one order of magnitude smaller than the magnetic field reached at the stellar center, B_c. If the magnetic field weakens for massive white dwarfs, we found that the magnetic dipole moments of such stars may increase (Figure 5), which is due to the fact that for a fixed baryon mass, the magnetic field is determined by the interplay between the magnetic dipole moment and the stellar radius. The situation is reversed for less massive white dwarfs, for which smaller magnetic fields imply smaller stellar magnetic dipole moments. The radii of massive (light) white dwarfs are found to increase (decrease) for decreasing central magnetic fields (Figure 6). This opens up the possibility that massive white dwarfs, with central magnetic fields greater than B ∼ 10¹³ G, increase their magnetic fields through continued compression. This phenomenology differs from previous studies carried out for magnetic fields less than ∼10¹³ G (Ostriker & Hartwick 1968; Suh & Mathews 2000), where an increase of the central magnetic field was found to make stars less dense and therefore bigger in size.

We note that stellar configurations that contain only poloidal magnetic fields (no toroidal component) are unstable (see, e.g., Braithwaite 2006; Armaza et al. 2015; Mitchell et al. 2015). Moreover, according to Goldreich & Reisenegger (1992), many different mechanisms can affect the magnetic fields and their distributions inside of white dwarfs. In this work, in the framework of a fully general relativistic treatment, we model the properties of magnetized white dwarfs with purely poloidal magnetic field components. Although this is not the most general magnetic field profile and a dynamical stability of these stars still needs to be addressed, magnetic fields considerably increase the masses of white dwarfs, even when microphysical instabilities are taken into account. As a consequence, such white dwarfs should be considered possible candidates of super-Chandrasekhar white dwarfs, thereby contributing to our understanding of superluminous type Ia supernovae.

Lastly, we note that for a typical magnetic field value of ∼10¹⁴ G and a density of ∼10⁹ g cm⁻³, one obtains an Alfvén velocity of v = 10⁹ cm s⁻¹, which, for a white dwarf with a typical radius of R = 1500 km, leads to an Alfvén crossing time of ∼0.1 s (Durisen 1973; Yakovlev & Urpin 1980; Cumming 2002). This is close to the hydrostatic equilibration time of white dwarfs. As a consequence, although magnetized white dwarfs seem to be short-lived stars, they might still be supported by magnetic fields. Our results represent magnetostatic equilibrium conditions. The stability analysis of such systems is beyond the scope of this study, which aims for a more complete discussion of the possible existence of super-Chandrasekhar white dwarfs. Studies that address those issues, such as the role of different (poloidal and toroidal) magnetic field configurations, stellar rotation, and different compositions of the stellar cores, will be presented in a series of forthcoming papers.

We express our deep appreciation to the late Prof. S. Schramm, whose contribution to this work was of great significance.

We acknowledge financial support from the Brazilian agencies CAPES, CNPq, and thank FAPESP for financial support under the thematic project 13/26258-4 B. Franzon acknowledges support from CNPq/Brazil, DAAD and HGS-HIRe for FAIR. G.A. Carvalho thanks CAPES for financial support process CAPES/PDSE # 88881.188302/2018-01. S. Schramm acknowledges support from the HIC for FAIR LOEWE program. F. Weber is supported by the National Science Foundation (USA) under grant PHY-1714068.