A Stellar Constraint on Eddington-inspired Born–Infeld Gravity from Cataclysmic Variable Binaries

In spite of the unprecedented success of Einstein's general relativity (GR), modifications thereof are important, and have been the focus of intense research over the last few decades. In fact, issues related to the observed cosmic acceleration and the cosmological constant indicate that such modifications are possibly necessary. Another important reason for studying such modifications is the fact that singularities, which are somewhat unpleasant but unavoidable in collapse processes in GR as well as in cosmology, are yet to be understood fully. Indeed, the singularity theorems of Penrose (1965) and Hawking (1966) establish and quantify their formation in terms of geodesic completeness and energy conditions (see, e.g., the textbooks by Hawking & Ellis 1973 and Wald 1984; for a more recent review, see Senovilla 1998). Namely, under certain reasonable assumptions on the energy-momentum tensor, collapse processes in GR lead to geodesic incompleteness, i.e., some classes of geodesics are restricted by an upper bound of an affine parameter, indicating a singular spacetime structure. One implication of this is, for example, the cosmological singularity at the formation of the universe.

Although it is commonly believed that quantum effects might smoothen these singularities, a consistent theory of quantum gravity has been elusive. An alternative then is to construct a singularity-free classical theory of gravity itself, i.e., to modify GR. Indeed, the current status of the literature suggests that in order to achieve geodesic completeness in a spacetime in which some of the scalar quantities diverge, one is possibly forced to consider extensions of GR, see, e.g., the recent works of Olmo et al. (2015), Olmo (2016), and Bazeia et al. (2017). In fact, one can envisage such an extension of GR in line with the celebrated works of Born (1933, 1934) and Born & Infeld (1933), who constructed a viable theory of electromagnetism free of the divergences associated with the more conventional Maxwell form.

One such extension of GR that has been the focus of attention in the recent past is the Eddington-inspired Born–Infeld (EiBI) theory of gravity, put forward by Bañados & Ferreira (2010, hereafter BF), building upon the work of Deser & Gibbons (1998) and Vollick (2004, 2005). Briefly put, BF used an alternative to the Einstein–Hilbert action of GR, proposed by Eddington, where the gauge connection is considered as a fundamental field as opposed to the metric tensor (Eddington 1924; Schrodinger 1985). In Eddington's formalism, the gravitational Lagrangian, apart from prefactors is taken to be $\sqrt{| \det {R}_{\mu \nu }| }$, with R_{μ ν} being the Ricci tensor. Variation of this action with the affine connection being considered as a dynamical variable (dubbed as the Palatini formalism) leads to Einstein's equation in the presence of a cosmological constant. This latter equation is otherwise obtained from a variation of the more conventional Einstein–Hilbert Lagrangian proportional to $\sqrt{| \det {g}_{\mu \nu }| }(R-2{\rm{\Lambda }})$. There, R = g^{μ ν} R_{μ ν} is the Ricci scalar, Λ is the cosmological constant, and in this latter case, the metric g_{μ ν} is the dynamical variable. BF considered a Born–Infeld type of "square root" action, with a minimal coupling of gravity with matter fields, with a Lagrangian proportional to $\sqrt{| \det [{g}_{\mu \nu }+\epsilon {R}_{(\mu \nu )}]| }-\lambda \sqrt{| \det {g}_{\mu \nu }| }$, apart from the matter contribution, and one considers the symmetric part of the Ricci tensor in the Lagrangian, denoted by the braces (for recent work on this, see Beltrán Jiménez & Delhom 2020). Here, 1/ is the Born–Infeld mass M_BI ≪ M_Pl, the Planck mass, and λ is a dimensionless nonzero parameter, related to the cosmological constant, with asymptotically flat solutions corresponding to λ = 1. BF showed that this theory led to singularity-free cosmology. Shortly afterwards, Pani et al. (2011) showed that EiBI gravity is indeed capable of avoiding singularities that arise due to a collapse process, see also Delsate & Steinhoff (2012).

