Structure and evolution of ultra-massive white dwarfs in general relativity

Context. Ultra-massive white dwarfs ( M (cid:63) (cid:38) 1 . 05 M (cid:12) ) are of utmost importance in view of the role they play in type Ia supernovae explosions, merger events, the existence of high-magnetic -ﬁeld white dwarfs, and the physical processes in the super asymptotic giant branch phase


Introduction
White dwarf stars are the most common end point of stellar evolution. Therefore, these old stellar remnants contain valuable information on the stellar evolution theory, the kinematics, and the star formation history of our Galaxy, and the ultimate fate of planetary systems (see Winget & Kepler 2008;Althaus et al. 2010;García-Berro & Oswalt 2016;Córsico et al. 2019a, for reviews). Furthermore, given the large densities that characterize the white dwarf interiors, these compact objects are considered reliable cosmic laboratories to study the properties of baryonic matter under extreme physical conditions (Isern et al. 2022). Of all the white dwarfs, of special interest are the so-called ultramassive white dwarfs, which are defined as those with masses larger than ∼ 1.05M ⊙ . Ultra-massive white dwarfs play a key population in our Galaxy, subject to the condition that white dwarf masses do not exceed 1.29 M ⊙ .
In recent years, observations of ultra-massive white dwarfs have been reported in several studies (Mukadam et al. 2004;Nitta et al. 2016;Gianninas et al. 2011;Kleinman et al. 2013;Bours et al. 2015;Kepler et al. 2016;Curd et al. 2017;Kilic et al. 2021;Hollands et al. 2020;Caiazzo et al. 2021;Torres et al. 2022). In particular, Gagné et al. (2018) derived a mass of 1.28 ± 0.08 M ⊙ for the long-known white dwarf GD 50. The number of ultra-massive white dwarfs with mass determinations beyond 1.29 M ⊙ is steadily increasing with current observations. Pshirkov et al. (2020) discovered a rapidly rotating ultra-massive white dwarf, WDJ183202.83+085636.24, with M = 1.33 ± 0.01 M ⊙; meanwhile, Caiazzo et al. (2021) reported the existence of a highly magnetized, rapidly rotating ultra-massive white dwarf, ZTF J190132.9+145808.7, with a mass of ∼ 1.327 − 1.365 M ⊙ . Kilic et al. (2021) studied the most massive white dwarfs in the solar neighborhood and concluded that other 22 white dwarfs could also have masses over 1.29 M ⊙ , if they had pure H envelopes and CO cores. Furthermore, Scholz (2022) confirms the existence of a branch of faint blue white dwarfs in the Gaia color magnitude diagram, with some of them also reported in Kilic et al. (2020), which is mainly composed of ultra-massive white dwarfs more massive than 1.29 M ⊙ .
In addition to all these observations, gravity(g)-mode pulsations have been detected in at least four ultra-massive white dwarfs (Kanaan et al. 1992;Hermes et al. 2013;Curd et al. 2017;Rowan et al. 2019). Although these stars have masses slightly below 1.29 M ⊙ , we expect more massive pulsating white dwarfs to be identified in the coming years with the advent of huge volumes of high-quality photometric data collected by space missions such as the ongoing TESS mission (Ricker et al. 2015) and the Cheops (Moya et al. 2018) mission, as well as the future Plato space telescope (Piotto 2018). This high volume of photometric data is expected to make asteroseismology a promising tool to study the structure and chemical composition of ultramassive white dwarfs (De Gerónimo et al. 2019;Córsico et al. 2019b). In fact, several successful asteroseismological analyses of white dwarfs have been carried out using space data thanks to the Kepler/K2 mission (Borucki et al. 2010;Howell et al. 2014;Córsico 2020) and TESS (Córsico 2022).
The increasing number of detected ultra-massive white dwarfs with masses beyond 1.29 M ⊙ as well as the immediate prospect of detecting pulsating white dwarfs with such masses, demand new appropriate theoretical evolutionary models to analyze them. Recently, Schwab (2021a) has studied the evolution of white dwarfs more massive than 1.29 M ⊙ with the focus on neutrino cooling via the Urca process, showing that this process is important for the age determination of ONe-core white dwarf stars. These models were calculated employing the set of standard equations to solve the stellar structure and evolution under the assumption of Newtonian gravity. However, the importance of general relativity for the structure of the most massive white dwarfs cannot be completely disregarded. This was recently assessed by Carvalho et al. (2018), who solved the general relativistic hydrostatic equilibrium equation for a completely degenerate ideal Fermi electron gas. They demonstrate that for fixed values of total mass, large deviations (up to 50%) in the Newtonian white dwarf radius are expected compared to the general relativistic white dwarf radius. The impact of a nonideal treatment of the electron gas on the equilibrium structure of relativistic white dwarfs was studied by Rotondo et al. (2011) and Mathew & Nandy (2017), who derived the mass-radius relations and critical masses in the general relativity framework for white dwarfs of different core chemical compositions. These studies conclude that general relativistic effects are relevant for the determination of the radius of massive white dwarfs. de Carvalho et al. (2014) and, more recently, Nunes et al. (2021) investigated the general relativity effects in static white dwarf structures of nonideal matter in the case of finite temperature. While de Carvalho et al. (2014) focused their work on the effects of finite temperature on extremely low-mass white dwarfs, Nunes et al. (2021) studied the stability of massive hot white dwarfs against radial oscillations, inverse β−decay, and pycnonulcear reactions. They find that the effect of the temperature is still important for determining the radius of very massive white dwarfs.
Despite several works devoted to the study of the effects of general relativity on the structure of white dwarfs, none have calculated the evolution of such structures. Moreover, in all of the works mentioned above, the white dwarf models are assumed to be composed solely of one chemical element. The exact chemical composition determines both the mass limit of white dwarfs and the nature of the instability (due to general-relativity effects or to β-decays; e.g., Rotondo et al. 2011). In this work, we computed the first set of constant rest-mass ultra-massive ONe white dwarf evolutionary models that fully take into account the effects of general relativity on their structural and evolutionary properties. Furthermore, we considered realistic initial chemical profiles as predicted by the progenitor evolutionary history. We employed the La Plata stellar evolution code, LPCODE, to compute the full evolutionary sequence of 1.29, 1.31, 1.33, 1.35, and 1.369 M ⊙ white dwarfs. The standard equations of stellar structure and evolution solved in this code have been modified to include the effects of general relativity. For comparison purposes, the same sequences have been computed for the Newtonian gravity case. We assessed the resulting cooling times and provide precise time dependent mass-radius relations for relativistic ultra-massive white dwarfs. We also provide magnitudes in Gaia, the Sloan Digital Sky Survey, and Pan-STARRS passbands, using the model atmospheres of Koester (2010); Koester & Kepler (2019). This set of cooling sequences, together with the models calculated in Camisassa et al. (2019) and Camisassa et al. (2022), provide a solid theoretical framework to study the most massive white dwarfs in our Galaxy. This paper is organized as follows. In Sect. 2 we describe the modifications to our code to incorporate the effects of general relativity. In Sect. 3 we detail the main constitutive physics of our white dwarf sequences. Sect. 4 is devoted to the impact of general relativity effects on the relevant evolutionary properties of our massive white dwarfs. In this section we also compare and discuss the predictions of our new white dwarf sequences with observational data of ultra-massive white dwarfs, in particular with the recently reported faint blue branch of ultra-cool and ultra-massive objects revealed by the Gaia space mission. Finally, in Sect. 5 we summarize the main findings of this work.

