© I. D. Saltas and J. Christensen-Dalsgaard 2022

1. Introduction

The need to explain the acceleration of the Universe, as well as the theoretical shortcomings of the cosmological constant paradigm have led cosmologists to question the validity of general relativity (GR) at cosmological scales. The simplest, and widely explored, scenario has been the extension of GR through a new dynamical scalar degree of freedom, which propagates a long-range gravitational force operating on cosmological scales, while it is typically suppressed at local scales to evade the stringent local gravity tests. These are known as scalar-tensor theories and go back to the pioneering work of Brans and Dicke. Their theoretical structure and phenomenology has been at the forefront of research in astrophysics and cosmology within the last decade. Past years saw the remarkable construction of general classes of scalar-tensor theories known as Beyond Horndeski (Zumalacárregui & García-Bellido 2014; Gleyzes et al. 2015) and DHOST theories (Ben Achour et al. 2016; Langlois et al. 2017), respectively. Very recently, these have been extended to the so–called U-DHOST (De Felice et al. 2022).

On the relativistic level, general families of scalar-tensor theories have an impact on the formation and dynamics of large-scale structures in the Universe (see e.g., Amendola et al. 2018; Traykova et al. 2019; Peirone et al. 2019; Hiramatsu 2022), or the structure of relativistic compact objects (see e.g., Babichev et al. 2016; Bakopoulos et al. 2022; Baake et al. 2021). Furthermore, the general class of DHOST theories has been shown to leave a characteristic imprint within Newtonian massive objects, opening up an opportunity to search for dark-energy imprints at stellar scales. It should be noted that a similar signature within Newtonian stars is expected in a broader context within geometrical extensions of GR which resemble the phenomenology of scalar-tensor theories (see e.g., Olmo et al. 2020, and references therein).

As we explain in more detail below, this local fifth-force effect has been tested in a variety of environments, from white or brown dwarf stars to neutron stars and galaxy clusters. The first work, which used the Sun to test this class of theories, was presented in Saltas & Lopes (2019), who argue that helioseismic observations have the potential to improve upon previous stellar physics constraints by two orders of magnitude. In this work, we make an important step further and study the solar equilibrium structure in the presence of the fifth force, accounting for the delicate solar microphysics, towards accurate and precise helioseismic analyses. We later show how the fifth force leaves a strong signal in the solar interior, which is detectable with helioseismic inferences of solar interior profiles.

Within the general scalar-tensor extensions of GR aiming to explain dark energy, known as DHOST theories, the weak-field limit of the theory predicts that the hydrostatic equilibrium of massive objects is modified through a new fifth-force term (Kobayashi et al. 2015; Crisostomi & Koyama 2018; Dima & Vernizzi 2018; Koyama & Sakstein 2015) as

$$\begin{array}{c}\hfill \frac{\mathrm{d}P(r)}{\mathrm{d}r}=-\rho (r)\frac{G_{\mathrm{eff}}(r)m(r)}{r^2},\end{array}$$(1)

where m(r) is the mass enclosed within a radius r from the centre and ρ(r),P(r) are the density and pressure, respectively. Within Newtonian gravity, one obviously has G_eff = G₀, with G₀ being the bare Newton’s constant measured in the Solar System. As soon as the fifth force operates, the effective Newton’s coupling on the level of hydrostatic equilibrium is promoted to a function of the radius as

$$\begin{array}{c}\hfill G_{\mathrm{eff}}\equiv (1+\frac{\mathcal{Y}}{4}{r}^2\frac{m^{″}(r)}{m(r)})G_0,\end{array}$$(2)

with ′ ≡ d/dr, and 𝒴 being the dimensionless fifth-force coupling. The coupling 𝒴 relates, parametrically, to those parameters governing the theory’s dynamics around a cosmological background. Since m″(r) < 0 within most of the star, relation (2) suggests that gravity tends to weaken (enhance) for Y > 0 (Y < 0). A plot of Eq. (2) for indicative values of 𝒴 is shown in Fig. 1. The effect of the fifth force acquires its maximum value around r ∼ 0.3 R_⊙, and it switches off at the solar surface where m″(r) = 0.

Fig. 1. Effective Newton’s coupling induced by the fifth force as a function of radius according to Eq. (1). The peak of the fifth-force effect around r ≃ 0.25 R⊙ in the radiative zone will be key for the phenomenology of solar interior profiles at hydrostatic equilibrium. The curves were evaluated according to Eq. (2) using the density profile of the reference model GS98 of Table 2, computed without the fifth force.

It is interesting to note that a similar modification of gravity’s strength is predicted in different theoretical setups (Olmo et al. 2020; Wojnar 2022). Indeed, modifications to the Newtonian law from some extended model of gravity is expected to typically involve gradients of density. Therefore, the results of our analysis can be appropriately adopted to apply to more general theory setups.

Earlier works constraining the fifth-force coupling 𝒴 were carried out in the context of white dwarfs (Jain et al. 2016; Saltas et al. 2018), red and brown dwarfs (Sakstein 2015a,b; Kozak et al. 2022), main-sequence stars (Sakstein et al. 2017), relativistic compact objects (Babichev et al. 2016), galaxy clusters (Sakstein et al. 2016; Pizzuti et al. 2019, 2022; Laudato et al. 2022; Haridasu et al. 2021), the Hulse-Taylor pulsar (Dima & Vernizzi 2018), gravitational waves (Creminelli & Vernizzi 2017; Ezquiaga & Zumalacárregui 2017; Baker et al. 2017; Sakstein & Jain 2017; Creminelli et al. 2018), and the Solar System (Crisostomi et al. 2019)¹. All constraints from stellar scales have thus far imposed 𝒴 < 𝒪(0.1) at 2σ.

A first exposition on the potential of the Sun and helioseismology to constrain this general class of dark-energy theories has been presented in Saltas & Lopes (2019). In that work, a polytropic EoS was used to predict the solar pulsation frequencies within the Cowling approximation, showing that helioseismology has the potential to tightly constraint the fifth-force coupling at the 10⁻³ level. Our present goal is to provide an accurate, quantitative description of how the fifth force affects the solar interior towards precision constraints on the fifth force with helioseismology. The work splits into two parts. In the first part (Sect. 2), we study the effect of the fifth force on solar observables by means of intuitive arguments and analytic relations. In the second part (Sects. 3 and 4), we numerically investigate the evolution of the present Sun with numerical simulations, and analyse the interplay between the fifth force and solar input physics. The focus is on the effect on the solar interior profile, most notably the sound speed, which is the key to the precise science of helioseismology. We later show that the fifth-force effect gives a characteristic effect deep in the interior zone, at a point where uncertainties from opacity or diffusion become negligible, thus making it distinct within helioseismic analyses. In Sect. 5 we derive an order of magnitude constraint for the fifth-force coupling based on its effect on the solar sound speed, and in Sect. 6 we draw our conclusions.

2. Fifth force in the solar interior: Theoretical considerations

The dynamics of stellar evolution are modelled by the equations describing hydrostatic equilibrium, energy conservation and transport, coupled to the time evolution of element abundances in the star. The system of equations closes under an assumption on the equation of state (EoS) and other microphysics, in particular, the opacity, diffusion, and the treatment of convective energy transport. Although assumptions on microphysics do not affect gravity directly, they may implicitly interfere with it through their impact on the stellar interior profiles such as density and pressure. Let us recap the basic equations and discuss the way they are affected by the presence of a fifth force.

Overall, a weakening (strengthening) of the gravitational constant G in the solar interior will lead to a change in pressure, as can be easily seen using the equation of hydrostatic equilibrium to estimate $P_{\text{central}}\sim G{M}_{\odot }^2/R_{\odot }^4$. The fractional change in pressure, δP/P, is accompanied by a fractional change in the temperature, δT/T, of similar magnitude. Obviously, a weakening (strengthening) of gravity will tend to make the star less (more) dense. In order to retain hydrostatic equilibrium at the same radius, the star will inevitable adjust its structure, which affects the energy generation from nuclear reactions, as we discuss later.

The energy equation relates the spatial gradient of luminosity to the energy generation, and will be not explicitly affected by the change in gravity. The same holds for the equation describing the change of element abundances with time. Both equations, however, will be affected implicitly through the change in the core’s temperature, density and initial element abundance.

The equation describing the evolution of element abundances with time will not receive explicit corrections, but will be implicitly affected through the change in the energy generation rate, which for the pp chain is ${ϵ}_{\text{pp}}\propto X_{\text{H}}^2\rho T^6$. The change in the central temperature, density and hydrogen abundance will induce a change in the energy generation. The change in the initial hydrogen abundance (X_H) follows from the need to maintain the same luminosity at the present solar age, as we discuss later on.