In the nonrelativistic limit, EiBI gravity gives rise to a modified Poisson's equation, with the modification of the low energy limit of Einstein gravity being characterized by a coupling term that is nonzero only in the presence of matter. Since the Poisson's equation is used as a basic input in many of stellar observables, it is then natural that EiBI theories can thus be tested by stellar physics. Indeed, there has been a variety of works in the recent past in this direction. Casanellas et al. (2012) proposed tests for the theory using solar constraints. Avelino (2012b) studied such constraints using cosmological and astrophysical scenarios. Avelino (2012a) obtained bounds on EiBI theories by demanding that electromagnetic forces dominate gravitational ones in nuclear reaction. More recently, Banerjee et al. (2017) have put constraints on the theory from an analysis of white dwarfs. Beltrán Jiménez et al. (2017) studied gravitational waves in nonsingular EiBI cosmological models. Further recent studies on constraining EiBI theories appear in Delhom et al. (2018), Feng et al. (2019), Delhom et al. (2020), and Beltrán Jiménez et al. (2021). A recent comprehensive review of EiBI gravity and related phenomenological tests appear in the work of Beltrán Jiménez et al. (2018).

The starting point of phenomenological studies of EiBI gravity in astrophysical scenarios is the modified Poisson's equation that the theory yields in the Newtonian limit. With the speed of light denoted by c and the gravitational constant by G, expanding the field equations up to first order in , one obtains (Bañados & Ferreira 2010; Beltrán Jiménez et al. 2018),

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}\phi =\displaystyle \frac{8\pi G}{2{c}^{4}}\rho +\displaystyle \frac{8\pi G\epsilon }{4{c}^{4}}{{\rm{\nabla }}}^{2}\rho ,\end{eqnarray} \tag{ 1 }$$

where = 1/M_BI. In units c = G = 1, we will write this equation as

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}\phi =4\pi \rho +\displaystyle \frac{{\kappa }_{g}}{4}{{\rm{\nabla }}}^{2}\rho ,\end{eqnarray} \tag{ 2 }$$

with κ_g = 8π . Observables in stellar physics, in which the Poisson's equation is a crucial input, therefore gets modified in EiBI gravity, and observational data can be used to constrain κ_g . In this paper, we will use data from tidal forces in cataclysmic variable (CV) binary systems (for a comprehensive introduction to CV systems, see Warner 1995) to put such a constraint. As we will discuss in detail in the next section, the fact that the donor star in such a binary fills up its Roche lobe provides us with a way to compute all the observables (like critical mass, radius, etc.) and Equation (2) then implies that these are dependent on the EiBI parameter. Then, using cataloged data on these binaries allows us to put bounds on κ_g . We are able to provide 5σ bounds on this parameter.

Importantly, an attractive feature of Equation (2) is that the modification to the Poisson's equation does not require us to assume spherical symmetry. This is in contrast with many other known theories where such modifications take place. For example, in the beyond-Horndeski class of modified gravity theories (for a recent review, see Kobayashi 2019), only the radial part of the modified pressure balance equation is known, and it is imperative in such examples to assume spherical symmetry of stellar objects. This assumption has to be carefully dealt with, in CV systems where the stellar structure is more complicated due to effects of tidal forces as well as rotations, as discussed by Banerjee et al. (2021). EiBI theories on the other hand provide a much neater picture where one can work explicitly in Cartesian coordinates.

In this paper, we will be interested in studying EiBI gravity in the context of low mass main-sequence stars that are well described by a polytropic equation of state that relates the pressure to the density, and is of the form $P=\kappa {\rho }^{1+\tfrac{1}{n}}$, with κ being the polytropic constant and n the polytropic index. These are the secondaries in CV binary systems with a white dwarf primary, and we will focus on stars of mass ∼0.4 M_⊙. For such low mass stars, the polytropic index n = 1.5 is a good approximation (see, e.g., Rappaport & Joss 1982; Renvoize et al. 2002), a fact that is also borne out by the observation these stars have a rotation time period of less than 6 hours. Here we use the observational data on CV binaries, with the underlying theory being EiBI gravity. This latter fact is inbuilt in our analysis which uses the modified Poisson's equation, Equation (2), and via this, stellar observables are numerically obtained as functions of κ_g , thus offering ready comparison with data, which allows us to constrain the possible values of κ_g . In Section 2, we will briefly elaborate on the numerical recipe that we use. In Section 3, we discuss the methodology of constraining the EiBI parameter κ_g followed by Section 4 which contains our main results. The paper ends with a summary of the results in Section 5.

