Towards precise constraints in modified gravity: bounds on alternative coupling gravity using white dwarf mass-radius measurements

Tests have to be performed to rule out proposals for gravity modification. We propose a new idea for constraining alternative theories of gravity using temperature-dependent white dwarf (WD) mass-radius (MR) observational data. We have shown that several alternatives to general relativity (GR), which modified GR only within matter, might be reduced to the well-known Poisson equation similar to that of Eddington-inspired Born Infeld (EiBI) and Minimal Exponential Measure (MEMe) gravity. Retaining EiBI notation, we constrain the value of the coupling constant, κ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document}, using a high-precision model-independent measurement of WD MR observations. We have demonstrated that the WD model should include detailed physics to achieve good precision. The model should include their temperature and evolutionary aspects, which may be computationally expensive. To overcome this issue, we construct a semi-analytical surrogate model based on Mestel’s model, calibrated with tabulated, detailed realistic models, to correct the zero-temperature radius. We have shown that the best-fit value of κ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} depends on the WD model, with the ’thick’ envelope models more consistent in describing data. The tightest bound obtained from the most precise MR measurement, QS Vir, with -0.19\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-0.19$$\end{document}≲κ≲\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lesssim \kappa \lesssim $$\end{document} 0.22 in 103\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^3$$\end{document} m5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^5$$\end{document}kg-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{-1}$$\end{document} s-2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{-2}$$\end{document} for 2σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\sigma $$\end{document}(∼95%)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\sim 95\%)$$\end{document} credibility. Overall, we assert that the recent precise WD MR measurements, combined with our current understanding of WD structure, are insufficient to see the deviation from the one predicted by GR. Both more precise observation data and detailed WD modelling are required to rule out gravity modification.