2.1. Density profile

Let us start by understanding the behaviour of the solar interior profiles from a quantitative point of view, through the derivation of approximate, and whenever possible, analytical solutions. An often useful EoS to model a star in a simple manner is the polytropic relation defined through

$$\begin{array}{c}\hfill P=K{\rho }^{{\mathrm{\Gamma }}_1}\equiv K{\rho }^{1+\frac{1}{n}},\end{array}$$(3)

where K is the adiabatic constant, Γ₁, n the (constant) adiabatic and polytropic index respectively, and P, ρ the pressure and density profiles. The solar radiative zone is crudely described by n = 3 and here the polytropic equation needs to be solved numerically. In the convective zone, Eq. (3) is approximately satisfied, with Γ₁ being the adiabatic compressibility, Γ₁ = (∂log p/∂log ρ)_ad; in the bulk of the Sun matter is approximately in the form of a fully ionised ideal gas, with Γ₁ ≃ 5/3, corresponding to n = 1.5. In this case an approximate analytic solution can be found, as discussed below.

Polytropic solutions for n = 3 in a solar context were previously discussed in Saltas & Lopes (2019). As was shown there, the effect of the fifth force on the sound-speed profile in the radiative zone grows from the centre until it peaks around r ≃ 0.3 R_⊙, similar to the scaling exhibited by the effective gravity strength in Fig. 1. A similar behaviour is expected for the fractional change of the density profile under the fifth force, and confirmed by the evolutionary simulations we performed in this work, to be presented later on. To understand the aforementioned behaviour of the sound speed in the radiative zone quantitatively, we model the radiative interior as an ideal gas, P ∝ ρk_BT, such that

$$\begin{array}{c}\hfill c_{\mathrm{s}}^2(r)\equiv {\mathrm{\Gamma }}_1\frac{P}{\rho }\simeq \frac{{\mathrm{\Gamma }}_1{k}_{\mathrm{B}}T}{\mu m_{\mathrm{u}}},\end{array}$$(4)

where k_B is the Boltzmann constant, μ the mean molecular weight, and m_u the atmomic mass unit. This approximation will be useful later on when we interpret the results of our numerical simulations. The adiabatic index Γ₁ intimately relates to the properties of the EoS, and the intuition suggests that the effect of the fifth force on it will be minor, and only indirect through the change in the density and pressure. Therefore, the main effect on $\delta c_{\mathrm{s}}^2/c_{\mathrm{s}}^2$ in the radiative zone will be driven by the change in the temperature and mean molecular weight,

$$\begin{array}{c}\hfill \frac{\delta c_{\mathrm{s}}^2(r)}{c_{\mathrm{s}}^2}\simeq \frac{\delta {\mathrm{\Gamma }}_1}{{\mathrm{\Gamma }}_1}+\frac{\delta T}{T}-\frac{\delta \mu }{\mu }.\end{array}$$(5)

We discuss this in detail in Sect. 4.

Let us now turn our focus on the derivation of analytic solutions for solar profiles in the solar convective zone. The mass enclosed in a spherical shell up to a radius around the base of the convective zone is fairly close to the total solar mass, so we may write for the mass profile in this region

$$\begin{array}{c}\hfill m(r)\simeq M_{\odot }+\mathcal{O}\left(\frac{\delta m}{M_{\odot }}\right),\delta m=M_{\odot }-m(r>R_{\mathrm{cz}}),\end{array}$$(6)

with the neglected terms sufficiently small, given that the convective zone accounts for 2.5% of the solar mass. We first examine the solution for 𝒴 = 0. Assuming m(r)≃M_⊙, we approximate the r.h.s of the equation of hydrostatic equilibrium as

$$\begin{array}{c}\hfill \frac{\mathrm{d}P}{\mathrm{d}r}=K{\mathrm{\Gamma }}_1{\rho }^{{\mathrm{\Gamma }}_1-1}\frac{\mathrm{d}\rho }{\mathrm{d}r}=-\frac{G_0{M}_{\odot }}{r^2}\rho [1-\mathcal{O}\left(\frac{\delta m}{M_{\odot }}\right)],\end{array}$$(7)

with the second equality due to the polytypic EoS. Assuming Γ₁ to be constant, as is approximately true except in the near-surface ionisation zones of hydrogen and helium, and neglecting the small correction on the r.h.s, we can straightforwardly integrate the above equation to find

$$\begin{array}{c}\hfill {\rho }_{(0)}(r)\equiv B{(\frac{R_{\odot }}{r}-1)}^{\frac{1}{{\mathrm{\Gamma }}_1-1}}B\equiv {\left(\frac{{\mathrm{\Gamma }}_1-1}{{\mathrm{\Gamma }}_1}\frac{G{M}_{\odot }}{K{R}_{\odot }}\right)}^{\frac{1}{{\mathrm{\Gamma }}_1-1}},\end{array}$$(8)

with B a constant with dimension of density. This solution describes approximately the density profile in the solar convective zone. Notice that in principle this solution needs be matched with conditions in the radiative zone at the convective boundary.

We now consider the extension of the previous analytic expression to the case where the fifth force is turned on, 𝒴 ≠ 0. Since the fifth-force coupling is small, we adopt a perturbative approach. In particular, we start from the general hydrostatic Eq. (7), and assume a solution for the density profile in the form

$$\begin{array}{c}\hfill \rho (r;\mathcal{Y})\simeq {\rho }_{(0)}(r)+ϵ{\rho }_{(\mathcal{Y})}(r),\end{array}$$(9)

where ϵ ≡ π𝒴 is a small parameter ≪1, and ρ_(𝒴)(r) the small contribution on the density profile from the fifth force. A similar perturbative solution will hold for the sound speed, which we derive later. To ease the computation, instead of deriving the solution for any Γ₁, we now restrict ourselves to Γ₁ = 5/3, which corresponds to the choice of a polytropic index n = 3/2. Plugging Eq. (9) back into the hydrostatic Eq. (7), the ϵ-dependent part of the equation takes the form

$$\begin{array}{c}\hfill z(1-z){\rho }_{(\mathcal{Y})}^{\prime }-{\rho }_{(\mathcal{Y})}+P(z)=0,z\equiv \frac{r}{R_{\odot }},\end{array}$$(10)

with ′ ≡ d/dz, and P(z)∝1/2 − 3z + (9/2)z² − 2z³, modulo an overall normalisation constant. The homogeneous solution is easily found first setting P(z) = 0, and then integrating through separation of variables. This yields

$$\begin{array}{c}\hfill {\rho }_{(\mathcal{Y})}^{\mathrm{homog}.}=C\frac{z}{1-z},\end{array}$$(11)

with C a constant of integration. Since the full solution will be given by the sum of the homogeneous and a inhomogeneous solution, to find the latter, we seek a solution of Eq. (10) in the form ${\rho }_{(\mathcal{Y})}^{\mathrm{inhomog}.}={\rho }_{(\mathcal{Y})}^{\mathrm{homog}.}Q(z)$, where Q a function to be determined. Substituting this ansatz in Eq. (10), one arrives at the differential equation Q′= − P(z)/z², which can be easily integrated to find Q, since the explicit form of the polynomial P(z) is known. After choosing the integration constant C in the homogeneous solution so that the solution for Eq. (10) is well-behaved on the solar surface, that is no divergence occurs as z → 1, one finally finds

$$\begin{array}{c}\hfill ϵ·{\rho }_{(\mathcal{Y})}(r)\simeq c(-1-7z+2z^2+\frac{6zlog(z)}{z-1}),\end{array}$$(12)

with the constant c defined as

$$\begin{array}{c}\hfill c\equiv \mathcal{Y}\left(\frac{3\pi }{10}\frac{B^{4/3}G{R}_{\odot }^2}{{\mathrm{\Gamma }}_{1}K}\right),\end{array}$$

and we also re-introduced our expansion parameter ϵ ≡ π𝒴, appearing in Eq. (9). The solution (12) is understood as a small correction on top of the 𝒴 = 0 solution, ρ₍₀₎, and it can be seen that it vanishes on the solar surface, that is to say in the limit of z → 1. One further notices the different scaling with radius z ≡ r/R_⊙ compared to the one derived in standard gravity (=3/2), Eq. (8). We now define the fractional change of density as

$$\begin{array}{c}\hfill \frac{\delta \rho }{\rho }\equiv \frac{\rho (\text{standard model})-\rho (\text{modified model})}{\rho (\text{standard model})},\end{array}$$(13)