We will now discuss the procedure to constrain the EiBI parameter using a set of observed data points and comparing it with that obtained numerically. Our numerical procedure involves solving the modified Poisson equation and the Euler equation for hydrostatic equilibrium inside a polytropic star. We find the star's deformed shape under the influence of tidal forces, such that the star is maximally deformed and fills its Roche lobe. Under a stronger tidal field, the star will be tidally disintegrated. Thus, the Roche lobe filling condition gives the critical mass that the star must have, so that it is not tidally disrupted. This is a condition that we can use as a source of information to find the masses and radii of the secondary stars in a CV system. Comparing them with the observed data allows us to constrain the EiBI parameter. This will be the approach that we will follow in this paper. Thus we can consider a polytropic fluid star in flat spacetime throughout its complete trajectory. The effect of the background curvature now comes in the form of an additional force field, namely the tidal field. To begin with, we numerically create a polytropic star, which remains in hydrostatic equilibrium under the influence of a tidal field. The tidal field is calculated in a locally flat Fermi–Normal (FN) frame (Manasse & Misner 1963). The FN frame is particularly useful, as its local flatness allows us to deal with fluid equations in the Newtonian limit, and the local inhomogeneity of gravity is already incorporated in the tidal field.

In the FN frame, the tidal potential can be written as Ishii et al. (2005):

$$\begin{eqnarray}\begin{array}{rcl}{\phi }_{\mathrm{tidal}} & = & \displaystyle \frac{1}{2}{C}_{{ij}}{x}^{i}{x}^{j}+\displaystyle \frac{1}{6}{C}_{{ijk}}{x}^{i}{x}^{j}{x}^{k}\\ & & +\displaystyle \frac{1}{24}\left[{C}_{{ijkl}}+4{C}_{\left({ij}\right.}{C}_{\left.{kl}\right)}-4{B}_{\left({kl}| n| \right.}{B}_{\left.{ij}\right)n}\right]\\ & & \times {x}^{i}{x}^{j}{x}^{k}{x}^{l}+O({x}^{5}).\end{array}\end{eqnarray} \tag{ 3 }$$

Here, the FN coordinates are denoted by xⁱ = {x⁰, x¹, x², x³} and the coefficients are given in terms of the rank 4 Riemann curvature tensor as given by

$$\begin{eqnarray}&&\begin{array}{l}{C}_{{ij}}={R}_{0i0j},\ {C}_{{ijk}}={R}_{0\left(i| 0| j;k\right)},\\ {C}_{{ijkl}}={R}_{0\left(i| 0| j;{kl}\right)},\,{B}_{{ijk}}={R}_{k\left({ij}\right)0},\end{array}\end{eqnarray} \tag{ 4 }$$

where, covariant derivatives are indicated by the symbol ";" and i, j, k, ... run from 1 to 3. The tensorial notation ${R}_{0\left(i| m| j;{kl}\right)}$ denotes a summation over all the possible permutations of i, j, k, and l with m fixed at its position, divided by the total number of such permutations. In the presence of the tidal force field, the deformation is obtained by numerically solving the Euler equation in the FN frame,

$$\begin{eqnarray}&&\begin{array}{l}\rho \displaystyle \frac{\partial {v}_{i}}{\partial \tau }+\rho {v}^{j}\displaystyle \frac{\partial {v}_{i}}{\partial {x}^{j}}=-\displaystyle \frac{\partial P}{\partial {x}^{i}}-\rho \displaystyle \frac{\partial (\phi +{\phi }_{\mathrm{tidal}})}{\partial {x}^{i}}\\ \ \ +\ \rho \left[{v}^{j}\left(\displaystyle \frac{\partial {A}_{j}}{\partial {x}^{i}}-\displaystyle \frac{\partial {A}_{i}}{\partial {x}^{j}}\right)-\displaystyle \frac{\partial {A}_{i}}{\partial \tau }\right].\end{array}\end{eqnarray} \tag{ 5 }$$

Here, as discussed before, P = κ ρ^1+(1/n) is the pressure inside the fluid body of the star. The last term on the right-hand side comes from gravitomagnetic forces where ${A}_{i}=\tfrac{2}{3}{B}_{{ijk}}{x}^{i}{x}^{j}$ is the corresponding vector potential. vⁱ is the velocity field of a fluid element, and ϕ_tidal is the tidal potential as experienced by the star inside the FN frame. The self-gravity of the star is calculated from the modified Poisson equation of Equation (2).

For nonzero κ_g , the second term in Equation (2) quantifies the difference in the self-gravity of the star, as compared to GR. It is important to note that κ_g is expressed in SI units as m⁵ kg⁻¹ s⁻² = [[G]][[R²]], where G is the gravitational constant and R is a length that is usually related to the size of the star. This implies that κ_g will have a higher value as one considers a star of a larger radius, to remain significant in Equation (2) and hence any bound on κ_g will be dependent on the particular star that one considers. A better procedure is then to define the dimensionless quantity ${\bar{\kappa }}_{g}={\kappa }_{g}/({{GR}}^{2})$, which does not depend on the size of the star and therefore can be constrained using any astrophysical object. Here, the only input in our analysis is the polytropic equation of state, and hence the bound that we obtain on the EiBI parameter is universal for all low mass main-sequence stars.

