Evolution of Highly Magnetic White Dwarfs by Field Decay and Cooling: Theory and Simulations

We investigate the luminosity suppression and its effect on the mass–radius relation and cooling evolution of highly magnetized white dwarfs. Based on the effect of magnetic field relative to gravitational energy, we suitably modify our treatment of the radiative opacity, magnetostatic equilibrium, and degenerate core equation of state to obtain the structural properties of these stars. Although the Chandrasekhar mass limit is retained in the absence of magnetic field and irrespective of the luminosity, strong central fields of about 1014 G can yield super-Chandrasekhar white dwarfs with masses ∼2.0 M ⊙. Smaller white dwarfs tend to remain super-Chandrasekhar for sufficiently strong central fields even when their luminosity is significantly suppressed to 10−16 L ⊙. Nevertheless, owing to the cooling evolution and simultaneous field decay over 10 Gyr, the limiting masses of small magnetized white dwarfs can fall to 1.5 M ⊙ over time. However, the majority of these systems still remain practically hidden throughout their cooling evolution because of their high fields and correspondingly low luminosities. Utilizing the stellar evolution code stars, we obtain close agreement with the analytical mass limit estimates, which suggests that our analytical formalism is physically motivated. Our results argue that super-Chandrasekhar white dwarfs born as a result of strong-field effects may not remain so forever. This explains their apparent scarcity, in addition to making them hard to detect because of their suppressed luminosities.