with ρ = ρ(r) understood as in Eq. (9). The definition (13) will serve as a definition for the fractional change of interior profiles under any change of physics. It is straightforward to check that δρ/ρ, like ρ_(𝒴), also vanishes on the solar surface. To understand the qualitative behaviour of the fifth-force term on the fractional change of density, we may expand to first order in z⁻¹ ≡ R_⊙/r around the base of the convective zone which we conventionally choose as R_cz ≃ 0.71 R_⊙,

$$\begin{array}{c}\hfill \frac{\delta \rho }{\rho }(r\gtrsim R_{\text{cz}})\simeq -0.57\mathcal{Y}-0.74\mathcal{Y}\left(\frac{r-R_{\mathrm{cz}}}{R_{\odot }}\right).\end{array}$$(14)

For the typical value of 𝒴 = 0.01, the above solution implies a very slow change in the convective zone, which is confirmed by our detailed simulations, to be discussed in Sect. 3.2.

2.2. Sound speed

We now look at the sound-speed profile defined through

$$\begin{array}{c}\hfill c_{\mathrm{s}}^2(r)\equiv {\mathrm{\Gamma }}_1\frac{P}{\rho }\approx K{\mathrm{\Gamma }}_1{\rho }^{{\mathrm{\Gamma }}_1-1},\end{array}$$(15)

with the last approximate relation is valid for a polytrope. Using the analytic solutions for the density derived earlier, it is easy to derive an expression for the sound speed in the convective zone and at zero fifth force (𝒴 = 0) as

$$\begin{array}{c}\hfill c_{s(0)}^2(r)\simeq \frac{G_0{M}_{\odot }}{R_{\odot }}({\mathrm{\Gamma }}_1-1)(\frac{R_{\odot }}{r}-1),\end{array}$$(16)

reminding us that we neglected corrections due to the non-constancy of the mass profile in the solar convective zone. It is also understood that this relation needs be matched with the polytropic one in the radiative zone, which will fix the polytropic constant K – we omit doing this here, since we are only interested in understanding the qualitative scaling. The solution (16) was first derived in Christensen-Dalsgaard et al. (1991).

If we switch on the fifth force, the relevant sound-speed profile is then straightforwardly found to be

$$\begin{array}{c}\hfill c_{\mathrm{s}}^2(r;\mathcal{Y})\simeq K{\mathrm{\Gamma }}_1{({\rho }_{(0)}(r)+ϵ{\rho }_{(\mathcal{Y})}(r))}^{{\mathrm{\Gamma }}_1-1},\end{array}$$(17)

which applies, under our assumed approximations, in the convection zone. Thus, here the speed of sound will inherit the density’s approximate constancy, which is also seen in the full numerical solutions of calibrated solar models, to be presented later on. A key observation from Eq. (16) is that the sound-speed profile is independent of the adiabatic constant K, making it rather insensitive to change of input physics such as opacity and diffusion, in the convective zone. Similar to what we did before, we expand to first order around the base of the convective zone (≃0.71 R_⊙), assuming again Γ₁ = 5/3, to find

$$\begin{array}{c}\hfill \frac{\delta c_{\text{s}}^2}{c_{\text{s}}^2}(r\gtrsim R_{\text{cz}})\simeq 0.1\mathcal{Y}-0.5\mathcal{Y}\left(\frac{r-R_{\text{cz}}}{R_{\odot }}\right).\end{array}$$(18)

This explains the mild variation of $\delta c_{\text{s}}^2(r)/c_{\text{s}}^2$ for r ≳ R_cz. The gradient of the sound speed forms a critical diagnostic for solar structure in the convective envelope, and it is usually defined through the quantity W(r) (e.g., Gough 1984; Däppen & Gough 1986) as

$$\begin{array}{c}\hfill W\equiv \frac{r^2}{G_{0}m(r)}\frac{\mathrm{d}{c}_{\text{s}}^2}{\mathrm{d}r}\simeq \frac{r^2}{G_0{M}_{\odot }}\frac{\mathrm{d}{c}_{\text{s}}^2}{\mathrm{d}r},\end{array}$$(19)

with the last approximate equality true within the convective zone. In the convective zone, W ≃ 1 − Γ₁, with Γ₁ ≃ 5/3 and therefore W ≃ −2/3. This in turn can be used to infer the thickness of the convective envelope if the sound-speed profile is known through helioseismic observations (Christensen-Dalsgaard et al. 1991). The applicability of this procedure in the presence of the fifth force relates to whether Γ₁ receives sizable corrections from the fifth force, such that the value of W would depart significantly from −2/3. Our simulations, to be presented below, show that under the fifth force δΓ₁/Γ₁ ∼ 𝒪(10⁻⁵), and therefore the effect is negligible. In Fig. 2 we show W(r) for calibrated solar models within standard and modified gravity. It is seen that a stronger gravity in the convective zone (𝒴 < 0) will act as to shift the bottom of the convective zone to a lower point compared to standard gravity, and the opposite is true for 𝒴 > 0. In other words, a stronger gravity implies a steeper sound-speed profile.

Fig. 2. The quantity W(r), defined in Eq. (19). This is a measure of the sound-speed gradient and a key diagnostic for the solar structure. Here, it is shown for the following three cases: for standard gravity (𝒴 = 0), weaker (𝒴 = 0.01) and stronger (𝒴 = −0.01) gravity respectively. The plots are based on full numerical simulations of calibrated models for radius, luminosity and surface metallicity (see Table 2 for a summary). Top: the solution of W(r) = − 2/3 yields an approximate estimate of the position of the radius of the base of the convective zone (r ≃ Rcz), which is shown in the upper panel. The reference model at standard gravity is model GS98 summarised in Table 2. Bottom: the peak of the quantity W(r) right below the solar surface due to the ionisation of helium. Although barely visible by eye, the position of the peak is shifted inwards (outwards) for stronger (weaker) gravity.

3. Solar calibration method in the presence of fifth force

3.1. Numerical method

We use the one-dimensional MESA stellar evolution code (Paxton et al. 2011, 2013, 2015)². The code solves the system of stellar equations discretely at each time step on a one-dimensional spatial grid, assuming quasi-hydrostatic equilibrium and spherical symmetry. Solutions to the equations are then sought using a variant of the Newton-Raphson method known as the Henyey method (Henyey et al. 1964). The equations of stellar structure form a boundary value problem with boundary conditions at the centre and surface of the star (e.g., Kippenhahn et al. 2012). The former are defined by the requirement that m = 0, L = 0 at r = 0, while for the latter, an atmospheric model is required. For the latter, we use a grey atmospheric model, based on a T − τ relation described by the Hopf function modelled on Vernazza et al. (1981). For the hydrostatic equilibrium in the atmosphere we assume standard gravity, since the effect of the fifth force in that region is irrelevant.

The fifth-force term in Eq. (1) depends on the second derivative of the mass profile with respect to radius, which we can express as

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^{2}m(r)}{\mathrm{d}{r}^2}=\frac{\mathrm{d}m(r)}{\mathrm{d}r}(\frac{\mathrm{d}log\rho }{\mathrm{d}r}+\frac{2}{r}),\end{array}$$(20)

which is the form which we also adopt for the numerical implementation. Here we note that in the above equation, dm/dr = 4πρ(r)r². The derivatives in Eq. (20) are implemented in a discreet manner as

$$\begin{array}{c}\hfill \frac{\mathrm{d}\rho (m)}{\mathrm{d}m}\simeq \frac{{\rho }_{k+1}-{\rho }_k}{m_{k+1}-m_k},\end{array}$$(21)

where k labels the kth cell on the grid and all physical variables are evaluated either at the centre or the edge of the kth cell. The fifth force is then evaluated at each time step (t_k) next to the standard Newtonian term, from the outermost (k = 1) to the innermost cell corresponding to the stellar surface and core respectively. For a typical run, we use about ∼1000 spatial grid points and a time step of ∼10⁷ yr.

Our microphysics assumptions are as follows. For our reference EoS we use a blending of OPAL (Rogers & Nayfonov 2002) and Free EoS³, whereas the opacities are based on the tabulations from OP (Seaton 2005) with low- and high-temperature extensions from Ferguson et al. (2005) and Buchler & Yueh (1976) respectively. Nuclear and weak reaction rates are from JINA REACLIB (Cyburt et al. 2010) and Fuller et al. (1985), Oda et al. (1994), Langanke & Martínez-Pinedo (2000) respectively. We also choose to set overshooting to zero, an assumption which does not restrict the generality of results. In particular, some of our key results presented here we confirmed allowing for a small (exponential) overshooting. An important input is the initial fraction of metals, which affects the computation of opacities too. Unless otherwise stated, we choose to work with the so–called GS98 fraction based on Grevesse & Sauval (1998), due to their good agreement with helioseismic data. Later on, we discuss the significance of this choice in the light of the revised solar metallicity of Asplund et al. (2009, A09) and the abundance problem.