We solve Equations (2) and (5) to find the central density, mass, and volume equivalent radius of the deformed star. The star is corotating in the FN frame with velocity (Ishii et al. 2005; Banerjee et al. 2021)

$$\begin{eqnarray}&&{v}^{i}={\rm{\Omega }}\left[-\{{x}^{3}-{x}_{c}\sin ({\rm{\Omega }}\tau )\},0,\{{x}^{1}-{x}_{c}\cos ({\rm{\Omega }}\tau )\}\right],\end{eqnarray} \tag{ 6 }$$

where Ω is the corotational frequency, and x_c is a constant that arises as the rotational axis deviates from the x²-axis due to the star's deformation. It is convenient to convert these equations into dimensionless form since we do not know the amount of deformation beforehand. First, we rewrite the equations in coordinate ${\tilde{x}}^{i}$, which is defined as the coordinate of the star such that the star has no rotation in the tilde frame and the tilde frame rotates in the FN frame with angular frequency Ω along the x²-axis. Such a choice allows us to remove the τ dependent parts from the above equations. Moreover, we need to write ϕ_tidal and ϕ_mag in terms of the tilde coordinates. Next, we convert the coordinates into dimensionless form ${\tilde{x}}^{i}={{pq}}^{i}$, where p is a constant to be found iteratively as it converges to a prescribed precision and qⁱ is the dimensionless coordinate. Now, Equations (2) and (5) can be written using the dimensionless coordinates as,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{\rm{\Omega }}}^{2}}{2}{p}^{2}\left[{({q}^{1}-{q}_{g})}^{2}+{({q}^{3})}^{2}\right] & = & \kappa (n+1){\rho }^{1+1/n}\\ & & +\phi +{\phi }_{\mathrm{tidal}}+{\phi }_{\mathrm{mag}}+C\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{{p}^{2}}{{\rm{\nabla }}}_{q}^{2}\phi =4\pi \rho +\displaystyle \frac{{\kappa }_{g}}{4{p}^{2}}{{\rm{\nabla }}}_{q}^{2}\rho \end{eqnarray} \tag{ 8 }$$

where ${{\rm{\nabla }}}_{q}^{2}\equiv {p}^{2}{{\rm{\nabla }}}^{2}$ is the dimensionless Laplacian, and in Equation (7) we have used the polytropic equation. As discussed earlier, in units of c = G = 1, κ_g has the dimension of length squared whereas p has the dimension of length, i.e., [[κ_g ]] = [[p²]]. Therefore, the right-hand side of Equation (8) is well defined in dimensionless form. Now, to solve Equations (7) and (8), we need to fix the boundary conditions in order to find the constants C, p, and q_g . These are obtained by fixing the surface of the deformed star at (q_s , 0, 0) where the density is zero. Also, we provide the central density ρ = ρ_c and ∂ρ/∂q¹ = 0 at the center of the star, which is assumed to be the origin of the tilde coordinates. Although the origin may not coincide with the center of the star in the deformed shape, the error is negligible. We solve Equations (7) and (8) together to find the solution in hydrostatic equilibrium. An initial density distribution, along with an initial p is used in Equation (8) to obtain ϕ, which is then used in Equation (7) to find the updated value of p. The updated p is used in Equation (8) and the iteration continues until p converges to a desired precision. This procedure is performed for different central density ρ_c until the Roche lobe filling condition is satisfied. At this critical condition, ∂ρ/∂q¹ at the surface of the star at (q_s , 0, 0) becomes smoothly zero (for more details on the numerical procedure, see Banerjee et al. 2019).