Introduction
Observations of more than a dozen overluminous Type Ia supernovae (SNe Ia; Howell et al. 2006;Scalzo et al. 2010) strongly indicate the existence of their very massive progenitors with masses M > 2 M e . SNe Ia are widely studied because of their utility as standard candles to estimate cosmological distances. While binary evolution of accreting/rapidly differentially rotating white dwarfs (WDs ;Hachisu 1986;Yoon & Langer 2004) has already been used to explain such high-mass progenitors, none of these models could account for progenitor masses as high as 2.8 M e , which are inferred from the observations. Therefore, an alternate but equally exciting proposal relates to the highly magnetized super -Chandrasekhar WDs (hereafter B-WDs). In addition to being the super-Chandrasekhar mass progenitors of overluminous SNe Ia, B-WDs have also been proposed as promising candidates for soft gamma-ray repeaters and anomalous X-ray pulsars that have significantly lower magnetic fields and ultraviolet luminosities (Mukhopadhyay & Rao 2016).
Previous studies (Das & Mukhopadhyay 2012Subramanian & Mukhopadhyay 2015) already showed that strong magnetic fields can appropriately modify the equation of state (EOS) of the electron-degenerate matter to yield super-Chandrasekhar WDs with M ≈ 2.6 M e , with or without sufficiently rapid rotation. Interestingly, a misalignment between the rotation and magnetic axes can generate, along with dipole radiation, significant gravitational radiation, which can be detected by space-based gravitational wave detectors, leading to a direct detection of super-Chandrasekhar WDs (Kalita & Mukhopadhyay 2019). Observations confirm that magnetized WDs are indeed more massive than other nonmagnetized WDs (Ferrario et al. 2015). Furthermore, data from the Sloan Digital Sky Survey (SDSS) suggest that, in addition to having marginally higher masses, magnetized WDs span a similar effective temperature range to their nonmagnetic counterparts (Vanlandingham et al. 2005).
The effect of strong magnetic fields on the WD mass-radius relation has already been explored in some detail previously by our group, for both Newtonian (Das & Mukhopadhyay 2012) and general relativistic formalisms (Das & Mukhopadhyay 2014Subramanian & Mukhopadhyay 2015). These investigations were carried out for various different magnetic field configurations and were in good agreement with the results from independent studies that indicate the existence of super-Chandrasekhar WDs (Boshkayev et al. 2013;Franzon & Schramm 2015;Carvalho et al. 2018). Otoniel et al. (2019) recently studied potential matter instabilities for these B-WDs in a general relativistic framework, particularly in relation to the pycnonuclear and electron capture reactions. While they found that the limiting mass for nonrotating B-WDs is about 2.14 M e , and potentially larger if rotation is considered, the pycnonuclear reactions are likely to destabilize the star once the central density exceeds about 10 10 g cm −3 .
It should be noted that magnetized WDs have many important implications apart from their apparent link to overluminous SNe Ia, and hence their other properties should also be explored (Mukhopadhyay & Rao 2016;Mukhopadhyay et al. 2017aMukhopadhyay et al. , 2017b. Apart from increasing the limiting mass of WDs, strong magnetic fields can also influence the thermal properties, such as luminosity, temperature gradient, and cooling rate of the star. Bhattacharya et al. (2018a) explored the luminosities and cooling rates of B-WDs with the theoretical model proposed by Shapiro & Teukolsky (1983) for nonmagnetic WDs. Assuming that the interface properties are similar for both B-WDs and nonmagnetized WDs for a given stellar age and a nonzero temperature gradient across the surface layers, they showed that the luminosity for B-WDs can be suppressed significantly up to about 10 −9 L e for large fields B  10 12 G. They also computed the cooling rates for these B-WDs with suppressed luminosities and showed that their cooling evolution is significantly slowed down at large B. Indeed, Valyavin et al. (2014) analyzed the optical data for cool WD 1953−011 and found that strong fields suppress convection over the entire B-WD surface and thereby attenuate the cooling rate.
The initial exploration into the cooling of nonmagnetized WDs started with attempts to model the degenerate core as the primary source of thermal energy, which is then radiated away as the observed luminosity from the surface layers as the star gradually evolves over time (Mestel 1952;Mestel & Ruderman 1967). Tutukov & Yungelson (1996) calculated the cooling curves for low-mass WDs starting from a luminous star and evolving to the crystallization stage after about 10 Gyr, while Fontaine et al. (2001) discussed some limitations of the Mestel (1952) model in the context of WD cosmochronology. However, it is important to note that these studies either did not consider the effects of strong fields or assumed that the underlying fields were just too weak to have any practical effects on the WD cooling evolution.
Recently, Gupta et al. (2020) revisited the physics of luminosity suppression and the mass-radius relation in the context of B-WDs. In contrast to Bhattacharya et al. (2018a), they relaxed the assumption of fixing interface parameters and a preassigned mass or radius across all B-WDs and nonmagnetized WDs. They further extended their analysis to compute the radial profiles of the B-WD thermal properties for both the nondegenerate surface layers and the electrondegenerate B-WD cores. However, they ignored the effect of strong fields on the EOS of the degenerate core electrons, as well as the correction to the total B-WD mass by general relativistic effects.
Here, in a considerably more generalized framework, we model the B-WD structure properties from the center to the surface by solving the magnetostatic equilibrium, mass conservation, and photon diffusion equations simultaneously. We investigate the effect of the temperature gradient (directly related to the luminosity) on the mass-radius relation by considering both radiative and convective cooling. While the total pressure is the sum of the contributions from the electrondegenerate, ideal gas and magnetic pressures, the interface location is taken to be the radius where the degenerate pressure is roughly comparable to the ideal gas pressure.
In order to distinguish between the weakly and strongly magnetized cases, we compare the relative energy densities in the magnetic and gravitational fields. Accordingly, we modify the radiative opacity, magnetostatic balance equation, and EOS for the electron-degenerate core. This paper is organized as follows. In Section 2, we provide a brief overview of the method used to obtain the radial profiles of the WD properties and the mass-radius relation for a given luminosity. In Section 3, we discuss the physical effects of the magnetic field for both the weak-and strong-field limits, based on a comparison with the gravitational energy of the B-WD. Subsequently, in Section 4, we evaluate the suppressed luminosity of strongly magnetized B-WDs, after including the effects of cooling evolution and magnetic field decay by dissipative processes over long timescales. Next, we explore a set of numerical models produced using the stellar evolution code STARS in order to validate our analytical approach in Section 5. Finally, we conclude with a summary of the main results in Section 6.

White Dwarf Structure Properties
In this section, we describe the physical considerations used to self-consistently obtain the structure properties of magnetized WDs. We formulate a method to solve the magnetostatic equilibrium, mass conservation, and photon diffusion equations for a given luminosity and magnetic field configuration. We model the total pressure by including the contributions from the degenerate electron gas (dominant within the isothermal core), ideal gas (dominant in the surface layer), and magnetic pressures. The interface is defined to be at the WD radius, where the contributions from the inner electron-degenerate core and outer ideal gas pressures are equal. For our study, we consider radially varying magnetic fields that are realistic (see Deb et al. 2021). The presence of strong fields inside compact stars gives rise to additional pressure P B = B 2 /8π and density is the strength of the magnetic field (Sinha et al. 2013;Bhattacharya et al. 2018aBhattacharya et al. , 2018b. We assume the B-WDs to be approximately spherical. The assumption of spherical symmetry is justified if the central field does not exceed about 10 14 G (Subramanian & Mukhopadhyay 2015), which is toroidally dominated, so the model equations characterizing magnetostatic equilibrium, photon diffusion, and mass conservation can be written within a Newtonian framework as where the magnetic tension terms are ignored for a radially varying B. In these equations, P deg and P ig = ρkT/μm p are the electron degeneracy pressure established by Chandrasekhar (1935) and the ideal gas pressure, respectively, ρ is the matter density, k is Boltzmann's constant, T is the temperature, μ ≈ 2 is the mean molecular mass per electron, m p is the proton mass, G is Newton's gravitational constant, m(r) is the mass enclosed within radius r, a is the radiation constant, c is the speed of light, κ is the radiative opacity, L r is the luminosity at radius r, and γ is the adiabatic index of the gas. The second term on the right-hand side of Equation (2) is the convective cooling contribution, which can be dominant over the radiative photon cooling at some combinations of T and P = P deg + P ig + P B . Although strong central magnetic fields can potentially impede convection (Canuto & Mazzitelli 1992;Solanki 2003), a large fraction of the total energy flux can still be efficiently transported across the B-WD. Hence, we consider its effect on the stellar temperature profile in addition to the radiative cooling. The opacity in the surface layers of a nonmagnetized WD can be approximated by Kramers's formula, κ = κ 0 ρT −3.5 , where κ 0 = 4.34 ×10 24 Z(1 + X) cm 2 g −1 and X and Z are the mass fractions of hydrogen and heavy elements (other than hydrogen and helium) in the stellar interior, respectively. Assuming helium WDs for our purposes here, we set the helium mass fraction to Y = 0.9 and Z = 0.1 for simplicity. The radiative opacity in the surface layers is primarily due to the bound-free and free-free transitions of electrons (Shapiro & Teukolsky 1983). The radiation conduction typically dominates over the electron conduction in these regions, and hence the same goes for the corresponding opacities (Potekhin & Yakovlev 2001). However, in the presence of strong magnetic fields, the radiative opacity depends strongly on B, as we shall discuss in the next section.
It should be noted that there will be residual currents generated owing to the variation of magnetic field within the star. However, the entire degenerate B-WD core is essentially isothermal owing to its very large thermal conductivity, which also leads to the frozen magnetic flux. As a consequence, the heating ability is minimal given the near-superconducting nature of the degenerate core, and magnetic dissipation turns out to be negligible.
While a large number of B-WDs with surface fields up to 10 9 G or so have already been discovered by SDSS , it is likely that the central fields are several orders of magnitude larger (Fujisawa et al. 2012;Das & Mukhopadhyay 2014;Subramanian & Mukhopadhyay 2015). This is expected for residual field arising from the fossil field of the original star, which had a stronger core field in addition to dynamo effects that can replenish those fields (see Potter & Tout 2010;Mukhopadhyay et al. 2021). In order to capture the variation of field magnitude radially within the B-WD, here we adopt a profile used extensively to model magnetized neutron stars (NSs) and B-WDs (Bandyopadhyay et al. 1997;Das & Mukhopadhyay 2014;Deb et al. 2021), Here B s is the surface magnetic field, B 0 is a fiducial magnetic field, and η and γ are dimensionless parameters that determine how the magnetic field changes from the core to the surface. As ρ → 0 close to the WD surface, B → B s , whereas ρ → ρ c near the B-WD core that leads to B → B 0 . For our analysis here, we set ρ 0 = 10 9 g cm −3 , η = 0.8, and γ = 0.9 for all calculations, following Bhattacharya et al. (2018aBhattacharya et al. ( , 2018b. The profile in Equation (4) essentially indicates the magnitude of the field at various density points within the star and hence radial coordinates. Here we neglect effects such as offset dipoles and magnetic spots that can arise with more complex field configurations (Maxted & Marsh 1999;Vennes et al. 2003). Hereafter, we denote the magnetic field as B = (B s , B 0 ) in order to specify both the surface and core fields.
It should be noted that the model field profile given by Equation (4) is not a unique choice and alternate profiles have been explored in the literature, especially for magnetized NSs. In particular, Dexheimer et al. (2017) showed that the magnetic field profile can be well approximated by a quadratic polynomial in the baryon chemical potential instead of an exponential function of matter density. They presented a realistic distribution for a poloidal magnetic field in the polar direction, for which the magnetic field distribution selfconsistently satisfies the Einstein-Maxwell field equations.
There is no hydrogen burning or other nuclear fusion reactions taking place within the WD core, so we assume that the radial luminosity is constant, L r = L, where L is the surface luminosity. The degenerate electrons in the WD core generally have a large mean free path owing to the filled Fermi sea, and therefore their high thermal conductivity leads to a uniform temperature throughout the region (Shapiro & Teukolsky 1983). In the case of magnetized NSs, the thermal conduction can be suppressed along directions transverse to the magnetic field lines (Hernquist 1985;Potekhin et al. 2007). However, these changes in conduction rates are unlikely to affect the cooling process in B-WDs because the insulating region is nondegenerate and thermal conduction occurs only in the stellar interior (Tremblay et al. 2015). Furthermore, the average magnetic fields considered here are much weaker than those in magnetized NSs, and so we assume that the core is perfectly isothermal for B-WDs. Gupta et al. (2020) considered speculative cases with dT/dr ≠ 0 below the interface to show that even nonmagnetized WDs can have super-Chandrasekhar masses for sufficiently large luminosities. However, for the purpose of our work here, we ignore this possibility because there is no existing observational evidence to indicate that nonmagnetized and/or nonrotating WDs can have super-Chandrasekhar masses.

Nonmagnetic Results
First, we explore some basic features of nonmagnetic WDs considering B = (0, 0). We solve the set of differential Equations (1)-(3) for model WDs with the Runge-Kutta method, by providing the surface density, mass, and surface temperature as the boundary conditions. For a given WD surface luminosity L and its corresponding radius R, the surface temperature is obtained with the Stefan-Boltzmann law as T L R 4 s 2 1 4 ( ) p s = , where σ is the Stefan-Boltzmann constant. We set the surface density to ρ s = 10 −4 g cm −3 as representative for all the cases considered here. Following Gupta et al. (2020), the total mass is obtained by iteratively solving Equations (1)-(3) until the integrated mass starting from the WD surface and extending up to 10 km from the center matches the mass that is obtained by solving the mass conservation Equation (3) independently using the solution for the density profile.
The left panel of Figure 1 shows the effect of luminosity and thereby the temperature gradient (see Equation (2)) on the mass-radius relation for nonmagnetized WDs compared to Chandrasekhar's results (Chandrasekhar 1935). We find that the Chandrasekhar mass limit is retained irrespective of the luminosity 10 −4  L/L e  10 −2 . However, an increase in surface luminosity leads to progressively higher masses for the larger WDs. This is expected, as larger L translates to more thermal energy, which results in higher ideal gas pressure, thereby allowing the WD to hold more mass. The right panel of Figure 1 shows the radial temperature profiles corresponding to the same luminosities for an R = 10,000 km WD. While the uniform T within the isothermal core increases with a corresponding increase in L, the temperature drops rapidly within the thin nondegenerate surface layers. With the increase in L, the interface shifts inward and the degenerate region shrinks in volume. As expected from Equation (2), the temperature gradient near the surface increases with L. Table 1 lists the central and interface properties for a range of WD radii 2000  R/km  20,000 and luminosities 10 −4  L/L e  10 −2 . It can be seen that an increase in L for a given R leads to somewhat larger central densities ρ c but a more compact degenerate core (smaller interface radius R * ) and therefore an increased capacity to retain more mass. This effect tends to be more pronounced for larger-radius WDs. In the case of smaller stars, the increase in central density due to larger effective thermal energy (or higher L) is not found to be significant such that the Chandrasekhar mass limit is always preserved. For the same L but a larger R, the core temperatures are approximately similar but ρ c decreases substantially.

Magnetic Field Effects
In the previous section, we described the formalism used to self-consistently obtain the properties of WDs by solving the structure equations for a given luminosity and radius and also obtained the solutions for nonmagnetized WDs. Here we discuss in detail the effects of magnetic field on the thermal properties and the mass-radius relations of B-WDs, for both the weak and strong magnetic field limits. Based on the strength of the magnetic field, we appropriately modify our treatment of the radiative opacity, magnetostatic pressure balance, and EOS for the degenerate core and thence the interface location.
We compare the energy density in the magnetic field with that of the gravitational field in order to distinguish between the weakly and strongly magnetized cases. We only evaluate the contribution to the respective energy densities from the degenerate core because the envelope is generally very thin with

Table 1
The Effect of Luminosity L on the Mass-Radius Relation for Nonmagnetized WDs Note. Here M is the total mass and R is the WD radius. The core temperature and density are T c and ρ c , respectively, while ρ * and R * denote the interface density and radius, respectively. a high magnetic field and theoretically vanishing matter density for the cases that we consider here. Assuming a polytropic stellar structure model with an index n = 3 (corresponding to γ = 4/3), the WD central density can be evaluated for a general solution of the Lane-Emden equation as ρ c = M/0.077R 3 ≈ 3.636 ×10 7 g cm −3 (M/1.4M e ) (R/10,000 km) −3 . For a spherically symmetric WD geometry, the average gravitational energy density is ρ m,avg ≈ 3M/4πR 3 6.688 10 g cm where B avg is obtained from the field profile (see Equation (4)) averaged over the matter density assuming ρ 0 ≈ 0.1ρ c . We consider the strong magnetic field limit to be valid provided that ζ = ρ B,avg /ρ m,avg  0.01 for a given WD mass and radius. Figure 2 shows the ratio between the magnetic field and gravitational energies within the isothermal degenerate B-WD core for three different magnetic field configurations and radius within 2000  R/km  20,000. We find that the effect of surface field on the ratio E B /E grav is insignificant provided that the central field is considerably larger than the surface magnetic field. As shown in Figure 2, the gravitational energy dominates over the magnetic field energy for the entire range of WD radius and for central fields as strong as ∼10 14 G. While assuming a spherically symmetric star is justified provided that B  10 14 G, WDs with such strong magnetic fields can often have significantly larger masses as compared to their nonmagnetized counterparts (see Subramanian & Mukhopadhyay 2015). This would then result in a higher gravitational energy density relative to magnetic energy density and consequently marginally smaller E B /E grav ratios than what we estimate for the B = (10 9 , 10 14 ) G case here. Therefore, the B-WD configurations considered in our study will always satisfy the structural stability criteria for both weak and strong magnetic field limits.

Weak-field Limit
For weaker magnetic fields with ζ = 1, the radiative opacity can be effectively approximated with Kramers's formula, κ = κ 0 ρT −3.5 , similar to the nonmagnetized WDs discussed in Section 2. Similarly, the general relativistic effects can also be ignored while considering the magnetostatic pressure balance (see Equation (1)) because the incremental change to the inferred B-WD mass is not found to be significant (Gupta et al. 2020). While the EOS for the B-WD core matter corresponds to that of a nonrelativistic degenerate gas, the surface layer EOS is given by nondegenerate ideal gas. In order to evaluate the interface location, the respective pressures can be equated on both sides to obtain (Shapiro & Teukolsky 1983) where μ e ≈ 2 is the mean molecular weight per electron and ρ * and T * are the density and temperature, respectively, for the interface between the degenerate core and the nondegenerate envelope.

Strong-field Limit
Strong magnetic fields can modify the radiative opacity for photon diffusion and the EOS of the matter therein. For sufficiently large magnetic fields, the variation of radiative opacity with B can be modeled similarly to NSs as κ = κ B ≈ 5.5 × 10 31 ρT −1.5 B −2 cm 2 g −1 (Potekhin & Yakovlev 2001;Ventura & Potekhin 2001). The magnetic-field-dependent Potekhin's opacity is generally used instead of Kramers's opacity if B/10 12 G T/10 6 K and if radiation dominates over convection, which is valid for the strong B cases that we consider here. The general relativistic effects on the inferred mass-radius relation of B-WDs have already been investigated in detail for various magnetic field and rotational configurations (Das & Mukhopadhyay 2014;Subramanian & Mukhopadhyay 2015). The effect of poloidally and toroidally dominated field configurations is seen as super-Chandrasekhar WDs with generally nonspherical shapes that depend on the magnetic field and geometry.
Here we explore the general relativistic effects on the B-WD with the Tolman-Oppenheimer-Volkoff (TOV) equation (Oppenheimer & Volkoff 1939) 6 d dr G P P c m r r P P c r Gm r r c dP d dP d where P = P deg + P ig is the sum of the electron degeneracy pressure and ideal gas pressure. The remaining terms are defined similarly to Equation (1). However, for simplicity, here we assume the WDs to be approximately spherical in shape and neglect the magnetic tension term that can lead to a potentially anisotropic pressure component. Hence, Equation (6)  In the presence of strong fields ζ  0.01, quantum mechanical effects turn out to be important, and Equation (5) is not strictly valid while obtaining the interface radius (Haensel et al. 2007). After including the Landau quantization effects, the degenerate core EOS depends on the magnetic field, whereas the nondegenerate envelope EOS remains unaffected Figure 2. Comparison between magnetic field energy and gravitational energy within the degenerate core used in order to determine the appropriate field regimes for B-WDs with different radii. The results here are shown for varying magnetic fields B = (10 9 , 10 12 ) G (blue circles), B = (10 9 , 10 13 ) G (red squares), and B = (10 9 , 10 14 ) G (green triangles). (Ventura & Potekhin 2001). The electron pressure can then be equated on both sides to give is the magnetic field at the interface radius r = r * . It has already been shown (Das & Mukhopadhyay 2012) that if B > B c ≈ 4.414 × 10 13 G, the electron Larmor radius becomes comparable to the Compton wavelength and the electron-degenerate matter EOS has to be adequately modified. This can potentially yield super-Chandrasekhar WDs with mass limit M ≈ 2.58 M e .

Results
Here we discuss the results related to the mass-radius relation and WD structure properties for both weak and strong magnetic field cases. Figure 3 shows the effect of magnetic field on the mass-radius relation for B-WDs with surface luminosity L = 10 −4 L e and compares them to the nonmagnetic Chandrasekhar results. We find that the magnetic field affects the mass-radius relation in a manner that is analogous to increasing L (see Figure 1), by shifting the curve toward higher masses for WDs with larger radii. While B s does not have any appreciable effect on the inferred mass, the magnitude of B 0 is vital in deciding the shape of the mass-radius curve. The massradius curves for B 0  10 13 G practically overlap with each other and retain the Chandrasekhar mass limit. However, for strong central fields with B 0 ∼ 10 14 G, we obtain super-Chandrasekhar WDs with masses as large as ∼1.9 M e . Table 2 lists the central and interface properties for different magnetic fields and B-WD radii for fixed luminosity L = 10 −4 L e . While the interface location R * and the interface density ρ * do not vary much with B, the core density ρ c and therefore the total mass M increase significantly with an increase in the magnetic field. The central density for stronger fields is expected to be larger in order to compensate for the additional magnetic pressure, as the relative contribution of magnetic field to the total pressure is larger than that to the total density. With an increase in stellar radius for a given B, ρ c reduces substantially whereas R * increases, resulting in smaller average mass density within the isothermal degenerate core. The left panel of Figure 4 shows the radial variation of the matter density for R = 10,000 km and different magnetic field profiles, while the right panel shows the temperature profiles for the same cases. As listed in Table 2, the interface densities ρ * tend to be roughly similar, but the central density ρ c increases with the effective magnetic field. Although the core temperature T c (primarily determined by L) is unchanged irrespective of the magnetic field, the radial temperature gradient dT/dr within the surface layers for B = (10 9 , 10 14 ) G is approximately half that for the other weaker-field cases.

Luminosity Suppression, Cooling, and Field Decay in B-WDs
We study the effects of strong magnetic fields on the observed luminosities for magnetized WDs, i.e., B-WDs. In particular, we compute the luminosities that are expected for B-WDs in order to obtain mass-radius relations that are similar to those for nonmagnetized WDs. We then discuss the cooling process of these B-WDs and also the decay of large magnetic fields by various dissipative processes such as ohmic dissipation, ambipolar diffusion, and Hall drift that can occur within the typical WD age evolution.

Suppression of Luminosity
Bhattacharya et al. (2018a) investigated the variation of surface luminosity as the magnetic field increases, particularly for fixed interface radii and/or temperatures. The motivation for fixing interface parameters between nonmagnetized and magnetized WDs was to better constrain the individual components (thermal, gravitational, and magnetic) of the conserved total energy for these stars. Here we relax the assumptions of fixed interface parameters between nonmagnetized and magnetized WDs. Nevertheless, in order to ensure structural stability for a B-WD, an increase in magnetic energy density has to be compensated by a corresponding decrease in the thermal energy and hence the luminosity, provided that the gravitational energy is not affected significantly. This effect is especially prominent for B-WDs with larger radii, where the magnetic, thermal, and gravitational energies are comparable to each other. Table 3 lists the initial luminosities (see rows with time t = 0) corresponding to field B = (10 9 , 10 14 ) G that yield masses closest to those obtained for nonmagnetized WDs with radii 2000 R/km 20,000. While a slight decrease in the luminosity (within the observable range) for R  12,000 km WDs leads to masses that are similar to those of their nonmagnetic counterparts, the smaller-radius B-WDs require a substantial drop in their luminosity (well outside the observable range) and still do not achieve masses that are similar to the nonmagnetized WDs. This is expected because the thermal pressure is subdominant in the case of WDs with small radii, and therefore a decrease in the thermal energy (or luminosity) does not significantly affect the total mass. As a result, for stars with 2000 R/km 10,000, even if the luminosity decreases substantially to 10 −16 L/L e 10 −12 , the resulting mass of the B-WD remains larger than its nonmagnetic counterpart. This leads to an extended branch in the mass-radius relation.

Magnetic Field Decay
In a strongly magnetized NS, ambipolar diffusion and beta decays can cause the magnetic energy release that is observed from magnetars. However, the magnetic fields inside an electron-degenerate WD generally undergo decay by ohmic dissipation and Hall drift processes with timescales that are given by Heyl & Kulkarni (1998) where ρ c,n = ρ c /10 n g cm −3 , R 4 = R/10 4 km, T c,7 = T c /10 7 K, B 0,14 = B 0 /10 14 G, and l = l 8 × 10 8 cm is a characteristic length scale of the flux loops through the outer core of the WD. In the case of isolated and cool WDs, theoretical calculations indicate that if the magnetic field is not very strong like that for the B-WD center, it typically decays owing to ohmic dissipation by a factor of two in 10 Gyr (Fontaine et al. 1973;Wendell et al. 1987). Cumming (2002) estimated a lowest-order decay time of 8 t Ohm /Gyr 12 for dipole fields and 4 t Ohm /Gyr 6 for quadrupole fields in the context of accreting WDs.
Ohmic decay is characterized by the induction equation, ∂B/ ∂t = −∇ × (η∇ × B), where η = c 2 /4πσ is the magnetic diffusivity and σ is the electrical conductivity. The field decay timescale can then be written as t Ohm ≈ 4πσL 2 /c 2 , where L is the length scale over which the field varies (Cumming 2002). The electrical conductivity is set by the collisions between the electrons and ions (see, e.g., Yakovlev & Urpin 1980;Itoh et al. 1983;Schatz et al. 1999). Unlike within the degenerate core, the electrical conductivity is dependent on the temperature in the surface layers, where the electrons are nondegenerate. As the B-WD mass increases, the increase in conductivity is offset by the decreasing radius.
Ohmic decay is the dominant field dissipation process for B  10 12 G, while for 10 12 B/G 10 14 the decay occurs via Hall drift and for B  10 14 G the principal decay mechanism is likely to be ambipolar diffusion (Heyl & Kulkarni 1998). We assume that Hall drift dominates within the B-WD degenerate core with B 0 ≈ 10 14 G, while ohmic dissipation is the main decay mechanism for surface fields B s ≈ 10 9 G. Although Equations (8) and (9) are primarily used to model the field decay over time for strongly magnetized NSs, here we appropriately adopt them for typical B-WD structure properties.
The magnetic field decay in magnetars with surface fields between 10 14 and 10 16 G was studied, using an appropriate cooling model by Heyl & Kulkarni (1998) and by solving the decay equation where t Amb denotes the ambipolar diffusion timescale. Because we consider strongly magnetized WDs with central fields of about 10 14 G, comparable to the surface fields of magnetars, we assume that the magnetic fields in magnetars and B-WDs both undergo similar decay mechanisms. Muslimov et al. (1995) showed that Hall drift is not expected to be a direct cause of magnetic field decay in WDs because it conserves the total magnetic energy. However, in the presence of magnetic turbulence, Hall drift can twist the field lines and thence enhance ohmic dissipation (Goldreich & Reisenegger 1992). To model and compare the relative contributions from these different processes, we consider two separate cases: (a) when only ohmic dissipation occurs for both the surface and central magnetic fields; and (b) while B s continues to evolve over t Ohm , Hall drift determines the B 0 evolution until the central field drops to about 10 12 G, below which ohmic dissipation sets in.

B-WD Cooling
The cooling evolution of nonmagnetized WDs has been investigated in detail with theoretical attempts to model the degenerate core as the primary energy source (Mestel 1952;Mestel & Ruderman 1967). The thermal energy is radiated Note. The surface field B s and fiducial field B 0 are defined as in Equation (4).
away gradually over time in the observed luminosity from the surface layers as the star evolves. Tutukov & Yungelson (1996) computed the cooling curves for low-mass WDs, starting as luminous stars, until their crystallization stage after about 10 Gyr. Because most of the electrons occupy the lowest energy states in a degenerate gas, the thermal energy of ions is the only significant energy source that can be radiated. The degenerate electrons in the interior of B-WDs have large mean free path, which leads to high thermal conductivity and therefore uniform temperature. The isothermal interior is covered by the nondegenerate surface layers, which transport the energy flux outward by photon diffusion. The interface location between the nondegenerate surface layers and the degenerate interior is obtained by equating the pressure from both sides, as done in Equations (5) and (7) for the weak-and strong-field cases, respectively. The luminosity is related to the interface temperature as L ∝ T 7/2 /κ (Shapiro & Teukolsky 1983), where the opacity κ depends on the field strength (see Section 3). The uniform interior temperature is then computed for a given luminosity, composition, and mass of the B-WD. The rate at which the thermal energy of ions can be transported to the surface and thence be radiated depends on the specific heat (Shapiro & Teukolsky 1983) and is given by where c v ≈ 3k B /2 is the specific heat at constant volume, m μ is the proton mass, and A is the atomic weight. Given an initial luminosity and temperature T 0 at time t 0 , the final temperature T at time t after cooling is obtained by directly integrating Equation (11) where τ = t − t 0 is the WD age. Although convection can aid faster cooling with a more efficient energy transport, its effect has been shown to be insignificant in a first-order approximation (Lamb & Van Horn 1975;Fontaine & Van Horn 1976). This is due to the fact that convection only influences the cooling time once the base of the convection zone reaches the degenerate thermal energy reservoir and couples the surface to the reservoir. However, this is the case only for significantly lower surface temperatures than we consider here. Tremblay et al. (2015) recently showed that convective energy transfer can be significantly impeded once the magnetic pressure dominates over the thermal pressure. It is important to note that, for simplicity in the calculations here, we have assumed self-similarity of the cooling process over the entire evolution of the WD. However, a more detailed calculation of nonmagnetized WD cooling has shown that this might not strictly be the case (Hansen 1999).
In the case of magnetized WDs, the relevant parameter affecting the state of the ionic core and therefore its thermodynamic properties is b = ω B /ω p , where B w =ZeB Mc and p w = Z e n M 4 2 2 p~are the ion cyclotron and plasma frequencies, respectively, Z is the atomic number, M is the mass of the nuclei, n is the ion number density, and e is the electronic charge. The effect of a magnetic field on the ionic core is expected to be strong when the cyclotron frequency is comparable to or larger than the lattice Debye frequency for which b  1. Baiko (2009) studied the effect of magnetic fields on ionic lattices and showed that there is an appreciable change to the specific heat only when b ? 1 except when T = θ D (Debye temperature). For the magnetic field configurations that we consider here b 1 and the interface temperature is comparable to θ D . Hence, we assume a specific heat model that is the same as that for nonmagnetized WDs despite the magnetic fields. The effect of magnetic fields on the phonon spectrum of ions in conventional systems has been studied and found to be generally weak (Holz 1972). Table 3 lists the luminosities and masses corresponding to WD radii 2000  R/km  20,000 for initial B = (10 9 , 10 14 ) G at time t = 0 and t = 10 Gyr. The top entry for each radius is for time t = 0, and the bottom two entries list the corresponding quantities after cooling for t = 10 Gyr, also accounting for simultaneous magnetic field dissipation. We consider two possibilities for the field decay: (a) ohmic dissipation dominates for the entire evolution, and (b) Hall drift is the primary field decay mechanism until B 0 drops below 10 12 G, when ohmic dissipation sets in. Based on the B-WD structure parameters (see Equations (8) and (9)), we find that the fraction of time t HO /τ, when the Hall drift dominates, falls significantly with increasing stellar radius. Consequently, the magnetic field decays considerably more because the faster ohmic dissipation process turns out to be critical for much of the cooling evolution. As a result of the field decay and simultaneous cooling over 10 Gyr, the luminosity-adjusted limiting B-WD masses are found to be substantially lower than their t = 0 counterparts-particularly for the smaller-radius stars. We find that the mass-radius relations practically merge for R  6000 km WDs. Furthermore, the inferred luminosities are also much less suppressed for the intermediate-radius WDs with 6000  R/km  12,000. Although the limiting masses for small B-WDs (with R ≈ 2000 km) drop to about 1.5 M e compared to 1.9 M e without evolution, the majority of these systems still remain practically hidden throughout their cooling evolution because of their strong fields and correspondingly low L, outside the observable range. Figure 5 shows the effect of the evolution of B-WDs on their mass-radius relations, including both the magnetic field decay and thermal cooling effects. The luminosities are varied with the magnetic fields such that the B-WD masses Note. The initial field at t = 0 is kept fixed at B = (10 9 , 10 14 ) G for all the radii listed here. The topmost entry for each radius is for the initial time t = 0, whereas the bottom two entries list the corresponding parameters for t = τ = 10 Gyr after including the cooling rate and magnetic field decay over time. While we evaluate the magnetic fields assuming that ohmic dissipation is the dominant process for the top entries of τ = 10 Gyr, for the bottom entries we assume that Hall drift is the primary process until the field parameter B 0 decays to ∼10 12 G, below which ohmic dissipation dominates.

Results
can match those obtained for the nonmagnetized WDs. We do not show the results for the Hall drift decay because they essentially overlap with the ohmic dissipation results. For the B = (0, 0) case, we find that the mass-radius relation is shifted slightly more toward the Chandrasekhar result, due to the cooling evolution, and the mass limit remains unchanged. In the case of B = (10 9 , 10 14 ) G, even though the limiting mass ∼1.9 M e at small radius turns out to be much larger than the Chandrasekhar limit of 1.4 M e (also see Table 3), we find that it is lowered considerably to ∼1.5 M e primarily as a result of magnetic field decay and also thermal cooling over t = 10 Gyr. The left panel of Figure 6 shows the radial variation of the matter density for the same cases as shown in Figure 5 and R = 10,000 km, while the right panel shows the temperature profiles for the corresponding cases. We find that the matter density at the core is slightly suppressed in the presence of strong B and also as a result of the evolutionary processes. As the total stellar energy is conserved at t = 0, an increase in the magnetic energy has to be compensated by a similar decrease in the gravitational energy and hence the central density ρ c . Once the field decays and L declines owing to cooling, the central density adjusts itself to be slightly lower in order to balance the loss of magnetic and thermal energies with time. From the right panel of Figure 6, we see that the temperature T c within the degenerate core is smaller for larger fields. This is expected, as the luminosity needs to be suppressed more in order to have the total stellar energy fixed. However, owing to the appreciable field decay over t = 10 Gyr, the luminosity increases by more than an order of magnitude. As the degenerate core volume is greater for magnetized WDs with a larger R * , the temperature gradient dT/dr within the envelope turns out to be considerably smaller for B-WDs than their nonmagnetized counterparts.

STARS Results
In this section, we describe the grid of numerical models that we have produced and analyzed in order to investigate, compare, and validate the analytical results described in Sections 2, 3, and 4. In Section 2, we have described the considerations made to self-consistently obtain the structural properties of magnetized WDs. These were used to produce a set of nonmagnetic results to validate our results in the context of the existing literature. In Section 3, we have described the effects of the magnetic field on the mass-radius relations and the thermal properties of B-WDs. We have also described in detail the differences in our methodology and results between the weak-and strong-field cases. In Section 4, we have described the effects of strong magnetic fields on the observed luminosities of magnetized WDs. In particular, we have described the phenomenon of luminosity suppression, cooling, and field decay in B-WDs.

Implementation and Method
In order to investigate the analytical prescription described in the previous sections, here we explore a set of numerical stellar evolution models using a modified version of the STARS stellar evolution code (Eggleton 1971). The EOS solving subroutine STATEF.F is modified appropriately to include the prescriptions of Gupta et al. (2020) (also see Bhattacharya et al. 2018a, who initiated this venture). This involves computing B at each calculation shell and then computing the magnetic contribution to the pressure and density. For large fields such that B/10 12 G T/10 6 K, the opacity is expected to be dominated by the field-dependent Potekhin's opacity rather than the usual Kramers's opacity. Hence, alongside our usual tabulated opacity κ tab , which includes both OPAL opacities for the envelope (Iglesias & Rogers 1996) and electron conduction (Itoh et al. 1983), we evaluate Potekhin's opacity κ B as discussed in Section 3. The overall opacity is then computed by summing the tabulated opacity and the Potekhin opacity in inverse as 1/κ tot = 1/κ tab + 1/κ B .
In order to compute the magnetic field in our calculations, we choose the profile described by Equation (4). It should be noted that the default field profile used here can be appropriately modified to achieve any desired parameterization within the STARS code. The numerical routine STATEF.F receives as input the density at a given mesh point in the model and computes the magnetic field. The magnetic contribution to the density is then computed as ρ B = B 2 / 8πc 2 and that to the pressure as P B = B 2 / 8π. The magnetic pressure and density contributions are added to the current model mesh point pressure and density, before the subroutine continues to compute the remaining thermodynamic quantities. The κ tab opacity computations are then completed before κ B is computed and added in inverse to the tabulated opacity. This allows for the opacity to be computed in a selfconsistent manner rather than switching from Kramers's opacity, or Itoh et al. (1983) style electron conduction opacity, to κ B at specific B and T.
To model strongly magnetized and super-Chandrasekhar WDs, we use the STARS code to generate a zero-age mainsequence star with M = 3 M e and Z = 0.02. The star is then evolved up to the asymptotic giant branch stage until its carbon-oxygen core has grown to about 0.6 M e . At this point, the chemical evolution of the star is halted and an artificial Figure 5. The effect of magnetic field on the WD luminosity set to match with the nonmagnetized mass-radius relation. The results are shown for B = (0, 0) at initial time t = 0 (green circles), B = (0, 0) at t = 10 Gyr (blue diamonds), B = (10 9 , 10 14 ) G at t = 0 (orange triangles), and B = (10 9 , 10 14 ) G at t = 10 Gyr (magenta crosses). The time evolution results are obtained after including B-WD cooling evolution and field dissipation processes. For the magenta curve, we only show the ohmic dissipation results because the results with Hall drift decay practically overlap with them (see Table 3 for the specific luminosities). mass-loss mechanism is enabled to strip the outer envelope until a CO WD with a thin helium atmosphere is formed. The computation of B is then enabled once the model in question has relaxed and been allowed to proceed along its cooling track. Setting a particularly large B (in terms of a large B s , a large B 0 , or both) typically requires the field parameters to be increased in stages to allow the model to relax. The residual mass can then be directly added or removed in order to produce a model WD of any desired mass.

Model Results
We use the STARS evolution code with the modifications described to create a grid of B-WD models with a range of masses and field parameters. We use our grid of models to investigate qualitatively the B-WD mass-radius relationship at different fields, with the objective of numerically validating our analytical models. In Table 4, we list the central temperature, central density, and radius for each of our fixed mass numerical models and for each set of field configurations. In all cases, the models are allowed to cool until the luminosity has reached L = 10 −4 L e . Furthermore, the stellar composition is held fixed and equivalent to the description in Section 2. Unlike in our analytical models, we do not consider separate core and envelope regions but instead allow for the cooling to occur naturally in the numerical models with no explicit prescription. It is not trivial to produce numerically stable B-WD models with arbitrary field configurations with STARS, so we limit the range of field configurations as opposed to our analytical treatment.
In Table 4, we present a number of trends in T c , ρ c , and R for a range of mass and magnetic field. As expected, the B-WD radius decreases with M, as shown in Figure 3 and Table 2. For models of the same mass, R increases very slightly as a function of the magnetic field, until field parameter B 0 = 10 14 G is reached, at which point R increases significantly, consistent with our analytical expectation from Section 3. Here the central density ρ c is not the density in the core but rather the density attained at the central calculation point in our model. This is equivalent to our core and envelope analytical approach described in Tables 2 and 3. The central temperature T c is also computed in a similar manner to that in Section 3. As expected, ρ c increases rapidly with an increase in M. However, ρ c drops marginally once the threshold field B 0 ≈ 10 14 G is reached. This results from a corresponding reduction in R. The comparison of ρ c for a given R with the analytical models presented in Table 2 clearly indicates that the central density does in fact increase as the field increases, with a sharp rise in the density as the critical B 0 is reached. This is in line with our Figure 6. Left panel: the variation of density is shown as a function of radius for the same cases as in Figure 5 and R = 10,000 km. Right panel: the variation of temperature with radius is shown for the same cases as Figure 5 and R = 10,000 km. The corresponding luminosities are listed in Table 3. Note. Here T c and ρ c represent the temperature and density, respectively, computed at the central calculation point at the time step where the model reaches L = 10 −4 L e on its cooling curve.
earlier expectation that ρ c must increase to compensate for the increased magnetic pressure as the field increases. In Figure 7, we validate the mass-radius relations described by our analytical model in Figure 3. As before, the effect of increasing B on the mass-radius relation is analogous to increasing L. This corresponds to an equivalent star earlier on its WD cooling curve. Rather than specifying a WD radius and then inferring its mass, it is more reasonable to base our numerical grid on fixed masses and then compute R. It should be noted that the numerical models here are not computed for a fixed luminosity but at the time when the models have cooled to the point where the code no longer converges. This is owing to the limitations of the EOS utilized in the STARS code at low temperatures and high densities. The radii of these models at low temperatures are almost entirely independent of luminosity. This is expected, as the thermal pressure support in all of these models is negligible compared to degeneracy pressure support and the magnetic pressure support. The three curves representing B 0 = 0, B 0 = 10 12 G, and B 0 = 10 13 G almost completely overlap throughout. This is due to the fact that the magnetic pressure contribution only becomes comparable to the degeneracy pressure contribution for the highest field strength. This reflects the same trend observed at fixed luminosity in Table 4.
We obtain results that are in good agreement with our analytical formalism, and the magnitude of B 0 dictates the shape of the mass-radius curve. For stronger fields, specifically with a larger B 0 , the mass-radius relation deviates from the zero/low-field relation, with the deviation increasing at larger masses. In particular, for B = (10 7 , 10 14 ) G, we obtain super-Chandrasekhar WDs with limiting mass ∼1.9 M e . As anticipated, the radii inferred from our numerical models are not exactly equal to those computed analytically. The EOS for our numerical models is computed with the standard solver in STARS, and this essentially differs from the purely analytical estimates. We also include the effects of neutrino losses in our numerical models. The neutrino cooling effect may cause nonnegligible energy losses from the cores of very hot and/or dense WD stars. Based on the prescription given by Itoh et al. (1983), we model the neutrino losses that become significant once T 10 7 K and ρ 10 10 g cm −3 in the stellar matter. While no L = 10 −4 L e model is hot enough for these losses to occur, as listed in Table 4, many of these models would have had sufficiently hot cores at other points on their cooling curves for neutrino losses to occur.
Next, we investigate the results presented in Section 3.3, which suggest the possibility of obtaining super-Chandrasekhar WDs provided that the central magnetic field and B 0 are sufficiently large. Using our modified STARS code, we compute the highest stable mass model for a range of field configurations. Our results are summarized in Table 5 and are consistent with the Chandrasekhar mass limit being retained for B 0  10 13 G, while allowing for the existence of super-Chandrasekhar B-WD models for larger B 0 . We consider a surface field B s = 10 7 G for the WD models with B 0 = 10 14 G. Regardless, the limiting mass obtained for (B s , B 0 ) = (10 7 , 10 14 ) G with STARS is in perfect agreement with M ≈ 1.865 M e for (B s , B 0 ) ≈ (10 7−9 , 10 14 ) G from the analytical results presented in Table 2. This supports our earlier finding that B s has no appreciable effect on the mass-radius relation. Table 2 and the left panel of Figure 4 have shown the results for analytical computation of the density variation as a function of B-WD radius. As stated, we expect the central density for WDs with stronger fields to be larger in order to compensate for the additional magnetic pressure. Figure 8 shows the numerical validation of the same trend for two mass models, i.e., for M = 0.55 M e and M = 1.6 M e , with varying field configurations. At first sight, the numerical results appear to be inconsistent with our earlier analytical prediction that ρ c should increase to compensate for increased P B as field increases. However, it should be noted that the analytical computations are performed with a fixed radius, whereas our equivalent numerical models are generated assuming a fixed mass. For the super-Chandrasekhar WD models in the right panel of Figure 8 with M = 1.6 M e , we cannot investigate the lower field cases as in the left panel of that figure, as there are no solutions for lower values of B 0 ; hence, we elect to investigate field configurations with higher values of B 0 instead, in an attempt to elucidate a trend in density profiles for these models. We find that for M = 1.6 M e and L = 10 −2.5 L e the radius and density profiles of the model are very strongly dependent on the value of B 0 . As B 0 increases from 10 14 to 10 14.3 to 10 14.5 G, the radius of the model expands from ≈3500 to ≈6000 to ≈8000 km, while the central density falls considerably from ≈3 × 10 8 g cm −3 to ≈2 × 10 7 g cm −3 over the range of values considered for B 0 . This is good confirmation that for these models the total pressure is dominated by the degeneracy pressure and magnetic pressure, with thermal support being  Table 4, we do not compute these numerical models at a fixed luminosity, but rather at the stage where they have already cooled down to the point such that the code fails to converge further. negligible, as expected, and that the density structure and radii of the models are largely functions of the value of B 0 alone. Hence, a model's radius is a function of its magnetic field parameters and its mass at a fixed luminosity. Therefore, the M = 0.55 M e model with a larger central field has a larger radius and hence a lower mean density. If we compare the central density of the model, rather than just the relative order, our numerical results are indeed consistent with the analytical results for R = 10,000 km WDs, with central density ρ c ≈ 2.2 × 10 6 g cm −3 obtained for B = (10 7 , 10 14 ) G. In relation to the results presented in the right panel of Figure 4, we show the variation of the temperature as a function of the radius in Figure 9. In particular, we choose the same mass models and magnetic field configurations as in Figure 8. The analytical models have suggested that the core temperature is primarily determined by the luminosity and is largely unchanged with variation in magnetic field. In good agreement with this prediction, we find here that the core temperatures of models with masses M = 0.55 M e and M = 1.6 M e are in fact largely unchanged with varying magnetic field. The small difference in central temperatures with magnetic field between these numerical models is a result of the difference in radius and hence in the mean density across the models. In the super-Chandrasekhar case, however, the radii, central temperatures, and central densities of the models are very highly dependent on the field configuration. This is, in part, a consequence of the pressure support no longer being dominated by degeneracy pressure, but rather by both magnetic and degeneracy pressures. Hence, we observe that in the super-Chandrasekhar case the density and temperature profiles, as well as the radii of . Each model has been allowed to cool to a luminosity of L = 10 −2.5 L e . This larger luminosity was required as a result of the simulation's EOS encountering difficulties at lower temperatures. Naturally, a certain minimum value of B 0 is required to produce a super-Chandrasekhar model; hence, we have no models corresponding to the lower field models in the 0.55 M e model. Figure 9. Left panel: the variation of temperature as a function of radius is shown for B-WDs with varying magnetic fields B = (0, 0) G (blue), B = (10 7 , 10 12 ) G (orange), B = (10 7 , 10 13 ) G (green), and B = (10 6 , 10 14 ) G (red). These models are essentially equivalent to the analytical models presented in Figure 4. The mass is fixed at M = 0.55 M e for each model, so the radius varies as a function of the magnetic field, as in Figure 7. For each field configuration, the model has been allowed to cool until it reaches L = 10 −4 L e . Right panel: the variation of temperature as a function of radius is shown as in the left panel, but for a super-Chandrasekhar model mass of M = 1.6 M e with varying magnetic fields B = (10 7 , 10 14 ) G (blue), B = (10 7 , 10 14.3 ) G (orange), and B = (10 7 , 10 14.5 ) G (green). As in Figure 8, each model has been allowed to cool to a luminosity of L = 10 −2.5 L e . the models, are highly dependent on the field structure. Our numerical models further demonstrate that the radial temperature gradient dT/dr within the surface layers of each model falls as the magnetic field increases.

Summary of Simulation Results
We have produced a novel set of modifications to the STARS code based on our magnetic field prescription (see Equation (4)) in order to compute a grid of numerical models of highly magnetized WDs. This methodology has allowed us to validate qualitatively the analytical results, from Section 3 in particular. A qualitative approach is necessary because we cannot trivially fix the B-WD radius within the STARS framework and are rather restricted to fixing the mass to generate models that are essentially analogous to our analytical models. We have determined, as we had analytically, that the effect of the surface field B s is not significant, while the central magnetic field and hence B 0 can significantly affect the B-WD mass-radius relation. This is confirmed because the mass M ≈ 1.87 M e computed analytically for (B s , B 0 ) = (10 9 , 10 14 ) G is in very close agreement with the mass M ≈ 1.89 M e inferred from our numerical models for (B s , B 0 ) = (10 7 , 10 14 ) G, i.e., despite the two-order-of-magnitude difference in the surface magnetic fields between the two cases. In Table 5, we have demonstrated that stable numerical models of highly magnetized super-Chandrasekhar WDs can be created provided that the central magnetic field is sufficiently large.

Summary and Conclusions
Sufficiently strong magnetic fields can alter the EOS of electron-degenerate matter to yield super-Chandrasekhar B-WDs with masses that can be as high as M ≈ 2.6 M e , even in the absence of rapid rotation (see Das & Mukhopadhyay 2012Subramanian & Mukhopadhyay 2015). Besides elevating the limiting mass of WDs, strong fields can also impact the thermal characteristics of the underlying compact star and thereby its observed properties. Bhattacharya et al. (2018a) studied the luminosity suppression in B-WDs with B  10 12 G, assuming that the interface properties are essentially similar to their nonmagnetized counterparts, to demonstrate that their cooling rates are significantly attenuated for such strong magnetic fields. Subsequently, Gupta et al. (2020) improved on this preliminary analytical model by removing the assumption of fixed interface parameters that was initially considered for WDs with a preassigned mass. They further derived the mass-radius relation and the stellar structure properties within both nondegenerate and degenerate regions of the B-WD.
In this paper, we have revisited the physics of luminosity suppression and its effect on the mass-radius relation for highly magnetized WDs. We have included the contributions from the electron-degenerate isothermal core, ideal gas surface layer, and magnetic field to model the B-WD structure properties by solving the magnetostatic equilibrium, photon diffusion, and mass conservation equations. In order to distinguish the strongly magnetized cases from the weakly/nonmagnetized cases, we have appropriately amended our treatment of the radiative opacity, magnetostatic pressure balance, and EOS for the degenerate core. Although an increase in surface luminosity results in higher total mass, especially for larger WDs, we have shown that the Chandrasekhar mass limit is retained for 10 −4 L/L e 10 −2 . The increase in luminosity for a given stellar radius leads to a larger ρ c and more compact degenerate interior and hence an increased capacity to hold more mass. Although B s has negligible effect on the B-WD mass, B 0 affects the shape of the mass-radius relation by shifting the curve toward higher masses for stronger fields for a given radius. In particular, strong fields with B 0 ≈ 10 14 G can raise ρ c significantly and yield super-Chandrasekhar WDs with masses as high as ∼1.87 M e .
We have computed the B-WD luminosity necessary in order to obtain mass-radius relations similar to those for nonmagnetized WDs. Provided that the gravitational energy does not vary significantly, an increase in magnetic energy density needs to be compensated by a corresponding reduction in thermal energy, in order to maintain the structural stability of the B-WD. Nevertheless, we have shown that, even with a significant reduction in their luminosities, the masses of smaller-radius B-WDs are not similar to their nonmagnetized counterparts. In particular, for stars with radii 2000 R/ km 6000, the inferred B-WD mass remains considerably larger, leading to an extended branch in the mass-radius relation, even for highly suppressed luminosities such that 10 −16 L/L e 10 −12 . To model the time variation of the B-WD structure properties, we have also considered the cooling evolution of these stars in the presence of strong-field dissipation processes, particularly ohmic dissipation and Hall drift. As the star gradually evolves over time, its thermal energy is radiated away in the observed luminosity from the surface layers. Owing to primarily field decay and also simultaneous cooling over the typical WD age τ ≈ 10 Gyr, the luminosityadjusted masses turned out to be significantly lower than their initial estimates. Even though the limiting B-WD mass is lowered significantly to 1.5 M e compared to 1.9 M e at the initial time, and corresponding initial field B = (10 9 , 10 14 ) G, the majority of these systems still remain practically hidden throughout their cooling evolution owing to their strong fields and consequently low luminosities.
We have also explored a set of stellar evolution models for B-WDs using a modified version of the STARS code with the objective of numerically validating our analytical formalism. For this, we have appropriately modeled the chemical evolution of the star along its cooling track and included the effects of energy losses due to neutrino cooling from the cores of very hot and/or dense WDs. In validation of our analytical approach, we have found that the limiting mass ∼1.8703 M e obtained with the STARS numerical models is in very good agreement with M ≈ 1.87 M e inferred from the analytical calculations for WDs with strong fields B = (10 6 −10 9 , 10 14 ) G. However, the results presented in this work argue that the young super-Chandrasekhar B-WDs may not sustain very large masses over the course of their entire cooling evolution, and this essentially explains their apparent scarcity even without the difficulty of detection owing to their suppressed luminosities.
We have adopted a considerably more generalized framework to improve on the analytical approach presented in our previous studies (Bhattacharya et al. 2018a;Gupta et al. 2020). In addition to radiative cooling, we have now incorporated the effect of convective energy transport through the photon diffusion equation. To account for the effects of strong magnetic fields on the magnetostatic pressure balance and Landau-quantized electron energy states, we have considered the general relativistic TOV equation and appropriately modified the interface condition. Furthermore, we have included typical neutrino cooling energy losses with STARS and the dissipation of strong magnetic fields over time.
We have nevertheless made some simplifying physical assumptions in order to improve our modeling of the B-WD structure properties. First, we have assumed that the B-WDs are approximately spherical for B 0  10 14 G, as demonstrated by Subramanian & Mukhopadhyay (2015) for toroidally dominated fields. Second, we have adopted a constant radial luminosity profile L r = L for our models, as no H burning or other nuclear fusion reactions occur within the degenerate WD core. Third, we have used a temperature-independent EOS for the electron-degenerate pressure (see Gupta et al. 2020 for physical justification) and uniform B-WD core temperature (only for the analytical model) owing to the high thermal conductivity of the degenerate electrons. Lastly, we have assumed the self-similarity of the cooling process for evolution over the typical WD age of τ ∼ 10 Gyr. Future work can address the effect of more complicated magnetic field geometry on the inferred mass-radius relation and the effect of varying stellar compositions on the B-WD cooling evolution. The data produced in this study will be shared on reasonable request to the authors.