3.2. Solar evolution simulations

It is important to first lay out our framework for testing gravity with the Sun. Solar observables such as luminosity and surface temperature are known to very high precision, whereas thanks to helioseismic observations the same holds for the interior solar profiles such as density and speed of sound. Our approach to test the fifth force is as follows: We choose to work with calibrated solar models, that is models which match the equilibrium observables of the Sun within their 1σ errors. Any effect of the fifth force will then impact the structure of the interior solar profiles, and the latter will be our key quantities for investigation.

The computation of calibrated solar models corresponds to finding that subspace of initial conditions which reproduce the observables of the present Sun. Our main present-day observables are the total stellar radius (R), surface element abundances ((Z/X)_surf.) and luminosity (L; see also Table 1). The latter can be traded for the effective temperature T_eff given the radius R. In addition, we may include results of helioseismic analyses, such as the radius R_cz at the base of the convection zone, the envelope helium abundance, or the inferred solar sound speed. Our initial parameters correspond to the mixing-length parameter (α_MLT), initial helium (Y₀) abundance and initial metallicity (Z₀). On top of the previous parameters, we need to define the value of the fifth-force coupling 𝒴. Therefore, our parameter space is four-dimensional. We choose to fix our initial conditions at the stage of the zero-age main sequence, from which point the star is left to evolve until today. A calibrated solar model corresponds to a model which simultaneously lies within the statistical tolerance of the observables.

We start with a solar calibration without the fifth force, that is 𝒴 = 0. This enables us to pin down the three-dimensional parameter space of solar parameters before we turn on the fifth force, and will produce our base model at standard gravity. We employ the simplex method, which successively seeks for the values of the initial parameters (p_i) minimising a χ² functional built out of the desired observables (O_i),

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _i}\frac{{[O_i(p_j)-O_i(\text{Sun})]}^2}{{\sigma }_i^2}.\end{array}$$(22)

Typically, for a solar calibration,

$$\begin{array}{c}\hfill p_i=\{{\alpha }_{\mathrm{MLT}},Y_0,Z_0\},O_i=\{T_{\odot }^{\mathrm{eff}},R_{\odot },{(Z/X)}_{\mathrm{surf}}\}.\end{array}$$

Though this calibration approach is conceptually different from the traditional method of linearisation, the two methods should coincide modulo small differences of numerical origin. We have verified this agreement between the two procedures for indicative solar calibrations. In our computations we employ the GS98 metallicity mixture, but for the sake of comparison we also calibrate a model based on the revised abundance A09 (Asplund et al. 2009; for a discussion on the solar abundance issue see Sect. 4.3). The set of representative calibrated solar models we use in our analysis are summarised in Table 2.

We proceed turning on the fifth-force term using our GS98 calibration as our reference model, which extends our parameter space to

$$\begin{array}{c}\hfill p_i=\{{\alpha }_{\mathrm{MLT}},Y_0,Z_0,\mathcal{Y}\}.\end{array}$$

Given that we work with calibrated models, the effect of the fifth force will manifest itself in the interior profiles. Some useful intuition on the effect of the fifth force on a calibrated model can be drawn from homology arguments. In this context, the luminosity is related to the total mass, radius and mean molecular weight through (Kippenhahn et al. 2012),

$$\begin{array}{c}\hfill L\propto \frac{G^{15/2}{\mu }^{15/2}{M}^{11/2}}{R^{1/2}},\end{array}$$(23)

assuming Kramers opacity and μ = 4/(3 + 5X − Z) for the mean molecular weight. At fixed mass and radius, and for 𝒴 > 0, if the luminosity L is to remain the same, the strengthening of gravity in the solar core (see Fig. 1) requires the compensating effect of a decreased μ. Figure 3 shows the effect of the fifth force on some solar observables, for uncalibrated models, while Fig. 4 shows the predicted sound-speed profile against the solar one for different metallicity mixtures. From Fig. 3 one sees that the fifth force predominantly affects the total radius and the thickness of the convective envelope, and to a smaller extent temperature.

Fig. 3. Variation of total radius, base of convective zone and effective radius under the fifth-force for uncalibrated solar models. Each observable is normalised with respect to its value at 𝒴 = 0 (no fifth force) at a fixed initial helium abundance (Y0), initial metallicity (Z0) and mixing-length parameter (αMLT). Whereas the fifth force has a significant effect on the radius and position of the convective zone, its effect is milder on the effective temperature.

Fig. 4. Fractional difference for the sound speed in the sense (Sun) – (Model). ‘Sun’ refers to the sound-speed profile as inferred from helioseismic observations at standard gravity according to the results of Rhodes (1997). Models GS98 and A09 correspond to theoretically computed calibrated solar models at standard gravity (𝒴 = 0) with the two distinct, and most popular metallicity mixtures for the Sun (see also Table 2 and Sect. 3). For comparison, we also show two computed models with the fifth force turned on (𝒴 = ±0.01) based on the GS98 mixture. We note that the fifth-force effect peaks in the radiative zone which will be key for constraining the fifth force with helioseismic observations.

4. Input physics systematics and fifth force

Searching for a departure from Newtonian gravity in the Sun is a rather delicate task, due to interplay between the well-observed solar structure and the small fifth-force effect. Therefore, an understanding of the systematics that could hinder a clean extraction of the fifth-force effect is necessary. In this section, we use our numerical results to attack this issue. We later show that the effect of the fifth force leaves a characteristic and clean signature on the solar sound-speed profile which allows it to be accurately distinguished from typical theoretical uncertainties in the modelling of the Sun.

4.1. Opacity and diffusion

Realistic opacity profiles for the Sun have a highly non-trivial dependence on the density, temperature and element abundances, and they contribute to one of the most important theoretical uncertainties in solar modelling. For example, an increase in opacity at standard gravity, will tend to increase the density in most parts of the star, an effect that can be also achieved by strengthening gravity at constant opacity. The uncertainty of opacity modelling in the Sun ranges from ∼2% near the centre of the Sun, to an increase of up to 7%±4% around the base of the solar convective zone, where element ionisation becomes significant. The 2% uncertainty close to the centre is typically computed as the averaged opacity difference between the two popular opacity tabulations, namely the OPAL (Iglesias & Rogers 1993, 1996) and OP opacity tables. The 7% uncertainty results from a comparison of opacity tables with a measurement of iron’s opacity at the laboratory. Here, as a conservative estimate, we adopt a constant error of in the opacity of

$$\begin{array}{c}\hfill \kappa ={\kappa }_0\pm 5\%,\end{array}$$(24)

with κ₀ the reference opacity profile based on the OP tabulation, and the chemical mixture GS98 (see also Table 2). We implement this uncertainty by increasing (decreasing) the opacity throughout the star by the same constant factor, 1 ± 0.05.

The diffusion error is defined as the theoretical uncertainty in the diffusion coefficients for the various elements. The latter govern the efficiency of diffusion in the star, and in the context of the diffusion treatment of Thoul et al. (1994) employed in our simulations, the comparison between different solutions for the diffusion equations can be used to place a conservative uncertainty of 15%. We therefore use for the uncertainty on diffusion,

$$\begin{array}{c}\hfill D_i=D_i\pm 15\%,\end{array}$$(25)

where D_i the diffusion coefficients related to element i. Notice that the uncertainty is applied the same to the coefficients of all elements. For a further discussion on the topic we refer to Vinyoles et al. (2017) and references therein.

The effects of variations in opacity and diffusion (at standard gravity), as well as that of the fifth force (at standard opacity and diffusion), are shown in the model differences of Fig. 5. Here the effect on the thermodynamic properties of the change in composition is illustrated in terms of the mean molecular weight, approximately calculated as

$$\begin{array}{c}\hfill \mu =\frac{4}{3+5X-Z}.\end{array}$$(26)

Fig. 5. Fractional model differences in the sense (reference model) – (modified model) for solar interior profiles. This illustrates the importance of the most significant uncertainties in solar modelling (opacity, diffusion) against the fifth-force effect. The legend for the curves is shown in the top left panel. The dot-dashed (dotted) curves correspond to an increase (decrease) of the respective input physics parameter. For the combined change of opacity and diffusion the dot-dashed curves correspond to increased opacity and decreased diffusion. The reference model is the calibrated solar model at standard gravity with GS98 metallicity and with unmodified opacity and diffusion coefficients, as explained in the text and Table 2. In all panels, the models with the fifth force turned on (black, solid curves) are computed on top of our reference model at standard gravity, and have a value of the fifth-force coupling |𝒴| = 10−2 for illustration purposes. As conservative estimates for the uncertainties in input physics we use ±5% for opacity and ±15% for diffusion, which are computed on top of our reference model (see Sect. 4). We also consider the combined effect of the simultaneous uncertainty of opacity and diffusion. The mild variation of the density and sound-speed profiles around the base of the convective zone (≃0.71 R⊙) has been explained through analytic relations in Sect. 2. (a) Density. (b) Speed of sound. (c) Temperature. (d) Adiabatic index. (e) Mean molecular weight. (f) Pressure.