To perform the numerical procedure, we still need to specify the polytropic constant κ, polytropic index n, and the EiBI gravity parameter κ_g . As was discussed earlier, the polytropic index is set to n = 1.5 in our case, whereas, κ_g is varied to put a constraint on it by comparing the numerical results with the observational data. On the other hand, κ can be found by equating the volume equivalent radius at the Roche limit, obtained numerically in GR (κ_g = 0), with the observed radius. The same κ is thereafter used to find the critical mass and radius due to other nonzero κ_g values. Such an assumption is motivated by the fact that the polytropic constant is determined entirely by stellar thermal properties, and thus should remain unchanged for different EiBI parameters. For example, in the case of white dwarfs (WDs), κ is obtained by equating electron degeneracy pressure with the carbon atom density. In our case, we use the Roche lobe filling condition to find κ for which GR is chosen as a reference. Such a procedure is essential when no other ways are known to find the polytropic constant beforehand. However, many other possible values of κ appear because we can use any other radius within the observed range to match the numerical result in GR. We will show in the next section how the best possible choice of the observed radius is made.

As discussed previously, we compare the numerical data with observational data of a set of cataclysmic variable systems, which are binary systems with a Roche lobe filling secondary star orbiting a WD primary. When a star fills its Roche lobe, it is at its critical mass below which it overfills the Roche lobe. The disrupted material is accreted by the primary. This Roche lobe filling condition can therefore be used to find the stellar parameters. We utilize this Roche fill condition to find the unknown polytropic constant κ for which the star's numerical size matches with the observed one.

A set of 13 CV systems is used in this paper. The orbital distances between the primary and the secondary stars in these CVs are large enough to safely neglect the rotation of the primary while calculating the tidal field around the secondary. It allows us to model the gravitational field of the primary as a Schwarzschild geometry. However, the radius of the secondary is about ∼0.1 times the orbital distance, which is too large to generate asymmetry in the deformed shape of the secondary. Hence, we take up to the fourth-order term in the tidal potential. Also, it is safe to assume that the secondary moves in a circular orbit since the type of the orbit does not make any significant difference in the tidal field when the orbital distance is large. Thus the observed orbital parameters of these CV systems can be used to compare with the numerical data for constraining the EiBI gravity parameter. We remind the reader that as already mentioned, the secondary stars in these CVs fall in the main-sequence category and have highly convective cores and small masses ∼0.4M_⊙ as can be interpreted from their orbital period, which is less than 6 hours.

We find the mass of the primary (M₁), mass (M₂) and radius (R₂) of the secondary, and the orbital distance (a) using the Monte Carlo error progression method from a set of observed input parameters such as orbital period (P) and inclination angle (i), mass ratio (q), binary phases at mid-ingress and mid-egress (Δϕ_1/2), radial velocities of the primary and the secondary (K₁ and K₂) and the rotational velocity (vsini) of the secondary star. A detailed description of the procedure can be found in Smith et al. (1998), Thoroughgood et al. (2005), and Horne et al. (1993). A list of the input parameters are shown in Table 1 and the output parameters M₁, M₂, a, and R₂, obtained from Monte Carlo error progression, are given in Table 2.

Table 1. List of Observed Binary Parameters Used for Calculating M₁, a, M₂, and R₂

Name	Observed Binary Parameters
	P	i	q	Δϕ1/2	K1	K2	$v\sin i$	M2
	(hr)	(deg)			(km s−1)	(km s−1)	(km s−1)	(M⊙)
V4140 Sgr a	1.467	80.2 ± 0.5	0.125 ± 0.015	0.0378 ± 0.0005	56 ± 7	⋯	⋯	⋯
V2051 Oph b	1.5	83.3 ± 1.4	0.19 ± 0.03	0.0662 ± 0.0002	91 ± 12	⋯	⋯	⋯
OY Car c	1.51	83.3 ± 0.2	0.102 ± 0.003	0.051 ± 0.004	⋯	470 ± 2.7	⋯	⋯
Ex Hya d	1.638	77 ± 1	⋯	0.017 ± 0.002	69 ± 9	356 ± 4	⋯	⋯
HT Cas e	1.77	81 ± 1	⋯	0.0493 ± 0.0007	58 ± 11	389 ± 4	⋯	⋯
IY Uma f	1.77	86 ± 1	0.125 ± 0.008	0.0637 ± 0.0001	⋯	383 ± 6	⋯	⋯
Z Cha g	1.79	81.78 ± 0.13	0.150 ± 0.004	0.0534 ± 0.0009	⋯	430 ± 16	⋯	⋯
DV Uma h	2.06	84.24 ± 0.07	0.151 ± 0.001	0.063604	⋯	⋯	⋯	0.15 ± 0.02
IP Peg i	3.797	81.8 ± 0.9	0.45 ± 0.04	0.0863	⋯	298 ± 8	⋯	⋯
UU Aqr j	3.93	78 ± 2	⋯	0.051 ± 0.002	121 ± 7	327 ± 31	⋯	⋯
Gy Cnc k	4.211	77 ± 0.9	⋯	0.060 ± 0.005	115 ± 7	283 ± 17	⋯	⋯
Ex Dra l	5.04	${85}_{-2}^{+3}$	0.72 ± 0.06	0.1085 ± 0.0006	⋯	210 ± 14	140 ± 10	⋯
V347 Pup m	5.566	87 ± 3	⋯	0.115 ± 0.005	⋯	198 ± 5	130 ± 5	⋯