Equations of stellar structure and evolution in general relativity
Our set of ultra-massive ONe white dwarf evolutionary sequences was computed with the stellar evolution code LPCODE developed by La Plata group, which has been widely used and tested in numerous stellar evolution contexts of low-mass stars and particularly in white dwarf stars (Althaus et al. 2003(Althaus et al. , 2005Salaris et al. 2013;Althaus et al. 2015;Miller Bertolami 2016;Silva Aguirre et al. 2020;Christensen-Dalsgaard et al. 2020). For this work, the stellar structure and evolution equations were modified to include the effects of general relativity, following the formalism given in Thorne (1977). Within this formalism, the fully general relativistic partial differential equations governing the evolution of a spherically symmetric star are presented in a way that they resemble the standard Newtonian equations of stellar structure (Kippenhahn et al. 2012). Specifically, the structure and evolution of the star is specified by the Tolman-Oppenheimer-Volkoff (TOV) equation of hydrostatic equilibrium, the equation of mass distribution, the luminosity equation, and the energy transport equation: where t is the Schwarzschild time coordinate, m is the rest mass inside a radius r or baryonic mass, that is, the mass of one hydrogen atom in its ground state multiplied by the total number of baryons inside r, and ϱ is the density of rest mass. Throughout the cooling process, the total baryonic mass remains constant. c is the speed of light, u is the internal energy per unit mass, and ε ν is the energy lost by neutrino emission per unit mass;H , G , V , and R are the dimensionless general relativistic correction factors, which turn to unity in the Newtonian limit. These factors correspond, respectively, to the enthalpy, gravitational acceleration, volume, and redshift correction factors and are given by where m t is the mass-energy inside a radius r and includes contributions from the rest-mass energy, the internal energy, and the gravitational potential energy, which is negative. ϱ t is the density of total nongravitational mass-energy, and includes the density of rest mass plus contributions from kinetic and potential energy density due to particle interactions (it does not include the gravitational potential energy density), that is, ϱ t = ϱ + (uϱ)/c 2 . Since the internal and gravitational potential energy change during the course of evolution, the stellar mass-energy is not a conserved quantity. Φ is the general relativistic gravitational potential related to the temporal metric coefficient. At variance with the Newtonian case, the gravitational potential appears explicitly in the evolution equations. We note that the TOV hydrostatic equilibrium equation differs markedly from its Newtonian counterpart, providing a steeper pressure gradient. We also note that the presence of V in that equation prevents m t from being larger than rc 2 /2G.
The radiative gradient ∇ rad is given by In Eq. (5), ∇ is the convective temperature gradient, which, in the present work, is given by the solution of the mixing length theory. In ultra-massive white dwarfs, the occurrence of convection is restricted to a very narrow outer layer 1 , which is mostly adiabatic. We followed Thorne (1977) in generalizing the mixing length theory to general relativity. In Eq. (3), we omit the energy generation by nuclear reactions since these are not happening in our models. However, they should be added when taking Urca processes into account.
To solve Eqs.
(1)-(5) we need two additional equations that relate m t and Φ with m. These two equations, which are not required in the Newtonian case, have to be solved simultaneously with Eqs. (1)-(5). These extra equations are given by (see Thorne 1977) ∂m For the boundary conditions, we proceeded as follows. The rest mass, total mass-energy, and radius of the star correspond, respectively, to the values of m, m t , and r at the surface of the star. We denote them as M G is the total gravitational mass, that is, the stellar mass that would be measured by a distant observer, which turns out to be less than the total baryonic mass of the white dwarf. Outer boundary conditions for our evolving models are provided by the integration of and assuming a gray model atmosphere. τ is the optical depth and g t is the "proper" surface gravity of the star (as measured on the surface) corrected by general relativistic effects and given by In addition, the general relativistic metric for space time in the star interior must match to the metric outside created by the  star (Schwarzschild metric). The match requires that Φ satisfies the following surface boundary condition: At the stellar center, m = 0, we have m t = 0, r = 0, and L = 0.

Initial models and input physics
We computed the full evolution of 1.29, 1.31, 1.33, 1.35, and 1.369 M ⊙ white dwarfs assuming the same ONe core abundance distribution for all of them. The adopted core composition corresponds to that of the 1.29 M ⊙ hydrogen-rich white dwarf sequence considered in Camisassa et al. (2019), which has been derived from the evolutionary history of a 10.5 M ⊙ progenitor star (Siess 2010). In this work, we restricted ourselves to ONecore massive white dwarfs, thus extending the range of ONe white dwarf sequences already computed in Camisassa et al. (2019) within the framework of the Newtonian theory of stellar interior. ONe core white dwarfs are expected as a result of semi-degenerate carbon burning during the single evolution of progenitor stars that evolve to the super asymptotic giant branch ONe core composition as a result of off-center carbon burning in the merged remnant when the remnant mass is larger than 1.05 M ⊙ (see Schwab 2021b). In particular, it is thought that a considerable fraction of the massive white dwarf population is formed as a result of stellar mergers (Temmink et al. 2020;Cheng et al. 2020;Torres et al. 2022). We note, however, that the existence of ultra-massive white dwarfs with CO cores resulting from single evolution cannot be discarded (see Althaus et al. 2021b;Wu et al. 2022). The adopted input physics for our relativistic white dwarf models is the same as that in Camisassa et al. (2019). In brief, the equation of state for the low-density regime is that of Magni & Mazzitelli (1979), and that of Segretain et al. (1994) was used for the high-density regime, which takes into account all the important contributions for both the solid and liquid phases. We include neutrino emission for pair, photo, and Bremsstrahlung processes using the rates of Itoh et al. (1996), and those of Haft et al. (1994) for plasma processes. The energetics resulting from crystallization processes in the core was included as in Camisassa et al. (2019), and it is based on the two-component phase diagram of dense ONe mixtures appropriate for massive white dwarf interiors, Medin & Cumming (2010). As shown by Blouin & Daligault (2021), 23 Na and 24 Mg impurities only have a negligible impact on the ONe phase diagram, and the two-component ONe phase diagram can be safely used to assess the energetics resulting from crystallization. We did not consider the energy released by a 22 Ne sedimentation process, since it is negligible in ONe white dwarfs (Camisassa et al. 2021).

General relativity effects on the evolution of massive white dwarfs
Here, we describe the impact of general relativity effects on the relevant properties of our constant rest-mass evolutionary tracks. We begin by examining Fig. 1, which displays the general relativistic correction factors H , G , V , and R (black, blue, red, and pink lines, respectively) in terms of the fractional radius for the 1.29, 1.33, 1.35, and 1.369M ⊙ white dwarf models at log L/L ⊙ = −3. Dashed lines in the bottom right panel illustrate the run of the same factors for a 1.369M ⊙ white dwarf model at log L/L ⊙ = −0.4 (log T eff = 5). We recall that these factors are unity in the Newtonian limit. As expected, the importance of general relativistic effects increases as the stellar mass is increased. We note that V is unity at the center and attains a maximum value at some inner point in the star. The relativistic factor R decreases toward the center, departing even more from unity; meanwhile, the other factors, G and H , increase toward the center of the star. The behavior of the relativistic correction factors can be traced back to curvature effects, as well as to the fact that the pressure and the internal energy appear as a source for gravity in general relativity. For maintaining hydrostatic equilibrium, then, both density and pressure gradients are steeper than in Newtonian gravity. This makes the factors G and H , which depend directly on density and pressure, to increase towards the center of the star. The relativistic factor V , which can be interpreted as a correction to the volume, would be unity at the center of the star where the volume is zero and increase because of the increasing density in general relativity respect to the density in Newtonian gravity. However, as the departures from the Newtonian case decrease toward the surface of the star, V decreases toward the outside, achieving a maximum value in between. We note that the relativistic factors depend slightly on the effective temperature.
The impact of relativistic effects on the mass-radius relation at two different effective temperatures can be appreciated in Fig. 2. We note that for the most massive white dwarfs, at a given gravitational mass, the radius is markedly smaller in the case that the general relativity effects are taken into account. At a stellar mass of 1.369 M ⊙ the stellar radius becomes only 1050 km, 25% smaller than predicted by the Newtonian treatment (see Table 1). As in the Newtonian case, the effect of finite temperature on the stellar radius is still relevant in very massive white dwarfs. We mention that general relativistic corrections become negligible for stellar masses smaller than ≈ 1.29 M ⊙ . In particular, for stellar masses below that value, the stellar radius results below 2 % smaller when general relativity effects are taken into account .
In our calculations, ONe white dwarfs more massive than 1.369 M ⊙ become gravitationally unstable (which occurs at a given finite central density) with respect to general relativity effects, in agreement with the findings for zero-temperature models reported in Rotondo et al. (2011) for a pure-oxygen white dwarf (1.38024 M ⊙ ) and Mathew & Nandy (2017) for white dwarfs composed of oxygen (1.3849 M ⊙ ) or of neon (1.3788 M ⊙ ), although their values are slightly higher 2 . We mention that for the 1.369 M ⊙ white dwarf model, the central density in the general relativity case reaches 2.11 × 10 10 g cm −3 (see Table 1). Such density is near the density threshold for inverse β−decays. We have not considered that matter inside our white dwarf models may experience instability against the inverse β−decay. O or Ne white dwarfs are expected to become unstable against the inverse β−decay process at a stellar mass near the critical mass resulting from general relativity effects, on the order of 1.37 M ⊙ (see Rotondo et al. 2011;Mathew & Nandy 2017).
The inner profile of rest mass and density of rest mass for the 1.369M ⊙ white dwarf model in the general relativity and Newtonian cases are shown in the upper and bottom panels of Fig. 3, respectively. For such a massive white dwarf model, general relativity effects strongly alter the stellar structure, causing matter to be much more concentrated toward the center of the star and the central density to be larger than in the Newtonian case. The impact remains noticeable toward lower stellar masses, although to a lesser extent, as can be noted for the case of the 1.35M ⊙ white dwarf model shown in the bottom panel of Fig. 3  radial coordinate for the general relativity case differs markedly from that resulting from the Newtonian case. This is shown in Fig. 4 for the 1.369M ⊙ , 1.35M ⊙ , and 1.29M ⊙ white dwarf models. In particular, the gravitational field in the general relativistic case as measured far from the star is given by Clearly, the gravitational field in the most massive of our models is strongly affected by general relativity. In the stellar interior, large differences arise in the gravitational field due to the inclusion of general relativity effects. We note that such differences do not arise from the relativistic correction factors G V 2 (see Fig. 1) to the Newtonian gravitational field g New = Gm/r 2 that appear explicitly in Eq. (18), but from the solution of the relativistic equilibrium instead, which gives a different run for m(r) compared to the Newtonian case.
Additionally, the surface gravity and stellar radius are affected by the effects of general relativity. These quantities are shown in Fig. 5 in terms of the effective temperature for all of our sequences for the general relativity and Newtonian cases, using solid and dashed lines, respectively. In the most massive sequences, general relativity effects markedly alter the surface gravity and stellar radius. In this sense, we infer that general relativity effects lead to a stellar mass value about 0.015M ⊙ smaller for cool white dwarfs with measured surface gravities of log g ≈ 10. The photometric measurements of Kilic et al. (2021) for the radius of the ultra-massive white dwarfs in the solar neighborhood are also plotted in this figure. For the more massive of such white dwarfs, the stellar radius results 2.8-4% smaller when general relativity effects are taken into account.
We note that most of our sequences display a sudden increase in their surface gravity at high effective temperatures. As noted in Camisassa et al. (2019), this is related to the onset of core crystallization (marked with blue filled circles in each sequence depicted in Fig. 5), which modifies the distribution of 16 O and 20 Ne. Specifically, the abundance of 20 Ne increases in the core of the white dwarf as crystallization proceeds, leading to larger Coulomb interactions and hence to denser cores, and, therefore, to higher surface gravities. This behavior can also be regarded as a sudden radius decrease (bottom panel of Fig. 5). In this context, we note that the density increase due to the increase in the core abundance of 20 Ne during crystallization eventually causes ONe white dwarf models with stellar masses larger than ≳ 1.36M ⊙ to become gravitationally unstable against general relativity effects. In order to explore the mass range of stable white dwarfs in the absence of this processes, the 1.369M ⊙ relativistic sequence was computed disregarding the effect of phase separation (but not latent heat) during crystallization.

General relativity effects on the white dwarf cooling times
The cooling properties of the ultra-massive white dwarfs are also markedly altered by general relativity effects, in particular the the most massive ones. This is illustrated in Fig. 6, which compares the cooling times of our models for the general relativity and Newtonian cases, with solid and dashed lines, respectively. The cooling times are set to zero at the beginning of cooling tracks at very high effective temperatures. Gravothermal energy is the main energy source of the white dwarfs, except at very high effective temperatures where energy released during the crystallization process contributes to the budget of the star. As noticed by Camisassa et al. (2021), ultra-massive ONe-core white dwarfs evolve extremely quickly into faint magnitudes. General relativity effects cause ultra-massive white dwarfs to evolve faster than in the Newtonian case at advanced stages of evolution. In particular, the 1.369M ⊙ relativistic sequence reaches log(L/L ⊙ )=-4.5 in only ∼ 0.5 Gyr, in contrast with the ∼ 0.9 Gyr needed in the Newtonian case. The larger internal densities inflicted by general relativity make the Debye cooling phase more relevant than in the Newtonian case at a given stellar mass, thus resulting in a faster cooling for the sequences that include general relativity effects. The fast cooling of these objects, together with their low luminosity and rare formation rates, would make them hard to observe. The trend in the cooling behavior is reversed at earlier stages of evolution, where white dwarfs computed in the general relativity case evolve more slowly than their Newtonian counterparts. This is because white dwarfs computed in the general relativity case crystallize at higher luminosities (because of their larger central densities), with the consequent increase in the cooling times at those stages. In the 1.369M ⊙ relativistic sequence, the full impact of crystallization on the cooling times is smaller due to the fact that we neglected the process of phase separation during crystallization in that sequence. We mention that we neglected the neutrino emission resulting from Urca process, which is relevant in ONe white dwarfs at densities in excess of 10 9 g cm −3 (Schwab 2021a) . In our modeling, such densities are attained at models with stellar masses of ≳ 1.33M ⊙ (see Table 1). Hence, the depicted cooling times for the sequences with stellar masses above this value may be overestimated at high and intermediate luminosities. A first attempt to include the Urca cooling process from a 23 Na-23 Ne Urca pair in our stellar code leads to the formation of a mixing region below the Urca shell, as reported by Schwab (2021a). Because of the temperature inversion caused by Urca process, our most massive white dwarf models develop off-centered crystallization. We find numerical difficulties in modeling the interaction of crystallization and the Urca process-induced mixing that prevent a consistent computation of white dwarf cooling during these stages. As recently shown by Schwab (2021a), the cooling of such massive white dwarfs is dominated by neutrino cooling via the Urca process during the first 100 Myr after formation. Our focus in this work is on the effects of general relativity on ultra-massive white dwarfs, so we leave the difficult treatment of Urca-process impacts on the structure of relativistic white dwarfs for an upcoming work.

Observational constrains on ultra-massive white dwarf models
The ESA Gaia mission has provided an unprecedented wealth of information about stars (see Gaia Collaboration et al. 2021a, and references therein). In particular, nearly ≈359,000 white dwarf candidates have been detected (Gentile Fusillo et al. 2021), with it being estimated that the sample up to 100 pc from the Sun can practically be considered complete (Jiménez-Esteban et al. 2018). The extreme precision of astrometric and photometric measures allow us to derive accurate color-magnitude diagrams with which to test our models. Some unexpected peculiar features have already been observed in the Gaia white dwarf colormagnitude diagram (Gaia Collaboration et al. 2018), for instance the Q branch, resulting from crystallization and sedimentation delays (Cheng et al. 2019;Tremblay et al. 2019;Camisassa et al. 2021). However, a new branch called the faint blue branch has been reported by Scholz (2022). This faint blue branch is formed by ∼60 ultra-cool and ultra-massive objects, which have been astrometrically and photometrically verified and cross-validated with the Gaia catalog of nearby stars (Gaia Collaboration et al. 2021b) and the white dwarf catalog of Gentile Fusillo et al. (2021). It is also important to mention that some of these objects that form this peculiar feature in the color-magnitude diagram have already been reported, see Kilic et al. (2020) and references therein. Most of these white dwarfs exhibit a near-infrared flux deficit that has been attributed to the effects of molecular collision-induced absorption in mixed hydrogen-helium atmo- spheres, Bergeron et al. (2022). Some issues still remain to be clarified under this assumption and not all the objects in Scholz (2022) are present in the analysis of Bergeron et al. (2022). Consequently, for our purpose here, which is not in contradiction with the analysis of Bergeron et al. (2022), we adopted hydrogen-pure atmosphere models for the analysis of the entire Scholz (2022) sample, where particular objects are treated individually.
In the left panel of Fig. 7, we show a color-magnitude diagram for the 100 pc white dwarf Gaia EDR3 population (gray dots) together with the faint blue branch objects from Scholz (2022) (solid red circles). The color-magnitude diagram selected is absolute magnitude G versus G BP −G, instead of G BP −G RP , thus minimizing the larger errors induced by the G RP filter for faint objects. We also provide the magnitudes for our relativistic and Newtonian models (black and cyan lines, respectively) in Gaia EDR3 passbands (DR2, Sloan Digital Sky Survey, Pan-STARRS and other passbands are also available upon request) using the non-gray model atmospheres of Koester (2010) and Koester & Kepler (2019). Isochrones of 0.25, 0.5, 1, and 2 Gyr for our relativistic model are also shown (dashed black line) in Fig. 7. An initial inspection of the Gaia color-magnitude diagram reveals that our new white dwarf sequences are consistent with most of the ultra-massive white dwarfs within 100 pc of the Sun. In addition, the relativistic white dwarf sequences are fainter than Newtonian sequences with the same mass. Therefore, general relativity effects must be carefully taken into account when determining the mass and stellar properties of the most massive white dwarfs through Gaia photometry. Not considering such effects would lead to an overestimation of their mass and an incorrect estimation of their cooling times. Finally, we check that faint blue branch objects do not follow any particular isochrone, thus ruling out a common temporal origin of these stars. A closer look at the faint blue branch is depicted in the right panel of Fig. 7. The vast majority of faint blue branch white dwarfs appear to have masses larger than ∼ 1.29 M ⊙ . Thus, this sample is ideal for testing our models, and particularly objects that present the largest masses or, equivalently, the smallest radii. Hence, for the analysis presented here and for reasons of com- pleteness we estimated the error bars for objects that lie on the left of the Newtonian 1.369 M ⊙ track. Errors are propagated from the astrometric and photometric errors provided by Gaia EDR3. Although correlations in Gaia photometry are very low, we assumed that some correlation may exist between parameters. This way, errors are added linearly and not in quadrature, thus obtaining an upper limit estimate of the error bars. The parameters corresponding to the 20 selected ultra-massive white dwarf candidates of the faint blue branch are presented in Table  2. In the first column we list the Gaia EDR3 source ID with a label for an easy identification in Fig. 7. Columns 2 to 5 present the parallax, apparent and absolute G magnitudes, and color G BP −G with their corresponding error, respectively. Columns 6 and 7 represent the observational distance within the color-magnitude diagram measured in σ deviations to the limiting 1.369M ⊙ cooling track when the general relativity model or the Newtonian model, respectively, is used. Finally, the last column is a fivedigit number flag. The first digit indicates if the relative flux error in the G BP band is larger than or equal to 10% (1) or smaller (0). The second digit indicates if the relative flux error in the G RP band is larger than or equal to 10% (1) or smaller (0). The third digit indicates if the β parameter as defined by Riello et al. (2021) is ≥ 0.1 (1) or < 0.1 (0); if it is 1 then the object is affected by blending. The fourth digit is set to (0) if the renormalized unit vector ruwe (Lindegren et al. 2018) is < 1.4 (indicative that the solution corresponds to a single object) or set to (1) if it is ≥ 1.4 (bad solution or binary system). The fifth digit indicates if the object passes (1) a 5σ cut on the corrected G BP and G RP flux excess or not (0) (C * ; Riello et al. 2021). An ideal case will show a 00000 flag.
The detailed analysis of the color-magnitude distance to the limiting 1.369M ⊙ relativistic and Newtonian tracks shown in the sixth and seventh columns, respectively, indicates that, on average, the selected faint blue branch objects are more compatible with the general relativistic model than with the Newtonian model. Six of them {a, b, c, f, g, m} lie below the limiting 1.369M ⊙ relativistic track when they are 1σ compatible with the Newtonian model. Moreover, up to four objects {h, j, n, s} are compatible with the relativistic model at the 1σ level but only marginally at a 2σ level with the Newtonian model. In particular, objects { j, s} are ideal candidates to confirm relativistic models given that they present a 00000 flag, which is indicative of a reliable photometry and astrometry. The rest of the objects {d, i, k, o, p, r, t} lie at a distance of 2σ or 3σ (the last two) for the relativistic model but at larger distances for the Newtonian model (up to 4σ). According to our study, these objects with such a small radius or larger masses would be unstable against gravitational collapse. However, any conclusion on this should be taken with caution. On one hand, although some of these objects belong to the sample analyzed by Bergeron et al. (2022; j, J1251+4403, also named WD1248+443 (Harris et al. 2008); o, J1136−1057; and s, J0416−1826) and some nearinfrared flux deficit has been reported for them, a more detailed spectroscopic analysis for all of our candidates is necessary for a precise mass and radius estimation. On the other hand, the presence of strong internal magnetic fields or a rapid rotation, not considered in this paper, could allow these objects to support the enormous gravity. It has been shown, in the general relativity framework, that including strong magnetic fields and/or a rapid rotation could lead to a smaller radius and/or a larger limiting mass for the most massive white dwarfs (for instance Boshkayev et al. 2013;Bera & Bhattacharya 2016;Subramanian & Mukhopadhyay 2015). The existence of super-Chandrasekhar white dwarfs, with masses of 2.1 − 2.8 M ⊙ has indeed been proposed as a possible scenario to explain the over-luminous Type Ia supernovae SN 2003fg, SN 2006gz, SN 2007if, and SN 2009dc (for instance Howell et al. 2006Hicken et al. 2007;Yamanaka et al. 2009;Scalzo et al. 2010;Silverman et al. 2011;Taubenberger et al. 2011). A detailed follow-up of these objects is, in any case, deserved, and, at the same time, general relativistic models such as the ones presented in this work, but for white dwarfs with carbon-oxygen cores are expected to play a significant role in our understanding of the true nature of these objects.

Summary and conclusions
In this paper, we present the first set of constant rest-mass ultramassive ONe white dwarf cooling tracks with masses of M ⋆ > 1.29M ⊙ , which fully take into account the effects of general relativity on their structural and evolutionary properties. Ultramassive white dwarfs are relevant in different astrophysical contexts, such as type Ia supernovae explosions, stellar merger events, and the existence of high magnetic field white dwarfs. In addition, they provide insights into the physical processes in the super asymptotic giant branch phase preceding their formation. In the last few years, the existence of such ultra-massive white dwarfs in the solar neighborhood has been reported in several studies, including the recent discovery of a branch of faint blue white dwarfs in the color-magnitude diagram (Kilic et al. 2020;Scholz 2022). Although some of these objects present an infrared flux deficit, it is also thought to be composed of ultramassive white dwarfs with masses larger than 1.29 M ⊙ . It should be noted that it is very likely that g-mode pulsating ultra-massive white dwarfs with masses of M ⋆ ≳ 1.29M ⊙ will soon be discovered thanks to space missions such as theTESS and Plato space telescopes, and it will then be possible to study them through asteroseismology.
We computed the complete evolution of 1.29, 1.31, 1.33, 1.35, and 1.369 M ⊙ hydrogen-rich white dwarfs models, assuming an ONe composition for the core. Calculations were performed using the La Plata stellar evolution code, LPCODE, for which the standard equations of stellar structure and evolution have been modified to include the effects of general relativity. To this end, we followed the formalism given in Thorne (1977). Specifically, the fully general relativistic partial differential equations governing the evolution of a spherically symmetric star are solved in a way that they resemble the standard Newtonian equations of stellar structure. For comparison purposes, the same sequences were computed, but for the Newtonian case. Our new white dwarf models include the energy released during the crystallization process, both due to latent heat and the induced chemical redistribution. We provide cooling times and time-dependent mass-radius relations for relativistic ultra-massive white dwarfs. We also provide magnitudes in Gaia, Sloan Digital Sky Survey and Pan-STARRS passbands, using the model atmospheres of Koester (2010); Koester & Kepler (2019). This set of cooling sequences, together with those calculated in Camisassa et al. (2019) and Camisassa et al. (2022) for lower stellar masses than computed here, provide an appropriate theoretical framework to study the most massive white dwarfs in our Galaxy, superseding all existing calculations of such objects.
As expected, we find that the importance of general relativistic effects increases as the stellar mass is increased. According to our calculations, ONe white dwarfs more massive than 1.369 M ⊙ become gravitationally unstable with respect to general relativity effects. When core chemical distribution due to phase separation on crystallization is considered, such instability occurs at somewhat lower stellar masses: ≳ 1.360M ⊙ . For our most massive sequence, the stellar radius becomes 25% smaller than predicted by the Newtonian treatment. The evolutionary properties of our ultra-massive white dwarfs are also modified by general relativity effects. In particular, at advanced stages of evolution, the cooling times for our most massive white dwarf sequence are about a factor of two shorter than in the Newtonian case. In addition, not considering general relativity effects when estimating the properties of such objects through photometric and spectroscopic techniques would lead to an overestimation of their mass of 0.015M ⊙ near the critical mass.
In the color-magnitude diagram, we compare our theoretical sequences with the white dwarfs composing the faint blue white dwarf branch (Scholz 2022). We conclude that, regardless of the infrared deficit flux that some particular objects may exhibit, several white dwarfs of this branch can present masses larger than ∼ 1.29M ⊙ and that it does not coincide with any isochrone nor with any evolutionary track. We found that seven of the white dwarfs in this branch should have a smaller radius than our most massive cooling sequence and should be gravitationally unstable against collapse. However, apart from the need for a more detailed spectroscopic study to accurately characterize the possible effects of the infrared flux deficit in some of these objects, the presence of strong magnetic fields and a rapid rotation, not considered in this study, could favor the stability of such objects, thus supporting the existence of super-Chandrasekhar white dwarfs, that, in the case of CO-core white dwarfs, would likely be the progenitors of the over-luminous Type Ia supernovae SN 2003fg, SN 2006gz, SN 2007if, and SN 2009dc. Consequently, a detailed follow-up of these seven objects is required within the framework of the general relativity models exposed here.
As discussed throughout this work, our new ultra-massive white dwarf models for ONe core-chemical composition constitute an improvement over those computed in the framework of the standard Newtonian theory of stellar interiors. Therefore, in support of previous studies, the effect of general relativity must be taken into account to ascertain the true nature of the most massive white dwarfs, and, in particular, to assess their structural and evolutionary properties.  Table 2. Ultra-massive white dwarf candidates selected from the sample of faint blue white dwarfs of Scholz (2022). Sixth and seventh columns indicate the distance within the color diagram of Fig. 7 measured in 1σ deviations form the selected objects to the limiting 1.369 M ⊙ cooling tracks for relativistic and Newtonian models, respectively. Objects j and s, marked in bold, are ideal candidates with no flags to confirm relativistic models. See text for rest of columns and details.