In addition to the individual effects of opacity and diffusion, we also show the combined effect resulting from the simultaneous variation of both. At standard gravity, an increase in opacity will generally tend to increase the density, temperature and mean molecular weight, and decrease the sound speed. A similar situation holds for an increase in the efficiency of diffusion, but in that case, the sound-speed profile decreases. Now, for a fifth-force coupling 𝒴 > 0, density and pressure will generally tend to decrease for most of the solar region, while the sound speed will decrease. The opposite situation holds for 𝒴 < 0. The increase of density, pressure and temperature before the turn over at r ≃ 0.2 R_⊙ in Fig. 5, is due to the increase of gravity in that region, as follows from the scaling in Eq. (1).

One also notices that the change in sound speed in the convection zone resulting from changes in opacity or diffusion is small compared to the changes in density or pressure. Within the polytropic solutions, Eqs. (8) and (16), the structure generally depends on the constant K, which is different for each modified model. This directly affects density (see Eq. (8)) and hence pressure, while Eq. (16) shows that the sound speed is independent of K. In contrast, Eq. (18) shows that the fifth force directly affects the sound speed in the convection zone, leading to the behaviour seen in Fig. 5.

The peak around r = 0.25 R_⊙ in the fifth-force effect on density, pressure and sound speed shown in Fig. 5 was first observed in Saltas & Lopes (2019) under the crude approximation of an n = 3 polytrope. However, the striking feature shown here is that it not only persists within our accurate simulations, but it occurs at a point where the systematics from opacity and diffusion become negligible compared to it. To understand this better, we can use the ideal gas approximation for the EoS (see Eq. (5)). Then,

$$\begin{array}{c}\hfill \frac{\delta c_{\mathrm{s}}}{c_{\mathrm{s}}}\simeq \frac{1}{2}(\frac{\delta T}{T}-\frac{\delta \mu }{\mu }),\end{array}$$(27)

where we neglected the insignificant variation of Γ₁ under model differences compared to the other terms (see also the model differences of Fig. 5). Under variations of opacity and/or diffusion at standard gravity, the contributions coming from the variation of temperature and mean molecular weight approximately cancel with each other around the point r ≃ 0.25 R_⊙. This is not anymore true for model variations under the fifth force (at unmodified opacity and diffusion), where the temperature traces the change of the effective strength of gravity with radius in that region (see also Eq. (1)). An illustration of this mechanism is shown in Fig. 6. This surprisingly clean signal of the fifth force allows for precision tests of the theory with helioseismic observations.

Fig. 6. Top: model differences for the sound-speed profile similar to those of Fig. 5, in the sense (reference model) – (modified model), but with a smaller fifth-force coupling 𝒴 = 10−3. The plot zooms into the solar region where our estimates of systematics from opacity and diffusion may become negligible compared to that of the fifth force, in order to clearly illustrate the peak of the fifth-force effect in the radiative zone. The error bar shows approximately the error of helioseismic inversions at that point, ∼10−4, to highlight the precision down to which helioseismic reconstructions of the solar sound speed can constrain the coupling 𝒴. The insignificance of the opacity and diffusion uncertainties compared to the fifth force for this solar region is discussed in Sect. 4. Bottom: an illustration of the mechanism behind the ‘clean’ peak of the fifth-force effect on the speed of sound profile around r ≃ 0.25 R⊙, where the systematics from uncertainties related to opacity and diffusion become negligible compared to the effect of the fifth force (see also Fig. 5). The reference model at standard gravity is the GS98 model (see Table 2 and text). The plots assume for illustration the ideal-gas approximation for the sound-speed profile, where $c_{\text{s}}^2\propto {\Gamma }_{1}T/\mathrm{\mu }$, and neglect the small variation of the adiabatic index Γ1 (see also Fig. 5). Right panel: fractional difference of the sound-speed profile and its components under an increase of the opacity by 5%, but the picture is similar for uncertainties in diffusion. Left panel: how the picture changes for the fifth-force effect. For a detailed discussion we refer to the text around Eq. (27).

Before we close this section, we want to challenge our assumption of a constant uncertainty in the opacity throughout the star. In particular, we ask what would be the effect of a radius-dependent uncertainty which would peak around the region where the fifth-force effect peaks. We approach this by looking at the difference between two popular opacity tabulations, namely OPAL and OP, which provides an estimate of the opacity error. The model differences between the latter opacity tabulations have been previously studied in the literature (at standard gravity), see Fig. 30 of Christensen-Dalsgaard (2021), for example, and it is also illustrated in Fig. 7. The computation of Christensen-Dalsgaard (2021) suggests that the difference in the sound-speed profiles computed with OPAL and OP, respectively, shows a peak at r ≃ 0.3 R_⊙, inducing a reduction of the sound speed which can reach up to 0.2%. Such a localised opacity uncertainty could interfere with the fifth-force effect and we comment on it further in Sect. 5 in light of new constraints on 𝒴. A detailed confrontation of different opacity tabulations against helioseismic data has been performed in Villante et al. (2014) where it was found that standard solar models computed with the OP opacity do not provide a good fit to the sound-speed profile as inferred from helioseismology in the region 0.3 < r/R_⊙ < 0.6. We stress though, that such analyses would have to be in principle revisited in the presence of the fifth force, since one now has one extra degree of freedom to be accounted for. Furthermore, opacity-related uncertainties are also related to uncertainties in the metallicity fraction. As illustrated in Fig. 4, the model prediction with the proposed metallicity mixtures A09 shows tension with helioseismology, which is known as the solar abundance problem. If the metallicity measurements of A09 are assumed to be correct, the discrepancy is typically attributed to a mismatch in the modelling of the opacities in the solar interior. We return to this in Sect. 4.3, but note already here that the very recent analysis of Magg et al. (2022) suggests that the tension is significantly reduced through their revised mixture, close to the GS98 one.

Fig. 7. Difference between solar profiles computed with different opacity tabulations providing an estimate of the error coming from opacity modelling. Here, we show the difference in the sound speed for a model computed with our reference opacity table (OP) and a model computed with an alternative opacity (OPAL). Both computations assume standard gravity (𝒴 = 0), and are based on the GS98 metallicity mixture. In our analysis of systematics, the uncertainty in the modelling of opacity was considered constant throughout the star, as explained in Sect. 4. However, one cannot exclude the possibility of a localised opacity variation such as the one emerging from the comparison between OP and OPAL opacity tabulations. For comparison purposes we also show the corresponding impact of the fifth force on the sound speed computed with our reference opacity (OP) and for two values of the fifth force coupling. We discuss the impact of this on inferences of the fifth-force coupling 𝒴 in Sect. 5.

4.2. Equation of state

The EoS affects the interior profiles, including the sound speed, in two different ways. One is the direct effect of different input physics on the adiabatic index Γ₁, as the latter provides a proxy on the way the pressure relates to density with radius. The effect of the fifth force on the adiabatic index is shown in Fig. 5. It is evident that fifth-force effect on Γ₁ is negligible throughout the Sun, reaching a maximum of about 0.002% in the convective zone, where it is comparable to the respective effect of diffusion. It should be also noted that for most of the radiative zone its effect is smaller than the effect of opacity and diffusion. However, the EoS also determines the relation between pressure, density and temperature and hence the overall structure of the model.

To test these effects, we produce two solar models at standard gravity and GS98 metallicity, one based on pure Free EoS, and another on pure OPAL EoS respectively. We remind that our reference EoS is based upon the blended OPAL/Free EoS. Figure 8 shows a comparison between the choice of EoS at standard gravity and that of the fifth force at fixed reference EoS for the sound-speed profile. For 𝒴 = 10⁻² the fifth-force effect on the sound speed tends to be larger than that of the EoS in the radiative zone by at least one order of magnitude. However, for 𝒴 = 10⁻⁴ the two effects become of similar order of magnitude. One also notices that for the solar region in the radiative zone where the fifth-force effect peaks, both alternative EoS choices, that is pure Free EoS and pure OPAL EoS respectively, yield a variation of the sound speed of the same order. Clearly, the choice of the EoS for sufficiently small values of 𝒴 matters. Furthermore, the good agreement between our reference EoS and the alternative EoS choices in the radiative region where the fifth-force effect peaks (r ≃ 0.25 R_⊙) suggests that our results are not expected to be biased by the choice of EoS. Finally, we should note that relativistic effects of electrons close to the solar centre could potentially become comparable to the fifth-force effect. However, the EoS’s we use already implement such effects, and therefore there is no concern about their interference with the fifth force.