Notes.

^a Borges & Baptista (2005), Mukai et al. (1988), and Baptista et al. (1989). ^b Baptista et al. (1998). ^c Wood et al. (1989), Copperwheat et al. (2012), and Littlefair et al. (2008). ^d Hellier et al. (1987), Gilliland (1982), Breysacher & Vogt (1980), Smith et al. (1993), Vande Putte et al. (2003), Hellier (1996), and Ishida et al. (1994). ^e Horne et al. (1991). ^f Steeghs et al. (2003), Patterson et al. (2000), and Rolfe et al. (2002). ^g Wade & Horne (1988), Wood et al. (1986), and Cook & Warner (1984). ^h Feline et al. (2004), Patterson et al. (2000), Mukai et al. (1990), and Henry et al. (1999). ⁱ Martin et al. (1989), Smak (2002), Marsh (1988), and Wood & Crawford (1986). ^j Diaz & Steiner (1991), Vande Putte et al. (2003), and Baptista et al. (1994). ^k Thorstensen (2000). ^l Billington et al. (1996), Baptista et al. (2000), and Fiedler et al. (1997). ^m Thoroughgood et al. (2005), Buckley et al. (1990), and Still et al. (1998).

Download table as:  ASCIITypeset image

We use the observed M₁ and a to find the tidal field in the FN frame of the secondary. Next, we numerically calculate the critical mass (${M}_{2}^{\mathrm{crit}}$) and volume equivalent radius (${R}_{2}^{\mathrm{crit}}$) of the star at Roche limit in the presence of various nonzero values of κ_g . It is found that both ${M}_{2}^{\mathrm{crit}}$ and ${R}_{2}^{\mathrm{crit}}$ increase as κ_g is increased. Hence, we get a range of critical mass and radius of the secondary from the numerical analysis, which is then compared with the observed ranges. In the next section, we analyze the results and find a constraint on the EiBI parameter.

As the EiBI parameter κ_g has a dimension of length squared (in units of G = c = 1), it is evident that the constraint depends on the size of the secondary star R₂. However, we can avoid this limitation by constraining the dimensionless parameter ${\bar{\kappa }}_{g}={\kappa }_{g}/({{GR}}_{2}^{2})$ with any astrophysical object, including the CV secondaries. To numerically find the secondary star's critical mass and radius, we need to calculate the polytropic constant κ. The procedure is as follows:

We choose a value of radius R₂ of the secondary star from its observed range. Now, we perform the numerical procedure keeping κ_g = 0 (i.e., in GR). We can numerically find the volume equivalent radius ${R}_{2}^{\mathrm{crit}}$ of the star at the Roche limit if any value of the polytropic constant κ is given. However, the desired value of κ is found when ${R}_{2}^{\mathrm{crit}}$ becomes equal to R₂. The same κ is then used to numerically calculate the critical masses and radii of the secondary for other nonzero values of the EiBI gravity parameter κ_g . Thus a set of ${M}_{2}^{\mathrm{crit}}$ and ${R}_{2}^{\mathrm{crit}}$ is obtained numerically for various values of ${\bar{\kappa }}_{g}={\kappa }_{g}/{({{GR}}_{2}^{\mathrm{crit}})}^{2}$. A similar numerical procedure is performed for all the CVs.

Now we define a quantity χ² as

$$\begin{eqnarray}&&{\chi }^{2}=\sum _{i=1}^{N}\displaystyle \frac{{\left({M}_{2,i}^{\mathrm{mean}}-{M}_{2,i}^{\mathrm{crit}}\right)}^{2}}{{\sigma }_{M,i}^{2}}+\displaystyle \frac{{\left({R}_{2,i}^{\mathrm{mean}}-{R}_{2,i}^{\mathrm{crit}}\right)}^{2}}{{\sigma }_{R,i}^{2}},\end{eqnarray} \tag{ 9 }$$