Introduction
Even though Einstein's general relativity (GR) has become the standard for explaining a wide range of phenomena related to gravity, alternative theories or extensions of GR are still widely studied mainly due to several unresolved GR problems [1]. From an observational perspective, GR still cannot precisely explain the dynamics of galactic and intergalactic objects, namely the dark matter (DM) problem. GR also cannot explain the universe's acceleration, namely the dark energy (DE) problem, and recently the tension between some Hubble constant measurements [2]. While from a theoretical perspective, it is challenging to develop a quantum description of gravity within the standard framework. Even though it might be made compatible with quantum field theory as an effective field theory in the low-energy regime.
Modification of gravity usually breaks the underlying assumptions in GR. Some of them are usually designed to address specific problems with GR. For example, Modified Newtonian Dynamics (MOND) [3] was introduced to modify the behavior of gravity in low-density/curvature regions to address the DM problem. Another example is the modification of gravity on short distances, such as in Hořava-Lifshitz gravity [4], which was proposed to address the quantum gravity problem. There are also generalized action theories such as the family of Horndeski theory [5]. However, even if gravity is modified to cure a particular aspect of GR, discrepancies may arise in any other aspects of it. For example, the measurement of the speed of gravitational wave from GW170817 multimessenger observations tightly restricts the value of a T in Horndeski family [6], eliminating a number of theories with exciting features within this family. The above example indicates that applying simultaneous observational constraints with high precision yields tighter bounds on theo-retical parameters. In consequence, many theories might be ruled out.
Even though GR has the above issues, on the other hand, GR is the only theory in almost all fundamental aspects that has been tested through numerous observations and experiments. GR has passed numerous tests from ground to cosmic scale, including recent tests such as the multimessenger observation of GW170817 [7] and the imaging of accretion disks in M87 [8], and Sgr A* [9]. Our understanding of gravity has also been rigorously tested on a Solar System scale [10], leading to the necessity for all theories to be reduced to Newtonian at the non-relativistic, weak gravity field. Extensions in this scale, for example, in terms of the parametrized post-Newtonian (PPN) formalism within GR, have been tightly constrained exceptionally in the vacuum from low to strong field [10][11][12].
Several alternative gravity theories are reduced to GR in a vacuum, making it compatible with most tests in regions of low curvature and low energy. One well-known example is theories with the auxiliary field [13] such as Eddingtoninspired Born Infeld (EiBI) gravity [14], Palatini f (R) [15], and Brans-Dicke theory with ω 0 = − 3 2 . These theories can also be interpreted as GR with the effective energymomentum tensor, as shown in [16]. Another example comes from minimally modified gravity (MMG) [17][18][19], a class of modified gravity that breaks 4-dimensional diffeomorphism invariance. This class of modified gravity is divided into two types. The first consists of the theories which exist in the Einstein frame and can be recast to GR (type-I), while the second consists of the theories that are not in the Einstein frame and cannot be recast to GR plus extra degree(s) of freedom (type-II). In the type-I theory, spacetime and matter are coupled in a non-trivial way. A more specific class for this type, generalized coupling theories (GCTs), couples tensor-energy momentum via a rank 4 coupling tensor in Einstein's field equation. One known example is the minimal exponential measure (MEMe) [20,21], the simplest theory in this class. Since these theories introduce deviations from GR in the presence of matter, it is essential to perform tests to constrain or validate these theories accurately.
The above theories, even though they have different fundamental principles and different relativistic descriptions, have similar characteristics in their Newtonian limit. They reduce to the standard Newtonian theory in a vacuum and conceal the modification of gravity inside the distribution of matter via an additional term in their Poisson equation, which depends on the spatial gradient of matter density. The astrophysical constraints of this Poisson equation have been well-studied in the context of EiBI [22][23][24][25][26]. Similar behaviour is also present in other theories, e.g., the Newtonian limit of beyond-Horndeski theories via the Vainshtein mechanism [27]. The gravitational potential in this theory's spherically symmetric Newtonian solution depends on the second derivative of the mass distribution. As a result, it also relies on the radial gradient of the density.
From a phenomenological perspective, measuring this bound allows us to learn about the parameter values of these theories and helps us understand the matter-gravity coupling behaviour in GR. In the Newtonian limit, we can gain insights into the extent to which gravity depends on the distribution gradient of matter. Although these characteristics are found in many gravity theories, the precise constraint on the mattergradient term, unlike parameters in PPN formalism, has yet to be extensively explored. In fact, by imposing more stringent constraints, we can more effectively rule out the alternative coupling behaviour between matter and gravity. The authors of Ref. [21] argue that this behaviour can be tested in a laboratory experiment and improve the constraint to 10 orders of magnitude tighter than from gravitational wave propagation. However, such constraints should be tested at different scales [28]. Therefore, searching for methods to measure stringent constraints on this term is crucial to provide insights into the coupling behaviour of gravity.
Most astrophysical objects, such as self-gravitating objects, have non-uniform matter distribution so that gravity modification of the corresponding theory will be discernible on these objects. Therefore, the astrophysical objects will be the ideal laboratory to confront the predicted matter structure in modified gravity with observational data. WDs are among the most popular objects to constrain the modified gravity in a non-relativistic scale [26,29]. It has been shown that modified gravity could control WD structure and thus control various aspects of WD, e.g., the MR relation. The parameter constraint of gravity theories can also be deduced using theoretical considerations. For example, we can use perturbative stability, causality bound, or analysis of polytropic model [30]. One can also obtain the constraint by statistical astronomical data analysis, e.g., mass-radius and cataclysmic variability [23]. One of the most popular methods to confront modified gravity with observation is fitting the parameters through WD MR observational data [24,[31][32][33][34][35][36][37].
As modified gravity controls the MR relation of WD, both upper and lower bounds on certain parameters in modified gravity could be determined from the probability range of observational data. Furthermore, studies of modified gravity would provide insight into the phenomenology for both the theory of gravity and the stellar structure model, respectively [33]. This fact implies that the accuracy of these bounds on modified gravity also depends on how accurately the model describes the actual object. However, most studies about the bound on the theory of gravity from WDs use the simplest non-interacting Fermi gas equation of state (EoS). It is known as the Chandrasekhar model [24,31]. It might be useful only for a rough estimation of the range of the parameter. However, one should consider more detailed effects such as WD temperature and their evolutionary history to dive into a more precise range of the parameter.
Several effects have to be considered while modelling WD MR, e.g., Coulomb interaction and other electromagnetic interactions [38], the presence of an envelope and diffusion [39], and most importantly, finite temperature effects [40]. For lower mass WDs (e.g., M 0.8M ), the radius is more sensitive to finite temperature and is higher than their Chandrasekhar radius (e.g., Ref. [40]).
On the other hand, constructing a realistic WD model requires complex calculations that include detailed calculations of the core, envelope, and atmospheric structure and considering their evolutionary history [41]. Unfortunately, in the context of the vast number of alternatives to standard gravity, there still needs to be publicly available code to calculate detailed WD evolution. Besides, the computational cost might arise when these simulation codes are utilized to calculate a statistical distribution with a tiny variation from precise observational data, i.e., it is necessary to recalculate the same procedure for a range of possible values of the parameters to obtain the likelihood distribution. To solve this, Saltas et al. [33] approach the problem by simply correcting the zero-temperature degenerate radius with the difference between their zero-temperature degenerate radius in standard Newtonian gravity to the temperature-dependent WD tabulated model. Saltas et al. studied beyond-Horndeski gravity theory in their work. However, the degenerate core and non-degenerate envelope are expanded differently with temperature. For example, the analysis using Mestel's model [42] for the simplest WD cooling model shows that the core and envelope radius relationship is not strictly linear. Therefore, refining Saltas' approach to achieve better accuracy is crucial.
On the other hand, to select the representative MR data, we must consider their observational method to deduce the result. However, not every observational WD MR data is independent of the structure model, i.e., the well-known MR data from Refs. [43,44] still depend on the evolutionary model calculated within the framework of standard gravity. It means that we need an independent measurement for both mass and radius to have relevant MR data independent of the gravity model. It can usually be done by measuring the radius from luminosity and parallax measurement together with spectral information [45,46]. The more precise measurements may combine several parameters in a single measurement using the Bayesian framework, e.g., measuring dynamical mass and radius from light-curve analysis of an eclipsing binary system and information from stellar photometry and spectral profile [47].
We aim to propose a new approach for constraining the modified gravity parameter using observational mass-radius data. We investigate the bounds of alternative coupling gravity concerning the observational WD MR data by considering the realistic WD models (e.g., temperature effects) and the compatibility between WD models. This paper is organized as follows. The next section will show how general relativity with a non-trivial coupling between spacetime and matter may reduce to the EiBI-like Poisson equation. Then, in Sect. 3, we will explain how we model WDs within this theory and our temperature correction approach. We will also explain the construction of a surrogate model to correct the zero-temperature radius and how we calibrate it with the tabular models. In Sect. 4, we will compare our model with observational data by analysing their likelihood distribution. Finally, we summarize our work in Sect. 5 and conclude in Sect. 6. Note that throughout this paper, we will retain the more familiar notation in EiBI formalism, κ.

General relativity with non-trivial matter coupling
This section will demonstrate how theories that modify gravity inside matter may be reduced to the well-studied EiBIlike Poisson equation in the Newtonian regime. We define the term alternative coupling gravity as a formalism to describe theories that deviate from GR in the presence of matter. This deviation occurs either through a non-trivial spacetimematter coupling assumption or through effective interpretations of an energy-momentum tensor equivalent to such coupling without introducing additional degrees of freedom.
In standard GR, the Einstein tensor G μ ν is directly coupled to the energy-momentum tensor in the matter (or Jordan) space, T with Ricci tensor R μ ν constructed from the connection Γ using the metric q, The matter-side of the above field equation is coupled via the effective tensor T μ ν . Similar to GR, Eq. (1) can be inverted to with T = T α α . With this expression, R μ ν only depends on T . The effective energy-momentum tensor, T μ ν , is typically written solely as a function of the energy-momentum tensor in the matter metric and an arbitrary constant K , . This effective tensor, T μ ν , however, must also be reduced to zero on T μ ν = 0 to revert to GR in vacuum. In some theories, e.g., EiBI and MEMe, it reduces to T μ ν when K = 0. For example, in generalized coupling theory [20], it is written as T μν = χ αβ μν T αβ with χ αβ μν being a coupling tensor which controls the response of spacetime and matter from T αβ . This tensor can be manipulated, or the parameter in that model can take a certain value so that it yields a GR vacuum solution.
As a consequence, if we retain the characteristics above, the metric q μν must be reduced to g μν in vacuum or when a certain parameter reaches a specific value. We can write the relation between these metrices via a deformation matrix, Ω α ν [48], i.e., where Ω α ν has to depend on T , e.g., Ω μ ν (K , T α ν , T δ α ν ) and Ω α ν → δ α ν in vacuum. If we consider T only depends on T , the above equation can be approached by with F α ν depending on the theory or gravity model, which F α ν → 0 in vacuum. Even if Ω α ν matrix is not separable into Eq. (4) (e.g., in MEMe which is proportional to an exponential factor), F α ν can also be written as a correction of δ α ν , e.g., via Taylor series as long as Ω α ν → δ α ν as T α ν → 0. Note that, to retain symmetric q μν , F α ν should be symmetric and contractable with g μα .
To approximate gravity in low curvature and density, we perturb both q μν ≈ η μν + f μν and g μν ≈ η μν + h μν around the Minkowskian metric η μν . Equations (3) and (4) are then approximated as where F (1) μν is the Taylor expansion of F μν up to the first order of T μν . In the Newtonian limit, the pressure component of the energy-momentum tensor is approximately insignificant relative to the energy density ρ, T μν ≈ ρu μ u ν . Then the time component of Eq. (5) becomes with A being a constant that depends on the expansion of F 00 . We define f 00 = −2φ and h 00 = −2φ whereφ and φ are the effective and physical Newtonian potential, respectively. The above equation can now be written as with the constant C = K A. In principle, it is possible to expand F 00 up to an arbitrary order which would lead to a polynomial form of the second term of the above equation, e.g., Cρ + Dρ 2 + Eρ 3 + ..., with D and E denoting certain constants. These constants depend on the exact form of F 00 and hence the expression of T μ ν . Note that this relation only applies to the T which depends only on T (e.g., with no geometrical term except the metric) and reduces to GR in vacuum.
The R 00 component can be approximately reduced to R 00 ≈ 4π Gρ in low curvature and density. Defining R 00 = ∇ 2 φ, now we have the complete Poisson equation of this class of theory In EiBI, C = κ 4 while in MEMe, C = 3q 16π G and the value depends upon the relation between the effective energymomentum tensor T and the physical energy-momentum tensor T (or between the metric q and g).
In the full relativistic description of some theories, the coupling parameter C relates to the cosmological vacuum energy density. In EiBI, it controls the characteristics of the universe's expansion, i.e., the singular behaviour in the early universe [14]. This parameter also affects the propagation speed of gravitational waves in the presence of matter for both EiBI and MEMe [20,49].
Note that from the expression in Eq. (3), the tensor Ω α ν is analogous to the coupling function, A(φ) in scalar-tensor theories [50] which relates the gravitational metric in the Einstein frame and the matter metric in the Jordan frame. Interestingly, in the Newtonian limit, we can approach the Poisson equation in this theory to an EiBI-like term. From Refs. [51][52][53], the field equation in scalar-tensor theory can be expressed as with the scalar field relation while c 0 , c 1 , c 2 and mφ are some constants or parameters. Please refer to the original article (e.g., Ref. [53]) for the exact form of these constants. We can treat 1 − 1 m 2φ ∇ 2 as an operator and approximateφ via Taylor expansion up to the first order of (m −2 φ ) n , then we havē The full Poisson equation in this scalar-tensor theory then can be calculated by inserting Eq. (11) into Eq. (9) to have which is similar to Eq. (8) with parametrized constants. This shows that an EiBI-like additional Poisson term might also exist in other types of theories, regardless of its fundamental assumptions. The discussion regarding the possibility to approach other theories will be discussed elsewhere. This formalism might provide us (at least) with the first-order estimation of how far the spacetime-matter coupling characteristics might deviate from standard Newtonian/GR via the C parameter.
In the next section, we will model the WD structure within this type of Poisson equation before we compare it with observational data. Throughout this paper, we will utilize the EiBI notation, κ, to represent the value of C for convenience and ease of comparison with previous measurements.

Hydrostatic structure in alternative coupling gravity
We model the WD as a static, spherically symmetric star. The star is assumed to be in the hydrostatic equilibrium, which satisfies ∇Φ = − 1 ρ ∇ p. Applying this to Eq. (8), using the EiBI notation, the radial pressure gradient can be written as This equation is similar to the standard Newtonian hydrostatic equilibrium equation with an additional force proportional to the density distribution gradient, dρ dr . With the standard mass-continuity relation, dm(r )/dr = 4πr 2 ρ(r ), and the equation of state (EoS) which relates pressure p and density ρ, one can solve the structure equation to obtain the mass M and radius of the star R. Dividing Eq. (13) with the masscontinuity relation and some manipulations using chain relation, e.g., dρ dr = dρ dp dp dr , the pressure structure equation can also be written in terms of the Lagrangian fluid description (not to be confused with Lagrangian action) as Note that the structure depends on the speed of sound term, ( dρ dp ) −1 . Now, the modification of gravity only depends on the properties of matter. The above equation can be integrated simultaneously with the inverse mass-continuity relation i.e., from m = 0 and p = p c in the center of the star to the surface with m = M and p = 0. The EoS typically depends on temperature T and composition X i , p = p(ρ, T, X i ) to calculate the detailed stellar model. These equations may also be coupled with energy and thermal structure equations for that model. However, we approach the stellar structure for our calculation by calculating temperature-independent con-ditions, i.e., solving the pressure structure and mass continuity equation with an appropriate barotropic EoS, p = p(ρ). We integrate these equations to obtain the WD mass M and degenerate radius R c given κ and p c as the input. The R c is then corrected to find the stellar radius R due to finite temperature with the approach we introduced in Sect. 3.3.
The details on the EoS we choose and how we model the finite-temperature correction will be discussed in the following subsection.

Zero-temperature model
The simplest yet straightforward description of the white dwarf structure is described by the Chandrasekhar model [54] as degenerate matter ruled by the Pauli exclusion principle. In this model, the WD stability is supported only by the degeneracy pressure of electron gas to counterbalance the gravitational force. It implies that the equation of state can be written as the electron's degeneracy pressure relation at zero temperature. In this degenerate matter, electrons can be described as an ideal Fermi gas. The electron density can be determined at zero temperature by integrating the Fermi distribution phase space. The pressure relation can be written as where p F and m e are the Fermi momentum and electron mass, respectively. The density of ions inside the matter dominates the total energy density. The mass density is simply given by ρ = Ax 3 , where the con- 2 3 9.738 × 10 8 μ e kg/m 3 , and it depends on the number of electrons per ion μ e . In the Chandrasekhar model, the EoS depends on the μ e regardless of its atomic number Z . For example, a WD with a He, C, or O core with μ ≈ 2 would have a similar EoS. Therefore, the corresponding atomic abundance yields no effect on the WD structure. The EoS can be solved with the hydrostatic equilibrium and mass conservation equation, as shown in Eqs. (14,15) to obtain the mass and radius of the WD. This solution is sufficient to predict the so-called Chandrasekhar mass limit concept to find the upper bound of a WD mass. This mass limit is an important concept in both WD studies and stellar evolution in general. This principle also had important applications in other fields, e.g., cosmology. The Chandrasekhar limiting mass is connected to Type Ia supernovae as a standard candle to measure the universe's expansion [56].
Even though the Chandrasekhar EoS could approximate several important features of a WD, e.g., the upper limit of mass, to give a more realistic picture of a WD, one should consider the interactions between ions and electrons inside the WD. In this situation, Hamada and Salpeter (HS) [38] modeled additional refinements by assuming that electrons The scattered points are model-independent measurements from [47]. Dark and light gray show limits from [24] for 1σ and 5σ credibility. Note that the GR case is shown in the black line (κ = 0) and ions have a perfect lattice configuration. The Coulomb interaction between electrons and ions has the most significant correction to the Chandrasekhar model, followed by the Thomas-Fermi correction, the electron exchange correction, and the electron correlation correction. These corrections depend on atomic number rather than μ e . These corrections' detailed and explicit expressions are discussed in, e.g., Refs. [57,58]. Overall, the HS correction reduces the typical critical value of a WD mass and radius to around 5%, depending on their exact composition [59].
Even though the zero-temperature model provides a general picture of the WD mass and structure, WDs are hot objects. White dwarfs have been observed to span a wide range of effective surface temperatures (T e f f ) of around 10 3 − 10 4 K , which would consequently give a higher core temperature of up to several orders of magnitude. On the other hand, WDs are not static; for example, they have residual energy from their previous evolutionary sequence. Young WDs tend to have a hot temperature from the pre-WD gravothermal process. A young WD can have a radius up to twice its zero-temperature degenerate radius. As a WD radiates energy, its radius and temperature will be reduced over time [41].
From the observational standpoint, WD spectra usually provide information about their T e f f and surface gravity, g. In principle, one can calculate the mass and radius of a WD solely from g, given a theoretical model, as done in several works of literature, e.g., Refs. [43,44]. However, to confront theoretical predictions related to the structure with observations, both mass and radius must be independently measured from the theoretical model. It requires a combination of measurements. For example, the measurement combination of g with the mass from the astrodynamics of binary stars or alternatively, g with the radius information from luminosity. The more precise measurement usually combines several observational parameters. This procedure has also been used to test WD MR models in standard gravity [47,60]. Figure 1 shows the zero-temperature WD mass-radius for various values of κ compared to the model-independent measurement from Ref. [47,[61][62][63][64][65][66][67][68] and references therein.It shows that neglecting the finite temperature effect would result in a higher value of κ to match the mass and radius values with the observed WDs. The 2σ zero-temperature WD bound of κ from the literature [24] also disagrees with most measurements, especially for lower masses. By visual inspection alone, it is clear that retaining the zero-temperature assumption would result in far higher uncertainty in the predicted mass-radius due to dispersed best-fit values of κ than the standard deviation of observations. This fact shows the necessity to consider the temperature effect when calculating the WD mass and radius to determine the universal bound of κ.
However, the challenge here is calculating the realistic WD within the modified gravitational potential. To be consistent with recent WD mass-radius data with high accuracy, one may need to calculate a detailed physics model of a WD with various values of the parameter κ of the theory. However, as reported in Ref. [33], there is currently no publicly available detailed WD structure evolution program in the context of modified gravities. On the other hand, even with the program in hand, the computational cost would be expensive for multivariate fitting purposes. Hence, the following subsection will explain how we approach the problem by estimating a realistic WD mass radius within modified gravity based on a surrogate model calibrated to the publicly available WD evolutionary grid.

Approach to the realistic WD model
The first step to incorporate the thermal effect into the WD structure is to model a non-zero temperature equation of state on the degenerate matter. The simplest examples are nonzero temperature Fermi gas EoS done in [40] and in [69] for Feynman-Metropolis-Teller EoS. An advanced EoS involves further physics considerations, including quantum and relativity, phase transitions, pair production, and mixing effects. These corrections are discussed in Ref. [70]. It is generally accepted [40,69,71,72] that internal WD temperature may be able to increase the degenerate core radius up to two times its Chandrasekhar radius. The typical temperature of the WD core is approximately three orders higher than its effective temperature, e.g., ∼ 10 6 − 10 8 K [40]. In general, He-core WD is more expanded than CO-core WD for the same mass due to the finite temperature [47,73]. However, to be consistent with WD observation and evolution, one should consider the existence of a non-degenerate envelope surrounding the degenerate sphere.
One of the simplest temperature-dependent WD evolution models was introduced by Mestel [42] by reducing the problem into two aspects [41]. First, the finite temperature effects on the degenerate core are proportional to the total energy content of the WD. Therefore, the degenerate matter in the core tends to be a very effective energy transporter, resulting in an isothermal core. We must be careful here because this assumption is not universal, i.e., it is invalid in very young and hot WDs, where neutrino radiation plays a significant cooling effect. The second is the thermodynamics of the nondegenerate envelope, which seals the internal system of the WD. The envelope seals the internal heat from the outside, imposing a steep heat gradient. In the Mestel model, the envelope is treated as an ideal gas with Kramer's opacity to calculate the radiation transport. This way, the Mestel model can describe the fundamental characteristics of WD evolution, such as cooling times sensitive to core composition, mass, and luminosity.
However, several essential physics on the WD envelope must be correctly included in the Mestel model beyond Kramer's opacity and the adiabatic convection effect. Effectively, Kramer's opacity only accommodates bound-bound, bound-free, and free-free absorption in the atmosphere, which dominates in the ∼ 10 4 − 10 6 K temperature range. Hence, it may fail to predict radiation transport for low T e f f WDs. In addition to radiative transport, another process known as convection can occur in the envelope. In radiative equilibrium, an expanded gas due to higher temperature in the lower layer rises and then shrinks and sinks back due to lower temperature in the upper layer. However, when the temperature of the upper layer is sufficiently high, the expanded gas will keep rising, and convection occurs. As a rule of thumb, convection sets in when the radiative temperature gradient is more significant than its adiabatic counterpart. Overall, the temperature gradient would be less steep if convection was considered.
More realistic models usually solve the structure and its evolution with detailed physical calculations. For example, a model from Ref. [74] computed WDs evolution and synthetic spectra from multiple codes for different aspects (e.g., pre-WD evolution, spectral fitting, and cooling) and had been verified with optical spectral data from SDSS DR12. This model utilized pre-WD profile sequences that evolved from the zero-age main sequence as the initial boundary before it was calculated using a more detailed cooling evolution model to make the track consistent with general stellar evolution.
Another factor that affects the WD radius is the proportion of the surface hydrogen layer to the total WD mass. The observation of the pulsation of ZZ Ceti type WD has given an uncertainty range between M H /M W D = 10 −10 to M H /M W D = 10 −4 [75,76]. However, for the case of the radius of a binary WD-M system as measured in Ref. [47], the "thick" atmosphere (M H /M W D 10 −4 ) realis-tic models are more consistent than their "thin" atmosphere (M H /M W D 10 −10 ) counterpart. We should mention that these studies were reported using Panei's [73] model for Hecore, Fontaine's [77], and Benvenuto's [78] for CO-core. These studies were done assuming a standard Newtonian (non-relativistic limit of GR theory).
The first discussion of modified gravity constraint with finite temperature WD was explained in [33] by comparing mass, radius, and effective temperature within Beyond-Horndeski theories with the data from [47]. They assumed that the radius of low-temperature WD could be approximated with the HS profile. At the same time, the temperature effects were regarded as a correction to the HS radius with , where Y is a free parameter of that theory. The correction, δ R, is taken from the difference between the HS radius at Y = 0 and publicly available WD models constructed within GR/Newtonian by default. However, Mestel's relation shows that the relation between the core radius and the surface radius is not linear, i.e., the expansion of the radius between the degenerate core and the envelope to temperature is not proportional. Hence, the correction has to be treated separately between the core and envelope radius, especially for hot WD. Please see the details of the Mestel model in Appendix A.
Our approach to calculate WD radius We formulate a novel semi-analytic approach as a surrogate of a detailed realistic model by parametrizing Mestel's radiative transport model to find the correction of the zero-temperature radius as a consequence of the finite temperature effect. This method is inspired by Soares's relation [79], which tries to connect the zero-temperature mass-radius with the measured WD mass-radius to obtain temperature-corrected WD, in that case, from the SDSS WD data set [80]. Here, we treat the relation between the isothermal core and envelope radius with the finite temperature treatment based on Mestel's relation. We calibrate the free parameters against the tabulated model's mass, radius, and effective temperature to include beyond-Mestel effects based on several more realistic models.
We define the WD radius R e f f as the photosphere radius where photons escape the surface with effective temperature T e f f . We set the envelope region from the surface to the isothermal core radius R c with temperature T c and apply this assumption to the Mestel model. We assume that the expansion of the isothermal core radius R c due to finite temperature is proportional to its R H S . Similar to [79] we define R c = ξ(T e f f )R H S where ξ(T e f f ) is a temperaturedependent variable which scales the HS radius. The corrected radius, R e f f , is written as The detailed derivation of this approach is shown in Appendix A. To obtain the value of T c , we adopt Koester's relation [81] T 4 which depends on the surface gravity g.
There are two free parameters in this approach, ξ(T e f f ) and μ(T e f f ), which encodes both core and envelope radius correction for a given mass and radius of zero-temperature WD. The μ(T e f f ) represents the envelope, e.g., composition, opacity departure from Kramers, and convection in the atmosphere. In contrast, the ξ(T e f f ) represents the isothermal radius which depends on core composition and temperature. The typical value of μ(T e f f ) should not be far from the value of the number of electrons per ion of ionized H or He (e.g., around 1 or 2), although it will depend on other factors mentioned above. For example, at low-temperature stars (e.g., T < 10 4 K ), the value of μ(T e f f ) would break down since Kramer's law would fail to approach a realistic opacity. Meanwhile, the value of ξ(T e f f ) would be around unity, e.g., the isothermal core radius would not be far from the radius of zero-temperature WD. This way, we approach R e f f from R H S given T e f f , M and the values of free parameters, ξ(T e f f ) and μ(T e f f ), which encodes detailed micro-physics of a realistic model.
We treat these parameters as a numerical approach and by no means interpret the internal structure of a WD. In principle, even though these parameters are derived from and are connected to a physical model, it requires further analyses involving a realistic structure profile to understand the exact physical behaviour of these parameters.
To calibrate our model, we parametrize the mass-radius relation from the above correction relation within standard Newtonian gravity, i.e., with κ = 0 to the mass-radius relation of a realistic WD model. In this way, we obtain the optimum value of ξ(T e f f ) and μ(T e f f ) in a specific value of T e f f . We use the curve_fit function within Scipy [82], a python package for scientific computation. These optimum points for each T e f f in the model grid will form a curve as in Figs. 2 and 3. The ξ(T e f f ) and μ(T e f f ) curves will be a unique representation of a specific WD model. We use these parameters to correct the radius predicted by the HS EoS in modified gravity to be consistent with the detailed WD mass-radius-temperature tabular model. These parameters were fitted into a publicly available tabulated realistic WD evolutionary grid model containing stellar mass, radius or surface gravity, and effective temperature information. We use the model from [74] for CO-core WD (M 0.5M ) for both thick and thin atmosphere, [39] for colder He-core WD (M 0.5M and T e f f 2 × 10 4 K ) and [73] for hotter He-core WD (M 0.5M and T e f f 2 × 10 4 K ). The model from Ref. [39] does not represent a "thin" atmosphere model. The closest representation is only a Fig. 2 The fitted μ and ξ for CO-core models from [74], are shown for both thick (red) and thin (blue) atmosphere models along with their corresponding standard deviation Fig. 3 The fitted μ and ξ for He-core models are shown for both [73] (red line) and [39] (blue dots) atmosphere models along with their corresponding standard deviation. The blue lines are interpolated parameters of available points in the [73] model grid "no" atmosphere, which is irrelevant to our approach. Hence, we only use the "thick" model for He-core WD. The fitting results are shown in Fig. 2 for CO models and Fig. 3 for He models.

Bound of alternative coupling gravity on white dwarf mass-radius relation
To test the theoretical WD MR relation, one should consider independent MR measurements with high precision, which do not rely on the evolutionary or atmospheric model, as the tightness of the κ bound strongly depends on data precision. This can be achieved by measuring photometric, spectral, and orbital parameters of eclipsing binary systems. Therefore, we test our bound using twenty-six precise observational measurements compiled from Ref. [47]. This paper measures 16 WD MR properties in detached WD-M eclipsing binary systems and includes 10 previous measurements to test the theoretical WD MR relation in standard gravity.
We implement the procedure from the previous section to generate values of the stellar properties, M, R, and T e f f of a particular WD model M in a single numerical function. The radius R = R e f f is corrected from R H S from the main structure equations, which are solved using the fourth-order Runge-Kutta method. This numerical function requires p c , T e f f , and κ as the inputs, where the internal composition is enclosed in the WD model M . We match the atomic composition on HS calculation with the fitted WD model, that is, helium for He-core models and carbon for CO-core models. Here, we will show how far the predicted mass-radius would deviate from observational data in the presence of κ.
To do this, we calculate the observables by determining the value of κ, p c and T e f f as the inputs for each M . We interpolate the fitted value of ξ and μ from tabulation using linear interpolation along T e f f .
We calculate the posterior distribution of the optimum value of κ, ln p c and T e f f independently for each object since every observed WD has a unique age, composition, and evolution history. We test different models as explained in the previous section by using the data from Table 2 of Ref. [47]. The best-fit value is determined via Gaussian likelihood L (θ ) [83] ln where θ represents free parameters, and T e f f 2×10 4 K ), we use the He-core model from Ref. [39] with M H /M = 10 −5 . For high temperature (T e f f 2 × 10 4 K ), we use the model from Ref. [73]. For higher mass WDs (M 0.5M ), we use the model from Ref. [74] which can cover up to T e f f ∼ 6 × 10 4 K . We test the model for both types of atmospheres: the "thin" one with M H /M = 10 −10 and the "thick" one with M H /M = 10 −4 . We also fit all objects to a pure HS mass and radius as a comparison to the above results, with carbon and helium composition for masses above and below 0.5M , respectively.
We calculate the probability distribution for each object by generating 24,000 samples of mass-radius structure using emcee [84], a python package for Markov-chain Monte Carlo analyses. We exclude the first 1200 samples in the chain as the burn-in phase to ensure the result is not skewed into the initial value. Here, we include T e f f in the probability distri-  The best-fit value of κ for all objects in [47] is presented with their 1σ credibility. The result for CO-stars (M 0.5M ) is shown in black for the 'thick' envelope and red for the 'thin' envelope of [74]. For He-stars (M 0.5M ), the dotted marks are for a model from [39] (T e f f 2 × 10 4 K ) and crossed marks are for a model from [73] (T e f f 2 × 10 4 K ). Dark and light gray shades show cumulative 1σ and 2σ bounds, respectively. The 'thin' model is excluded from the calculation of the cumulative bound bution calculation to accommodate the uncertainty of measured T e f f ;obs . These generated samples are then marginalized to represent local posterior distribution of parameters, i.e., κ, ln p c , and T e f f and would be different for each object.
The typical posterior distribution of the objects is similar to the ones shown in Fig. 4. It shows the parame- Fig. 6 The comparison between bounds discussed in this work with κ in units of m 5 kg −1 s −2 . The 1σ and 2σ bounds are shown as the inner and outer whiskers, respectively. The most precise measurement is represented by QS Vir, while SDSS J0138-0016 represents the coldest star. Both are calculated using the 'thick' CO-core model. The cumulative bound is obtained by combining He-core and 'thick' CO-core ters' distribution with their predicted observations of QS Vir (M ∼ = 0.78M ) for both the 'thick' and 'thin' envelope models. Both models converge to specific values of p c and κ. The 'thick' envelope model results in a higher p c and lower κ than the 'thin' envelope model, where the distributions of observables overlap between the two models. QS Vir represents the most precise model-independent measurement in this catalogue. Other objects would generally have a similar distribution but with a larger standard deviation. Then we define the bounds on the modified gravity parameter by marginalizing the κ distribution into 1σ or 2σ credibility. Marginalized κ parameters for all twenty-six objects are summarized in Fig. 5 and Table 1. This distribution is similar to the ratio of the observed and model radius.
Approximately, all values of κ in this study (including He-core and CO-core, both 'thin' and 'thick' models) fall within the range of −2 × 10 3 m 5 kg −1 s −2 κ 3 × 10 3 m 5 kg −1 s −2 (see Fig. 5). However, the exact cumulative bound is obtained by combining samples of all objects in this catalogue. To determine the overall bound, we use all He-core results and select the 'thick' model to represent the CO-core model, which is more consistent with other measurements. The value of this bound yields in 10 3 m 5 kg −1 s −2 . The comparison among the bounds with varied WD types is shown in Fig. 6.

Discussion
There are several observational and theoretical motivations for modifying gravity. Even if general relativity is the correct one, studying modifications of gravity will give us lessons to understand the theory itself. Several modified gravity theories retain GR characteristics in a vacuum and behave differently in the presence of matter, even with fundamentally different assumptions. These theories can be recast as Einstein field equations with effective energy-momentum tensors. We discussed the Newtonian limit of this type of theory in Sect. 2.
Considering that this effective energy-momentum tensor is only coupled with the matter, these theories would share a similar Poisson equation with Eddington-inspired Born Infeld gravity in Newtonian scale, which has been extensively studied in the context of compact objects. We investigate how far our understanding of white dwarf structure in this class of theories will be compatible with more recent precise observational data as we test the range of acceptable bounds of the free parameter. We have also shown that the temperature-independent white dwarf model and their modified gravity bound from the previous studies are insufficient to explain the mass-radius distribution of low-mass WD from precise measurement data from [47]. To explain this, we assert that the temperature effect must be addressed, especially to model a low mass WD. However, no publicly available code can calculate WD structure and evolution with detailed physics. Moreover, it requires recurring calculations to estimate the posterior distribution of the parameters concerning data, which requires a high computational cost. Therefore, we construct a surrogate model using the semi-analytical approach to estimate the WD mass-radius in modified gravity based on the Mestel model and Koester's relation. After that, the model is calibrated with the publicly available tabular WD evolutionary model [74] for high mass WD, and Panei et al. model [39,73] for low mass WD. Using this method, we estimate the photosphere radius of a WD given its zero-temperature degenerate mass-radius and effective temperature within the framework of modified gravity, which is now dependent on a free parameter, i.e., κ in EiBI representation.
Then we compare the predicted WD mass-radius obtained using the above method with high-precision mass-radiustemperature data from [47] using Markov-chain Monte Carlo. In this way, we obtain the best-fit probability distribution of κ and other object-specific parameters for each WD. The κ value should remain consistent and universal for all objects and measurements. However, in this study, we found that the observational value of κ will vary proportional to the discrepancy between the predicted and measured WD radius (e.g., in κ = 0), as in Fig. 5. Two aspects can be analysed from the distribution of κ. First, we check the consistency of κ for each WD evolutionary model throughout the mass range. If κ values are dispersed and inconsistent, we expect the model to be incompatible with WD data. It applies even in GR, which requires κ ∼ 0. The second aspect is how far the current most precise measurement could restrict the value of κ.
Our result is summarized in Figs. 5 and 6. WD with mass 0.5M corresponds to CO-core stars, and 0.5M corre-sponds to He-core stars. Our posterior calculation considers the uncertainties from mass, radius, and effective temperature. For CO-core WD, the mass-radius relation from COcore models [74] still undergoes the uncertainty from the thickness value of the envelope. We calculate the probability distribution of κ for both the 'thick' (M H /M = 10 −4 ) and 'thin' (M H /M = 10 −10 ) model of CO-stars and found that from Fig. 5, the distributions of κ for the 'thin' CO-stars distributed dispersedly around 0 < κ < 3 × 10 3 m 5 kg −1 s −2 . This bound is far broader than the individual stars' standard deviation and inconsistent between measurements. On the other hand, the 'thick' CO-stars converge differently around κ ∼ 0 with relative standard deviation with the majority of data. Even with different CO-core tabular models, our result is consistent with the one from [47] that the 'thick' model is more consistent and compatible with GR even in 1σ credibility.
For lower mass WDs ( 0.5M ), due to the availability of tabular models, we only use the 'thick' representation of the model. We use the He-core model from [39] to calculate relatively cold WD, i.e., T e f f 2 × 10 4 K WDs and [73] for hotter ones, i.e., T e f f 2 × 10 4 K WDs. From Fig. 5, the overall value of κ is distributed around κ = 0 and slightly skewed into the negative value. We found that the cumulative κ of He-core WD is still consistent with GR in the 1σ bound. The distribution of κ also showed a similar distribution of R obs /R model as in Fig. 10 of Ref. [47], which shows the consistency in our approach.
In principle, it may be possible to obtain a tighter bound by increasing the observational precision of MR data. However, from Fig. 5, it is apparent that the uncertainty from WD atmosphere thickness significantly affects the consistency of κ. Here, we adopt the 'thick' and 'thin' dichotomy to represent the value of M H as there remains uncertainty in the context of realistic WD modelling. From the astroseismology analysis of ZZ Ceti stars, M H /M is distributed dispersedly. The range is 10 −10 M H /M 10 −4 . It has a different mean with a canonical value of 10 −4 up to two orders (10 −6 ) [75,76]. In our overall result, the 'thick' model (M H /M = 10 −4 , canonical value) gives more consistent κ for each object and is still compatible with GR.
Overall, the bound for most objects coincides with that of GR. It yields a 2σ bound, with SDSS J1212-0123 as the only exception. Most objects even agree with GR (Fig. 6) in the 1σ bound. The coldest WD from the dataset (SDSS J0138-0016), which has a small temperature correction, gives a comparable bound (relative to past result), about −0.71 (−1.08) κ 0.1 (0.47) in 10 3 m 5 kg −1 s −2 with 1σ (2σ ) credibility. The discrepancy of κ calculated between HS and the 'thick' radius on this object is also relatively low compared to other objects in the dataset. The most precise measurement in the dataset, QS Vir, gives a rather tight 1σ (2σ ) bound, −0.09(−0.19) κ 0.17 (0.22) in 10 3 m 5 kg −1 s −2 . We found that this bound is the tightest compared to other WD bounds in the literature.
In the context of the accumulated value of κ, the distribution of individual values is relatively dispersed, resulting in the wide accumulated bound. Even if the individual measurement is reasonably precise, the consistency of individual κ affects the cumulative bound, which reflects how accurately the model predicts WD properties. Our result shows that even though the 'thick' CO-core WDs give more consistent values and tighter bounds of κ than their 'thin' counterparts, the discrepancy between measurements is still apparent. For low-mass WDs, the measurement from He-core WDs gives a relatively loose cumulative bound concerning the standard deviation from individual observation despite having more observations relative to high mass. It shows that the He-model must still be improved to describe low WD mass system as accurately as the 'thick' CO-model explains higher-mass WD. In other words, the bound and accuracy of κ depends on the observation accuracy/precision and how correctly the model depicts the WD structure. The overall bound is then calculated by combining the result from 'thick' CO-core and He-core WDs, which gives 1σ (2σ ) bound of −0.92(−1.99) κ 0.65(1.72) in 10 3 m 5 kg −1 s −2 .
The advantage of using "κ" notation to this gravity model is that we can quickly compare this bound with previous bounds from well-studied EiBI theory in Newtonian scale. One of the earliest bound comes from the solar physics [22], with |κ| < 3 × 10 5 m 5 kg −1 s −2 . From the analysis of contact WD binary Ref. [23] obtained the bound of −1.62 × 10 5 < κ < 1.13 × 10 7 m 5 kg −1 s −2 . Compared to the previous (zero-temperature) WD MR bound from [24] (−0.70 κ 1.66 in 10 3 m 5 kg −1 s −2 ), our overall bound is tighter by a factor of 1.5. The obtained bound is relatively wider because we included the uncertainty of T e f f in our calculations, compared to the measurements using only M and R. Banerjee et al. [24] also defined an upper bound (κ < 0.35 × 10 2 m 5 kg −1 s −2 ) simply by assuming the upper mass of super-Chandrasekhar WD (M ∼ 2.8M ) as the upper bound for Chandrasekhar mass. Note that this bound is obtained for non-relativistic cases. It might deviate significantly from relativistic ones. Another bound that is interesting to compare from brown dwarves investigated by [25], which gives a very tight 1σ constraint of −0.15 κ 0.08 in 10 3 m 5 kg −1 s −2 . This bound is even more stringent than the one from an object with the most precise measurement on the dataset, QS Vir.
It should be noted that our bound may deviate from the exact calculation due to slight inaccuracies in our approach. In GR (κ = 0), for the very high temperature (T e f f ∼ 6 × 10 4 K ), the 'thick' CO-core model deviates up to 1.5% of the radius relative to the tabular model (See Appendix B). The limitation of this approach was explained in Sect. 3.3.
For the exact calculation of the detailed WD model, it may be necessary to simulate the WD cooling evolution from their pre-WD stage within modified gravity to be accurate and self-consistent. This procedure could begin with a simple investigation of the WD cooling model in modified gravity prior to building a more detailed model in the context of modified gravity. For example, the corresponding discussion related to the cooling model within modified gravity can be found in Ref. [85]. On the other hand, the detailed result from this exact structure calculation (e.g., degenerate radius and temperature profile) in modified gravity can be utilized further to calibrate the surrogate model with higher accuracy to accelerate the statistical computation. It is also helpful to test WD MR in a broader context (i.e., not only to test modified gravity) and to quickly calculate WD profile (e.g., for rapid analysis of transient WD detection). This is worth further study.
Despite their rarity, WD in eclipsing binary systems allow us to measure stellar parameters with high precision. Improvements in telescopes or observational methods may enable us to detect more binary WD systems and possibly provide more precise model-independent mass-radius measurements. For example, the large-scale time-domain photometry survey from the Zwicky Transient Facility (ZTF) has discovered more WDs in faint eclipsing binary systems [86] and their follow-up observations have discovered several ZZ Ceti (both confirmed and candidates) in detached eclipsing binary systems [87]. This might enable us to analyse the WD envelope with their high precision mass and radius measurements to rule out envelope thickness uncertainty. Future observations such as the Legacy Survey of Space and Time (LSST) might detect thousands of new WDs in eclipsing binary systems [88] and together with parallax from Gaia, they enable direct characterization of mass and temperature with better than 5 percent uncertainty [89].
Suppose a high-precision measurement of WD massradius has approached the required precision to reveal the necessity to modify gravity. In that case, we still require more detailed physics calculations and accurate WD simulation, particularly for the envelope thickness and structure of low-mass WDs, which are more temperature-sensitive. Moreover, the origin of low-mass WDs is still an open question. It is generally thought that He-core WDs originate from a star that undergoes strong mass-loss episodes during the Red Giant Branch evolution phase. Thus, the accuracy of the He-core WD model will depend on how accurate the evolutionary history simulation of the star is, including the mass transfer of the binary star.
In MEMe, to probe vacuum energy modifications at the TeV scale, one may expect values on the order of |κ| ∼ 10 −33 m 5 kg −1 s −2 . To this end, the controlled null test could be the most robust and straightforward method to constrain this gravity at the Newtonian scale. It can be performed with the modified spherical steel ball experiment, as explained in Refs. [21,90]. This method claims to be able to constrain the parameter up to |κ| 6.75 × 10 −7 m 5 kg −1 s −2 .

Conclusions
Gravity theories with non-trivial matter coupling provide an interesting and useful tool to better understand the coupling behaviour of General Relativity. We have shown that the modification of gravity can be represented by a matter density gradient term in addition to the standard Poisson equation in the Newtonian scale.
Our work indicates that to obtain a precise prediction of the coupling parameter of modified gravity, κ, using the WD MR relation, it is crucial to utilize a realistic, temperaturedependent white dwarf model. This model should be used to calculate both mass and radius with a level of accuracy comparable to the current level of observational precision. By approaching the realistic white dwarf radius through a calibrated surrogate model, we found that the "thick" COcore model is more consistent than its "thin" counterpart and still complies with general relativity. The discrepancies in κ between models indicate that the theoretical white dwarf evolutionary model still needs improvement, particularly in the case of envelope thickness, to comply with precise white dwarf mass and radius measurements.
In conclusion, the current most precise white dwarf measurements and our detailed theoretical model are still insufficient to show deviation from general relativity inside matter. Therefore, the simultaneous improvement in both observational and theoretical aspects of white dwarfs is strongly required to probe the modification of gravity inside matter with high precision.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical work. All data used in this study are from publicly available sources listed in our references.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copy- Fig. 7 The calibrated mass and radius relation of our approach with various temperature for 'thick' (left) and 'thin' (right) CO-core model [74] along with residual between this approach and the corresponding tabulated model used for the callibration (marked with cross mark)

Appendix A: Mestel model
We refer to Mestel's model as the simple WD model introduced in Ref. [42], which describes the thermally isolated degenerate core by the existence of a gas envelope. The derivation of the WD structure within this model can also be found in the literature, e.g., Refs. [41,55]. In this model, the radiative transfer is calculated by photon diffusion relation with Kramer's opacity κ K = κ 0 ρT −3.5 , energy luminosity L r , and emissivity ratio a. Because the envelope only has a tiny portion of mass relative to the stellar mass (m ≈ M), the modified gravity has little to no effect on the envelope structure. We also assume that there is no energy generation and non-radiative dissipation on the surface, so L r = L. The with T tr and R tr are the temperature and radius of the transition zone between core and envelope, while T e f f and R are the effective temperature and radius of the envelope. We use this expression to approach a realistic model by parameterizing μ and R tr = ξ R H S . The explicit expression is shown in Eq. (17). However, in the canonical Mestel model, the surface temperature usually assumed to be T → 0, so the relation is independent with T e f f .