Fig. 8. Comparison between the effect of a modification of the EoS at standard gravity, and that of the fifth force at fixed EoS on the sound-speed profile. The variation δcs is in the sense (reference model) – (modified model), and the vertical axis is on natural logarithmic scale. Our reference model uses the mixed OPAL/Free EoS, whereas the alternative EoS’s at standard gravity are based on pure Free EoS (‘modified EoS 1’), and pure OPAL (‘modified EoS 2’) respectively. For more details we refer to Sect. 4. The models with fifth force in the figure are computed on top of our reference EoS. We remind that the fifth-force effect on the sound-speed profile within the radiative zone peaks around r ≃ 0.25 R⊙ (see also Figs. 5 and 6).

4.3. The solar metallicity problem in light of the fifth-force effect

Our analysis has been based on the GS98 choice for the solar metals mixture. Here, we ask how the solar metallicity problem might affect our results. The tension between solar models and helioseismic observations for models computed with the revised abundances, Asplund et al. (2009, A09), compared to those computed with the ‘old’ solar metallicity inference of Grevesse & Sauval (1998, GS98) is shown in Fig. 4. As seen there, the predicted sound-speed profile, based on the metal mixture A09, is in tension with the profile predicted by helioseismic data, compared to calibrated models based on the GS98 metallicity mixture. If one assumes that the revised metallicity mixture (A09) is correct, the problem has been commonly attributed to mismatch in opacity modelling based on the GS98 mixture.

According to Christensen-Dalsgaard et al. (2009), Christensen-Dalsgaard & Houdek (2010) the mismatch between the A09 and GS98 models can be bridged⁴ if we adjust the opacity profile of the A09 model κ through a radius-dependent opacity correction Δκ according to

$$\begin{array}{c}\hfill log{\kappa }^{\mathrm{new}}(\rho ,T,X_i)=log\kappa (\rho ,T,X_i)+\mathrm{\Delta }log\kappa ,\end{array}$$(28)

with the Δlog κ given by

$$\begin{array}{c}\hfill \mathrm{\Delta }log\kappa \simeq {\delta }_{T}log\kappa -{\left(\frac{\partial log{\kappa }^{\mathrm{GS}98}}{\partial log\rho }\right)}_{T,X}{\delta }_{T}log\rho .\end{array}$$(29)

Here, ${\delta }_{T}log\kappa \equiv {\kappa }^{\mathrm{GS}98}(\rho ,T,X)-{\kappa }^{\mathrm{A}09}(\stackrel{\sim }{\rho },T,\stackrel{\sim }{X})$ is computed at constant temperature, and a similar relation holds for δ_T log ρ. The result of this procedure is shown in Fig. 9, where it is seen that the required change climbs up to ∼30% around the base of the convective zone. For comparison, it also shows the induced opacity under the fifth force computed with the old, GS98 mixture. The latter, is orders of magnitude smaller than the former, and therefore does not interfere with the proposed opacity adjustment. Certainly, the metallicity problem does not constitute a systematic effect similar to opacity or diffusion, and should be rather seen as a choice we make upfront; the fifth-force effect on interior profiles would manifest in a similar manner, had we chosen to work with the opacity-adjusted A09 mixture, instead of the GS98 one.

Fig. 9. Illustration of the opacity adjustment to reconcile the two metallicity mixtures, as explained in Sect. 4.3. Solar models computed with the revised mixture (A09) are in tension with helioseismic inferences, when compared to those computed with the ‘old’ mixture GS98 (see Fig. 4, Table 2 and text). A proposed approach to resolve the tension is to consider an appropriate opacity adjustment (κ(new)) on top of the A09 profile (κ(A09)), as shown in the above (black, solid curve). The shaded blue band shows the typical uncertainty of ±5%, as discussed in Sect. 4. For comparison, the inset plot shows the induced change in the opacity profile under the fifth force computed with the GS98 metallicity mixture and unmodified opacity. The horizontal axis of log-temperature corresponds to a radial range between ∼0.0002 − 0.73 R⊙.

Recently, Magg et al. (2022) re-visited the solar metallicity problem. Their analysis suggests that the tension is reduced provided one adopts a newly proposed composition which is very close to the metal fractions of GS98. If this result is confirmed, it would imply that the solar metallicity problem is practically alleviated. For our analysis, these results support our choice of the GS98 metals fraction.

5. Constraints on the fifth-force coupling

Based on our previous equilibrium modelling, we now discuss how the two key solar regions, namely the radiative and convective zone, can be used to constrain generic gravity models. As we demonstrated earlier, each of these regions is impacted by input physics uncertainties in fundamentally different ways. This aspect is crucial, since it offers for complementarity of tests. Although our main interest is in the general scalar-field extensions of GR, our discussion easily extends to different theory setups, and highlight the significance of each solar region for precision tests of gravity.

5.1. Radiative zone

In this region, uncertainties from input physics tend to have a significant impact on interior profiles, most notably the sound speed as also seen from Fig. 5. However, as we showed, the existence of a narrow region where our relatively simple estimates of modelling systematics become minimal offers an ideal place to test gravity (see Fig. 6). What is more, helioseismic inversions yield highly accurate results in this region, which implies that helioseismically inferred interior profiles show similar level of accuracy.

Focussing on our particular theory setup parametrised by the fifth-force coupling 𝒴, we denote with c_s(⊙) the solar sound-speed profile inferred from inversions of helioseismic data (see also Fig. 10), and c_s(model) ≡ c_s(m)(𝒴) the profile predicted by the GS98 model with 𝒴 ≠ 0 at unmodified opacity and diffusion. To compare the model with the helioseismically inferred sound speed, we employ a simple likelihood functional defined as

$$\begin{array}{c}\hfill \mathcal{L}(\mathcal{Y})\propto e^{-{\chi }^2(\mathcal{Y})/2},\end{array}$$(30)

Fig. 10. Difference between the predicted (model) sound-speed profile and the profile derived from the inversion of helioseismic data (Sun), in the sense (Sun) – (model). The brown vertical bars show the errors on the inverted helioseismic data points, according to the results of Rhodes (1997). Errors have been multiplied by a factor of 10. The shaded band corresponds to the maximum uncertainty from our modelling of opacity and diffusion, while the continuous and dashed curves to a modified model computed with a fifth-force coupling of 𝒴 = 10−2 and 10−4 respectively. Fifth-force models are computed on top of our reference GS98 model at standard gravity (see Table 2). The fifth-force effect in the radiative zone is clearly seen for 𝒴 = 10−2, as expected from the model differences of Figs. 5 and 6.

with the χ² constructed with the fractional change of the sound speed as

$$\begin{array}{c}\hfill {\chi }^2(\mathcal{Y})={\displaystyle \sum _{r_i=r_1}^{r_2}}\frac{1}{{\sigma }_i^2}{(\frac{c_{\text{s}(\odot )}^{(}{r}_i)}{c_{\text{s}(\text{m})}^{(}{r}_i;\mathcal{Y})}-1)}^2.\end{array}$$(31)

We note that r_i denotes the radial points at which the χ² is evaluated, between the radial boundaries r₁, r₂. We choose not to interpolate the discreet inverted data for the sound-speed profile (c_s(⊙)), as derived from Rhodes et al. (1997). The helioseismic errors are denoted with σ_i. Indeed, given the strong correlation between solutions at neighbouring points (e.g., Howe & Thompson 1996) such interpolation would have little meaning; we neglect this correlation in the present preliminary analysis, but it should be taken into account in future more detailed work. We choose for the limits r₁, r₂ approximately the solar region in the radiative zone where our systematics of opacity and diffusion become minimal (see Figs. 6 and 10) (r₁, r₂) = (0.2, 0.3). In principle, any systematic uncertainties should be included in the statistical analysis and marginalised over in the end. In the particular solar region we are interested in though, the uncertainties from simple estimates of errors in the input physics are typically smaller than the helioseismic errors (see Fig. 5). For this solar region, we have eight data points for the inverted sound-speed profile available to us. Evaluating (30), we find that its maximum lies at 𝒴 ≃ −8 ⋅ 10⁻⁴. Integrating it as usual to find the 2σ confidence interval and assuming a flat prior on 𝒴, we find