where, ${M}_{2,i}^{\mathrm{mean}}$, σ_M,i and ${R}_{2,i}^{\mathrm{mean}}$, σ_R,i are the observed mean and standard deviation of mass and radius of the secondary star of the ith CV system, respectively. On the other hand, ${M}_{2,i}^{\mathrm{crit}}$ and ${R}_{2,i}^{\mathrm{crit}}$ are the numerically calculated critical mass and volume equivalent radius of the Roche lobe filling secondary of the ith CV system, respectively. Here, the total number of systems are taken to be N = 13. We calculate χ² for various values of ${\bar{\kappa }}_{g}$ using the 13 sets of ${M}_{2}^{\mathrm{crit}}$ and ${R}_{2}^{\mathrm{crit}}$ already obtained numerically. Now, to constrain ${\bar{\kappa }}_{g}$ with confidence levels, we need to know the degrees of freedom (dof) of the chi-square test. In our case, dof = 2N − 2 = 24. Finally, we find a constraint on ${\bar{\kappa }}_{g}$ using the χ²/dof versus ${\bar{\kappa }}_{g}$ plot.

However, this constraint is obtained for choosing a particular set of polytropic constants using a particular set of secondary radii from the CV systems. If we take any other set of secondary radii within the observed ranges, another set of polytropic constants as well as a different constraint on ${\bar{\kappa }}_{g}$ is found. The best possible constraint on ${\bar{\kappa }}_{g}$ is obtained only if the best set of secondary radii is chosen, which is done using the chi-square analysis as well. We take different sets of secondary stars' radii as ${R}_{2,i}={R}_{2,i}^{\mathrm{mean}}+\epsilon {\sigma }_{R,i}$, where −1 ≤ ≤ 1. Hence, we find different constraints on ${\bar{\kappa }}_{g}$ for different choices of R_2,i . The best choice of R_2,i is the one for which the minimum of χ²/dof versus ${\bar{\kappa }}_{g}$ plot has the lowest value. Figure 1 shows three plots between χ²/dof and ${\bar{\kappa }}_{g}$ for three different sets of R_2,i . The first one is for ${R}_{2,i}={R}_{2,i}^{\mathrm{mean}}+{\sigma }_{R,i}$ (denoted by red) meaning the maximum values of R_2,i , another is for ${R}_{2,i}={R}_{2,i}^{\mathrm{mean}}-{\sigma }_{R,i}$ (green) indicating the minimum values of R_2,i , and the third one stands for ${R}_{2,i}={R}_{2,i}^{\mathrm{mean}}+({\sigma }_{R,i}/4)$ (blue) for which the χ²/dof versus ${\bar{\kappa }}_{g}$ plot has the lowest minimum. Therefore, ${R}_{2,i}={R}_{2,i}^{\mathrm{mean}}+({\sigma }_{R,i}/4)$ is the best choice for calculating the desired set of κ for each of the CV secondaries and therefore constraining ${\bar{\kappa }}_{g}$. In Figure 2, we show the confidence levels to which ${\bar{\kappa }}_{g}$ is constrained. We find ${\bar{\kappa }}_{g}$ to be $0.005\leqslant {\bar{\kappa }}_{g}\leqslant 0.352$ within 1σ and $-0.315\leqslant {\bar{\kappa }}_{g}\leqslant 0.597$ within the 5σ confidence level.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 3, we show the M₂, R₂ ranges generated by Monte Carlo using the observed parameters as compared to numerically obtained ${M}_{2}^{\mathrm{crit}}$ and ${R}_{2}^{\mathrm{crit}}$, which are obtained using the best choice of κ. Numerical data are denoted by black lines, which are obtained for various values of ${\bar{\kappa }}_{g}$ within its 5σ limit. As already mentioned, both the critical mass and radius increase with ${\bar{\kappa }}_{g}$. It can be seen that numerical results tend toward nonzero positive values of ${\bar{\kappa }}_{g}$ to match the observed parameters as also evident from the χ² analysis.

Zoom In Zoom Out Reset image size

There are two important issues that we will discuss at this stage. First, we note that Banerjee et al. (2021) discussed constraining modified gravity theories of the beyond-Horndeski class. An important difference between the method followed there as compared to the present paper is as follows. In Banerjee et al. (2021), the polytropic constant is not kept fixed while varying the modified gravity parameter. As a result, numerical mass values are influenced by both the polytropic constant as well as the modified gravity parameter. For each value of this parameter, the polytropic constant is freely adjusted until the volume equivalent radius matches the observed radius. While this procedure is correct in its own merit, it suppresses the modified gravity parameter's effects somewhat, and makes it difficult to constrain. A better procedure is to consider the modifications coming solely from the modified gravity parameter. However, eliminating the effects of the polytropic constant requires its value to be known either by some other physical equations (like that in the case of WDs) or by a statistical best choice method. In the case of the CV secondaries, we do not have any other information to find their polytropic constant beforehand. Therefore, here we have resorted to the best choice of the polytropic constant statistically, and we keep it fixed for all values of the modified parameter.

Second, we note from Figure 3 that for a few CV systems such as OY Car, Z Cha, and DV Uma, the numerical data do not fall onto the Monte Carlo generated distributions. Such tension in the data appears from the fact that R₂ is obtained using Δϕ_1/2 which is independent of mass M₂. On the other hand, ${M}_{2}^{\mathrm{crit}}$ is numerically dependent on R₂ since it is matched with the volume equivalent radius of the secondary. Therefore, observed mass ranges can possibly deviate from the numerical data. We, however, note that both ${M}_{2}^{\mathrm{crit}}$ and ${R}_{2}^{\mathrm{crit}}$ individually fall within the observed ranges of M₂ and R₂, respectively, while extending both sides of the observed mean values. That is why we do not categorize these CVs as outliers.

Modifications to GR are becoming increasingly popular of late, as it is by now commonly believed that such theories might be essential to understanding the nature of gravity at cosmological scales. Apart from this, a significant unresolved issue is the inevitable singularity at the end of a gravitational collapse process as predicted by GR. EiBI theories of gravity are very attractive in this sense, as they give rise to nonsingular cosmologies as well as predict nonsingular collapse processes in the realm of classical gravity, i.e., without invoking quantum effects. Although the main effects of such modified theories of gravity are expected to set in at strong gravity scales, nonetheless, they often leave their imprint at low energy scales, in this case by a modification of the Poisson's equation in the Newtonian limit. This allows us to constrain EiBI theories via stellar structure tests, and put bounds on the parameter that determines the deviation from GR.

To this end, in this paper, we have studied EiBI theories of gravity in the context of CV binaries, and constrained the theory using available data. A total of 13 systems were chosen, with the secondary star orbiting a white dwarf primary and filling its Roche lobe, and being well described by an n = 1.5 polytropic equation of state. The Roche lobe filling condition was used to compute stellar observables numerically and comparing these with data, we have obtained a constraint on the EiBI parameter ${\bar{\kappa }}_{g}$ appearing in Equation (2) which gives, within the 5σ confidence level, $-0.315\leqslant {\bar{\kappa }}_{g}\leqslant 0.597$. It is useful to compare this with existing results on stellar bounds of the EiBI parameter in the literature, obtained by completely different methods. Casanellas et al. (2012) obtained the bound $-0.016\lt {\kappa }_{g}/({{GR}}_{\odot }^{2})\lt 0.013$ from solar physics constraints, and Avelino (2012b) obtained an upper bound ${\bar{\kappa }}_{g}\lt 4/\pi$ from the fact that the effective Jeans length in EiBI theories should be less than the solar radius. On the other hand, taking a typical white dwarf radius ∼ 10⁶ m, the results of Banerjee et al. (2017) who constrained EiBI gravity using the mass–radius relation of white dwarfs, gives $-0.239\lt {\bar{\kappa }}_{g}\lt 0.728$ at the 5σ confidence level.

In this context, we point out that the numerical treatment presented here can be improved with better modeling of the secondary stars. For example, a general velocity profile (instead of the corotational motion assumed here), as also a dynamical situation, should lead to better comparison with the observed data. Also, instead of using the volume equivalent radius, an accurate cross-section of the secondaries along the line of sight can be used to estimate the size of the secondary from the eclipse data to increase accuracy in the stellar properties. Further, we have used 13 CV systems, which somewhat limits the precision of statistical analysis. Observational data from more CV systems should constrain the EiBI parameter with higher precision.

An assumption that we have made to simplify the analysis here is that of a polytropic equation of state inside stellar matter. As we have discussed, this is an excellent approximation for low mass CV secondaries considered here, which are fully convective. Nonetheless, it might be interesting to relax this assumption and consider a model with a core-envelope structure. Recently, Chowdhury & Sarkar (2021) have explored such a model in the context of beyond-Horndeski class of models and studied how modified gravity affects stellar radius and luminosity in such models. It will be interesting to understand these issues in the context of EiBI gravity, and we expect to report on this in the near future.