$$\begin{array}{c}\hfill -{10}^{-3}\lesssim \mathcal{Y}\lesssim 5·{10}^{-4}.\end{array}$$(32)

This constraint on the fifth-force coupling strength for the general scalar-field theories improves previous bounds from stellar astrophysics by about three orders of magnitude. We also note that the fifth-force coupling 𝒴 relates to two dimensionless parameters associated with the dynamics of the original relativistic theory around a cosmological background, as 𝒴 = −2(α_H + β₁)²/(α_H + β₁), with β₁ quantifying the significance of higher-order scalar-metric kinetic operators in the cosmological action, and α_H the kinetic mixing between the scalar field and matter respectively (Dima & Vernizzi 2018). Therefore, constraints on 𝒴 translate to implications for the evolution of cosmological quantities and the formation of structures at large scales. This is, however, beyond the scope of present work.

It is now important to discuss and challenge our assumptions that led to the above result.

– We have assumed the use of the OP opacity tabulation, which along with the OPAL opacity table is one of the most advanced opacity tabulation to date, and modelled the opacity error to be constant throughout the star (= ± 5%). This is a popular choice in the literature and appears to be a somewhat conservative strategy. However, as we discussed in Sect. 4 (see also Fig. 7), one cannot exclude the scenario where the opacity uncertainty acquires a radial dependence which peaks around r ≃ 0.25 R_⊙, that is, in the solar region within the radiative zone where the fifth-force effect peaks. Such an effect would interfere with the fifth-force effect in the solar radiative zone. We would like to entertain this scenario here. In order to understand the consequences of it, we assume the existence of a localised opacity uncertainty between 0.2 − 0.3 R_⊙, which in turn induces a radially-dependent and localised change of the sound speed according to

$$\begin{array}{c}\hfill \frac{\delta c_{\mathrm{s}}}{c_{\mathrm{s}}}=\pm {10}^{-3}·e^{-\frac{{(r/R_{\odot }-0.25)}^2}{{(0.0005)}^2}},\end{array}$$(33)

with the σ = 5 × 10⁻⁴ chosen so that the uncertainty fades off sufficiently fast close to the edges at 0.2 and 0.3 R_⊙ respectively. The chosen amplitude is about an order of magnitude smaller than the average (constant) uncertainty of 5% we considered earlier, since it is already large enough to have a significant effect for our purpose. Notice that the + and − signs correspond to the case where the localised opacity uncertainty induces a positive (negative) effect on the sound speed. This is indeed a more general case that the one shown in Fig. 7. We proceed repeating our previous likelihood analysis by adding on top of our theoretically computed sound speed the correction (33). We find that for the ‘+’ and ‘−’ case in (33) the following constraints at 2σ; 1.8 × 10⁻³ ≲ 𝒴 ≲ 3.8 × 10⁻³, and −4 × 10⁻³ ≲ 𝒴 ≲ −2 × 10⁻³ respectively. Both constraints are in tension with Newtonian gravity (𝒴 = 0). This result is not surprising. A localised decrease (increase) of the sound speed, induced by a localised opacity uncertainty, requires a sufficiently positive (negative) fifth-force coupling 𝒴. Indeed, 𝒴 will tend to increase (decrease) the sound speed, compensating for the effect of the opacity change. In the presence of such a localised opacity uncertainty one would need to consider another solar region like the convection zone. We discuss this further below. From a statistical perspective, the opacity issue could be tackled through the exploration of the full opacity theory space along with that of the fifth force under a Monte-Carlo analysis. We comment further on this in Sect. 6.

– In this work, we have been interested in the solar equilibrium structure and the associated solar modelling systematics in the presence of the fifth force. In the absence of a proper helioseismic analysis in modified gravity, a key assumption we made is that the inverted solar sound-speed data we used as our ‘observations’ are computed at standard gravity. Indeed, helioseismic inversions infer the solar sound-speed profile from observed frequencies; however, being perturbative in nature they rely on the use of a background theoretical solar model. A known result in the literature is that the dependence of helioseismically inferred profiles have a weak dependence on the background standard solar model used (see e.g., Basu et al. 2000). We expect that this result would carry over to some extent to the case of the particular gravity model studied here, given the smallness of the fifth-force coupling. However, only a fully consistent helioseismic analysis would provide a robust confirmation of our results, towards self-consistent constraints. Such an analysis is certainly necessary, and it will require the computation of solar frequencies and the extension of helioseismic inversions in the presence of the fifth force⁵.

5.2. Convective zone

The solar convection zone offers yet another promising region to test deviations from Newtonian gravity. A key point here is that the scaling of systematics from input physics uncertainties is radically different compared to the radiative zone (see Fig. 5). As it was explained in detail in Sect. 2 using analytic relations, the reason behind this result is that the sound-speed profile in the convective zone is rather insensitive to the adiabatic structure of models with modified input physics. Therefore, the impact of variations under different input physics on the sound speed remains sufficiently small compared to the case of the radiative zone.

For our particular gravity model, the fifth force leaves a sizable signature in this region, as seen in Figs. 5 and 10. A drawback of testing gravity within the solar convective zone is that helioseismic inversions do not work as well as they do within the radiative zone due to contamination from surface effects. This, however, can be circumvented with updated helioseismic inversion analyses based on observed solar frequencies of increased precision. Such an effort would promote the solar convective zone to a powerful region for tests of alternative gravity theories. Nevertheless, we find it a useful exercise to derive a constraint on 𝒴 in the convective zone with our current helioseismic data which were derived assuming standard gravity. In this direction, we focus on the region 0.71 − 0.85 R_⊙ in order to minimise potential biases from surface effects. Repeating our previous analysis, we find a preference for stronger gravity and −10⁻³ ≲ 𝒴 ≲ 10⁻⁵ at 2σ. This constraint should not be taken seriously for the reasons explained earlier. A separate helioseismic analysis is needed in order to investigate the effect of the fifth force in the convective zone. However, the fact that the latter constraint is close to the one we found in the radiative zone provides some indication for the consistency of the analysis.

6. Summary and future work

We employed the Sun as our laboratory to test general extensions of GR. Although we focussed on the particular case of the DHOST scalar-field extensions of GR aiming to explain dark energy, the concept of our analysis easily applies to a broader set of gravity theories modifying the hydrostatic equilibrium structure. By means of solar evolution simulations aided by analytical results, we presented a quantitative description of how the fifth force affects the equilibrium solar structure of calibrated solar models, and explained its interplay with the delicate solar microphysics such as opacity, diffusion, EoS, and metallicity (see e.g., Fig. 5).

Despite the rather subtle competition between the fifth-force effect and the solar modelling systematics, we showed that the fifth force predicted by these general dark-energy theories leaves a sharp imprint on the solar sound-speed profile, with the peak of the effect located in a region where simple modelling uncertainties associated with opacity and diffusion tend to be insignificant. As we explained in Sect. 4 and showed in Fig. 6, the underlying mechanism for the latter result is due to the non-trivial scaling of the temperature and mean molecular weight in the solar radiative zone under the influence of the fifth force, providing a distinct phenomenological imprint on the sound speed. As illustrated in Fig. 5, the sound speed in the solar convection zone is weakly affected by the general uncertainties in standard solar modelling, while the fifth force has a direct and significant effect. This offers an interesting prospect for future robust tests of fifth forces stemming from extensions of the standard paradigm of gravity and particles. Current inversions for the solar sound speed in this region suffer from some uncertainty in this region owing to the effects of the inadequate modelling of the near-surface layers (e.g., Christensen-Dalsgaard et al. 1996), but this may be alleviated through the use of more extensive data or improved analysis techniques.

A comparison with a helioseismically inferred sound-speed profile allowed us to derive a bound on the fifth-force coupling 𝒴 which improves previous constraints from stellar physics by approximately 3 orders of magnitude (see Sect. 5). This result relied upon the use of a helioseismically inferred sound-speed profile assuming standard gravity. A full analysis of the fifth-force effects will involve consistent calculations of solar oscillation frequencies in the presence of the fifth force, followed by an inverse analysis using the resulting frequencies. In addition, to consistently take other uncertainties of solar modelling into account, a Monte-Carlo analysis should be carried out including these uncertainties. Such efforts will be reported in a future paper.

Our results highlight the power of the Sun as a laboratory for fundamental physics. In future, it would be very interesting to perform a helioseismic inversion analysis in the presence of the fifth force, extending the current machinery to the case of the fifth force. In turn, this will allow us to infer the sound-speed profile from observations in a fully consistent manner, and confirm the results presented in this work from a rather bottom-up helioseismic approach. Finally, it would be interesting to understand the impact of the fifth force on solar neutrinos and its confrontation with observations from Earth-based detectors.

Acknowledgments

Ippocratis D. Saltas is supported by the Czech Grant Agency (GACR), under the grant number 21-16583M. Funding for the Stellar Astrophysics Centre is provided by The Danish National Research Foundation (Grant DNRF106). It is a pleasure to acknowledge enlightening discussions with Earl Bellinger, Rasmus Handberg, Steen Hannestad, Stefan Ilic, Ilidio Lopes, Lorenzo Pizzuti. The authors would also like to thank the anonymous referee for constructive criticism and comments.

References

All Tables

All Figures

Fig. 1. Effective Newton’s coupling induced by the fifth force as a function of radius according to Eq. (1). The peak of the fifth-force effect around r ≃ 0.25 R⊙ in the radiative zone will be key for the phenomenology of solar interior profiles at hydrostatic equilibrium. The curves were evaluated according to Eq. (2) using the density profile of the reference model GS98 of Table 2, computed without the fifth force.

In the text

Fig. 2. The quantity W(r), defined in Eq. (19). This is a measure of the sound-speed gradient and a key diagnostic for the solar structure. Here, it is shown for the following three cases: for standard gravity (𝒴 = 0), weaker (𝒴 = 0.01) and stronger (𝒴 = −0.01) gravity respectively. The plots are based on full numerical simulations of calibrated models for radius, luminosity and surface metallicity (see Table 2 for a summary). Top: the solution of W(r) = − 2/3 yields an approximate estimate of the position of the radius of the base of the convective zone (r ≃ Rcz), which is shown in the upper panel. The reference model at standard gravity is model GS98 summarised in Table 2. Bottom: the peak of the quantity W(r) right below the solar surface due to the ionisation of helium. Although barely visible by eye, the position of the peak is shifted inwards (outwards) for stronger (weaker) gravity.

In the text

Fig. 3. Variation of total radius, base of convective zone and effective radius under the fifth-force for uncalibrated solar models. Each observable is normalised with respect to its value at 𝒴 = 0 (no fifth force) at a fixed initial helium abundance (Y0), initial metallicity (Z0) and mixing-length parameter (αMLT). Whereas the fifth force has a significant effect on the radius and position of the convective zone, its effect is milder on the effective temperature.

In the text

Fig. 4. Fractional difference for the sound speed in the sense (Sun) – (Model). ‘Sun’ refers to the sound-speed profile as inferred from helioseismic observations at standard gravity according to the results of Rhodes (1997). Models GS98 and A09 correspond to theoretically computed calibrated solar models at standard gravity (𝒴 = 0) with the two distinct, and most popular metallicity mixtures for the Sun (see also Table 2 and Sect. 3). For comparison, we also show two computed models with the fifth force turned on (𝒴 = ±0.01) based on the GS98 mixture. We note that the fifth-force effect peaks in the radiative zone which will be key for constraining the fifth force with helioseismic observations.

In the text

Fig. 5. Fractional model differences in the sense (reference model) – (modified model) for solar interior profiles. This illustrates the importance of the most significant uncertainties in solar modelling (opacity, diffusion) against the fifth-force effect. The legend for the curves is shown in the top left panel. The dot-dashed (dotted) curves correspond to an increase (decrease) of the respective input physics parameter. For the combined change of opacity and diffusion the dot-dashed curves correspond to increased opacity and decreased diffusion. The reference model is the calibrated solar model at standard gravity with GS98 metallicity and with unmodified opacity and diffusion coefficients, as explained in the text and Table 2. In all panels, the models with the fifth force turned on (black, solid curves) are computed on top of our reference model at standard gravity, and have a value of the fifth-force coupling |𝒴| = 10−2 for illustration purposes. As conservative estimates for the uncertainties in input physics we use ±5% for opacity and ±15% for diffusion, which are computed on top of our reference model (see Sect. 4). We also consider the combined effect of the simultaneous uncertainty of opacity and diffusion. The mild variation of the density and sound-speed profiles around the base of the convective zone (≃0.71 R⊙) has been explained through analytic relations in Sect. 2. (a) Density. (b) Speed of sound. (c) Temperature. (d) Adiabatic index. (e) Mean molecular weight. (f) Pressure.

In the text

Fig. 6. Top: model differences for the sound-speed profile similar to those of Fig. 5, in the sense (reference model) – (modified model), but with a smaller fifth-force coupling 𝒴 = 10−3. The plot zooms into the solar region where our estimates of systematics from opacity and diffusion may become negligible compared to that of the fifth force, in order to clearly illustrate the peak of the fifth-force effect in the radiative zone. The error bar shows approximately the error of helioseismic inversions at that point, ∼10−4, to highlight the precision down to which helioseismic reconstructions of the solar sound speed can constrain the coupling 𝒴. The insignificance of the opacity and diffusion uncertainties compared to the fifth force for this solar region is discussed in Sect. 4. Bottom: an illustration of the mechanism behind the ‘clean’ peak of the fifth-force effect on the speed of sound profile around r ≃ 0.25 R⊙, where the systematics from uncertainties related to opacity and diffusion become negligible compared to the effect of the fifth force (see also Fig. 5). The reference model at standard gravity is the GS98 model (see Table 2 and text). The plots assume for illustration the ideal-gas approximation for the sound-speed profile, where $c_{\text{s}}^2\propto {\Gamma }_{1}T/\mathrm{\mu }$, and neglect the small variation of the adiabatic index Γ1 (see also Fig. 5). Right panel: fractional difference of the sound-speed profile and its components under an increase of the opacity by 5%, but the picture is similar for uncertainties in diffusion. Left panel: how the picture changes for the fifth-force effect. For a detailed discussion we refer to the text around Eq. (27).

In the text

Fig. 7. Difference between solar profiles computed with different opacity tabulations providing an estimate of the error coming from opacity modelling. Here, we show the difference in the sound speed for a model computed with our reference opacity table (OP) and a model computed with an alternative opacity (OPAL). Both computations assume standard gravity (𝒴 = 0), and are based on the GS98 metallicity mixture. In our analysis of systematics, the uncertainty in the modelling of opacity was considered constant throughout the star, as explained in Sect. 4. However, one cannot exclude the possibility of a localised opacity variation such as the one emerging from the comparison between OP and OPAL opacity tabulations. For comparison purposes we also show the corresponding impact of the fifth force on the sound speed computed with our reference opacity (OP) and for two values of the fifth force coupling. We discuss the impact of this on inferences of the fifth-force coupling 𝒴 in Sect. 5.

In the text

Fig. 8. Comparison between the effect of a modification of the EoS at standard gravity, and that of the fifth force at fixed EoS on the sound-speed profile. The variation δcs is in the sense (reference model) – (modified model), and the vertical axis is on natural logarithmic scale. Our reference model uses the mixed OPAL/Free EoS, whereas the alternative EoS’s at standard gravity are based on pure Free EoS (‘modified EoS 1’), and pure OPAL (‘modified EoS 2’) respectively. For more details we refer to Sect. 4. The models with fifth force in the figure are computed on top of our reference EoS. We remind that the fifth-force effect on the sound-speed profile within the radiative zone peaks around r ≃ 0.25 R⊙ (see also Figs. 5 and 6).

In the text

Fig. 9. Illustration of the opacity adjustment to reconcile the two metallicity mixtures, as explained in Sect. 4.3. Solar models computed with the revised mixture (A09) are in tension with helioseismic inferences, when compared to those computed with the ‘old’ mixture GS98 (see Fig. 4, Table 2 and text). A proposed approach to resolve the tension is to consider an appropriate opacity adjustment (κ(new)) on top of the A09 profile (κ(A09)), as shown in the above (black, solid curve). The shaded blue band shows the typical uncertainty of ±5%, as discussed in Sect. 4. For comparison, the inset plot shows the induced change in the opacity profile under the fifth force computed with the GS98 metallicity mixture and unmodified opacity. The horizontal axis of log-temperature corresponds to a radial range between ∼0.0002 − 0.73 R⊙.

In the text

Fig. 10. Difference between the predicted (model) sound-speed profile and the profile derived from the inversion of helioseismic data (Sun), in the sense (Sun) – (model). The brown vertical bars show the errors on the inverted helioseismic data points, according to the results of Rhodes (1997). Errors have been multiplied by a factor of 10. The shaded band corresponds to the maximum uncertainty from our modelling of opacity and diffusion, while the continuous and dashed curves to a modified model computed with a fifth-force coupling of 𝒴 = 10−2 and 10−4 respectively. Fifth-force models are computed on top of our reference GS98 model at standard gravity (see Table 2). The fifth-force effect in the radiative zone is clearly seen for 𝒴 = 10−2, as expected from the model differences of Figs. 5 and 6.

In the text