Searching for dark energy with the Sun

General extensions of General Relativity (GR) based on bona fide degrees of freedom predict a fifth force which operates within massive objects, opening up an exciting opportunity to perform precision tests of gravity at stellar scales. Here, focusing on general scalar-tensor theories for dark energy, we utilize the Sun as our laboratory and search for imprints of the fifth-force effect on the solar equilibrium structure. With analytic results and numerical simulations, we explain how the different solar regions offer powerful ways to test gravity. Accounting for the delicate interplay between fifth force and solar microphysics such as opacity, diffusion, equation of state and metallicity, we demonstrate that the fifth force still leaves a sharp signature on the solar sound speed, in a region where simple estimates of input physics uncertainties become negligible. For general scalar-field extensions of GR, known as (U-)DHOST, based solely on the observational helioseismic errors, our analysis at the equilibrium level allows to place an approximate constraint on the fifth-force coupling strength of $-10^{-3} \lesssim \mathcal{Y} \lesssim 5\cdot 10^{-4}$ at $2\sigma$. This result improves previous stellar constraints by $\sim 3$ orders of magnitude, and should be confirmed and improved by future helioseismic inversions in modified gravity combined with an elaborate accounting of theoretical uncertainties. Our analysis can be applied to a wide set of theories beyond GR, and also paves the way for helioseismic analyses in this context. In this regard, we discuss how the solar radiative and convective zone can be employed as promising laboratories to test generic theories of gravity.


Introduction
The need to explain the acceleration of the Universe, and 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 impact on the formation and dynamics of large-scale structures in the Universe (see e.g. Amendola et al. 2017;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 is to 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 will 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 argued that helioseismic observations have the potential to improve on 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 will show how the fifth force leaves a strong signal in the solar interior, 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 weakfield 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) as where m(r) is the mass enclosed within a radius r from the centre and ρ(r), P (r) the density and pressure respectively. Within Newtonian gravity one obviously has G eff = G 0 , with G 0 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 with ≡ d/dr, and Y the dimensionless fifth-force coupling. The coupling Y relates parametrically to those param-Article number, page 1 of 14 arXiv:2205.14134v1 [astro-ph.SR] 27 May 2022 eters 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 Y is shown in Figure 1. The effect of the fifth force acquires its maximum value around r ∼ 0.3R , and it switches off at the solar surface where m (r) = 0.
It is interesting to note that a similar modification of gravity's strength is predicted in different theoretical setups (Olmo et al. (2020); ). 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.
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 −3 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 (Section 2), we study the effect of the fifth force on solar observables by means of intuitive arguments and analytic relations. In the second part (Sections 3, 4), we numerically investigate the evolution of the present Sun with numerical simulations, and analyse the interplay between 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 will 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 Section 5 we derive an order of magnitude constraint for the fifth-force coupling based on its effect on the solar sound speed, and in Section 6 we draw our conclusions.

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 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 1 See also Koyama & Sakstein (2015); Saridakis et al. (2021); . (2)), except for r 0.2, owing to the effects of sperical geometry (see also Eq. (22) below). The peak of the fifth-force effect around r 0.25R in the radiative zone will be key for the phenomenology of solar interior profiles at hydrostatic equilibrium. As expected, the fifth-force effect vanishes at the solar surface . The curves are evaluated according to Eq. (2) using the density profile of the reference model GS98 of Table 1, computed without fifth force. 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 central ∼ GM 2 /R 4 . This will be followed by a change in the temperature of similar magnitude, which can be estimated assuming a perfect gas EoS, P ∝ ρT . 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 will 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. To see this, let us write explicitly the equation of energy transport as with the gradient ∇ ≡ d log T /d log P . For radiative transport of energy, ∇ = ∇ rad ∝ κ(r) P (r)L(r)/(m(r)T 4 ), where κ is the opacity and L the luminosity. In view of the modified hydrostatic equilibirum (Eq. (1)), Eq.
(3) implies that i.e. a weakening of gravity (Y > 0, G eff /G 0 < 1) will tend to decrease the temperature gradient, and the opposite is true for stronger gravity (Y < 0, G eff /G 0 > 1). Therefore, a weaker (stronger) gravity will act as to make the energy transfer less (more) efficient. This argument serves only as a qualitative way towards gathering first intuition, since the fifth force affects hydrostatic equilibrium in a non-trivial, radius-dependent fashion. 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 pp ∝ X 2 H ρ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 will discuss later on.

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 where K is the adiabatic constant, Γ 1 , 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. (5) is approximately satisfied, with Γ 1 being the adiabatic compressibility, Γ 1 = (∂ log p/∂ log ρ) ad ; in the bulk of the Sun matter is approximately in the form of a fully ionized ideal gas, with Γ 1 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.3R , similar to the scaling exhibited by the effective gravity strength in Figure 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, shown in Figure 6. 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 B T , such that 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 Γ 1 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 δc 2 s /c 2 s in the radiative zone will be driven by the change in the temperature and mean molecular weight, We discuss this in detail in Section 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 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 Y = 0. Assuming m(r) M , we approximate the r.h.s of the equation of hydrostatic equilibrium as with the second equality due to the polytypic EoS. Assuming Γ 1 to be constant, as is approximately true except in the nearsurface ionization zones of hydrogen and helium, and neglecting the small correction on the r.h.s, we can straightforwardly integrate the above equation to find 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, Y = 0. Since the fifth-force coupling is small, we will adopt a perturbative approach. In particular, we start from the general hydrostatic Eq. (9), and assume a solution for the density profile in the form where ≡ πY is a small parameter 1, and ρ (Y) (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 Γ 1 , we now restrict ourselves to Γ 1 = 5/3, which corresponds to the choice of a polytropic index n = 3/2. Plugging Eq. (11) back into the hydrostatic Eq. (9), the -dependent part of the equation takes the form with ≡ d/dz, and P (z) ∝ 1/2 − 3z + (9/2)z 2 − 2z 3 , 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 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. (12) in the form ρ inhomog.
where Q a function to be determined. Substituting this ansatz in Eq. (12), one arrives at the differential equation Q = −P (z)/z 2 , 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. (12) is well-behaved on the solar surface, i.e. no divergence occurs as z → 1, one finally finds with the constant c defined as c ≡ Y 3π 10 and we also re-introduced our expansion parameter ≡ πY, appearing in Eq. (11). The solution (14) is understood as a small correction on top of the Y = 0 solution, ρ (0) , and it can be seen that it vanishes on the solar surface, i.e. 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. (10). We now define the fractional change of density as with ρ = ρ(r) understood as in Eq. (11). The definition (15) will serve as a definition for the fractional change of interior profiles under any change of physics. It is straightforward to check that δρ/ρ, like ρ (Y) , 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 −1 ≡ R /r around the base of the convective zone which we conventionally choose as R cz 0.71R , For the typical value of Y = 0.01, the above solution implies a very slow change in the convective zone, which is confirmed by our detailed solar evolutionary simulations shown in Figure 6.

Sound speed
We now look at the sound-speed profile defined through 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 (Y = 0) as reminding that we neglected corrections due to the nonconstancy 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 (18) 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 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 in Figure 6. A key observation from Eq. (18) 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/diffusion, in the convective zone. Similar to what we did before, we expand to first order around the base of the convective zone ( 0.71R ), assuming again Γ 1 = 5/3, to find This explains the mild variation of δc 2 s (r)/c 2 s 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 with the last approximate equality true within the convective zone. In the convective zone, W 1 − Γ 1 , with Γ 1 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 Γ 1 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 δΓ 1 /Γ 1 ∼ O(10 −5 ), and therefore the effect is negligible. In Figure 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 (Y < 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 Y > 0. In other words, a stronger gravity implies a steeper sound-speed profile.

Numerical method
We use the 1-dimensional MESA stellar evolution code (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015. The code solves the system of stellar equations discretely at each time step on a 1-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 will 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 ex-  (21), 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 (Y = 0), weaker (Y = 0.01) and stronger (Y = −0.01) gravity respectively. The plots are based on full numerical simulations of calibrated models for radius, luminosity and surface metallicity (see Table 1 for a summary). 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 1. As it can be seen from the upper panel, a weaker (stronger) gravity tends to make the convective zone shallower (deeper). 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.
press as which is the form which we also adopt for the numerical implementation. Note that in the above equation dm/dr = 4πρ(r)r 2 . The derivatives in Eq. (22) are implemented in a discreet manner as where k labels the k-th cell on the grid and all physical variables are evaluated either at the centre or the edge of the k−th 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 7 years. Our microphysics assumptions are as follows. For the EoS we use OPAL (Rogers & Nayfonov 2002), whereas the opacities  Table 1: A description of solar models used in our analysis, and their initial parameters. All are calibrated to reproduce the observed solar surface metallicity, radius and effective temperature within 1σ of the respective error (see also Table 2), and have standard (unmodified) opacity and diffusion coefficients. They are computed with our reference EoS (OPAL EoS), except the model computed with the alternative choice of Free EoS. The first model is our reference model at standard gravity and standard input physics (opacity, diffusion, etc.), and is based on the "old" metallicity mixture of Grevesse & Sauval (1998), which is in good agreement with helioseisimic observations. Model A09 uses the revised metallicity of Asplund et al. (2009), which is, however, in tension with helioseismic data (see Fig. 4). The latter model is not used for computations with the fifth force, but was computed for comparison with GS98. The models in the table computed with the fifth force turned on (Y = ±0.01) are used for illustration purposes throughout the text. More details about our choice of input physics are provided in Section 3.

Observable
Value ±1σ Solar radius [R ] 1 0.713 ± 10 −3 Table 2: A summary of the solar observables and their precision, used in this work. The quoted value for the radius of the base of the convective zone, Rcz, is the one inferred from helioseismic observations (e.g., Christensen-Dalsgaard et al. 1991;Basu & Antia 1997). The old (GS98) and the revised (A09) surface metallicity is according to Grevesse & Sauval (1998) and Asplund et al. (2009) respectively. are based on the tabulations from OP (Seaton 2005) with lowand 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.

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  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 Y = 0 (no fifth force) at a fixed initial helium abundance (Y0), initial metallicity (Z0) and mixinglength 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. One notices that the weakening of gravity for most part of the star, i.e Y > 0, leads to an increase in the total radius and vice versa. Similarly, a weakening of gravity further tends to make the convective zone shallower. . "Sun" refers to the sound-speed profile as inferred from helioseismic observations at standard gravity according to the results of Rhodes et al. (1997). Models GS98 and A09 correspond to theoretically computed calibrated solar models at standard gravity (Y = 0) with the two distinct, and most popular metallicity mixtures for the Sun (see also Table 1, and Section 3). The discrepancy between the predictions of the GS98 and A09 model, when compared to helioseismic observations, is what 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 opacities in the solar interior. For comparison, we also show two computed models with the fifth force turned on (Y = ±0.01) based on the GS98 mixture. We note that the fifth-force effect peaks in the radiative zone (see also Figure 6) which will be key for constraining the fifth force with helioseismic observations.
solar models, i.e 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 observ-ables are the total stellar radius (R), surface element abundances ((Z/X) surf. ) and luminosity (L) (see also Table 2). 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 0 ) abundance and initial metallicity (Z 0 ). On top of the previous parameters, we need to define the value of the fifth-force coupling Y. Therefore, our parameter space is 4dimensional. 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 fifth force, i.e. Y = 0. This enables us to pin down the 3-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 χ 2 functional built out of the desired observables (O i ), Typically, for a solar calibration, 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 section 4.3.) The set of representative calibrated solar models we use in our analysis are summarised in Table 1. We proceed turning on the fifth-force term using our GS98 calibration as our reference model, which extends our parameter space to 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), At fixed mass and radius, and for Y > 0, if the luminosity L is to remain the same, the strengthening of gravity in the solar core (see Figure 1) requires the compensating effect of a decreased µ. Figure 3 shows the effect of fifth force on some solar observables, for uncalibrated models. One see that the fifth force predominantly affects the total radius and the thickness of the convective envelope, and to a smaller extent temperature.  Figure 6, in the sense (reference model) -(modified model), but with a smaller fifth-force coupling Y = 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 Y. The insignificance of the opacity/diffusion uncertainties compared to the fifth force for this solar region is discussed in Section 4 and further illustrated in Figure 6. Bottom: An illustration of the mechanism behind the "clean" peak of the fifth-force effect on the speed of sound profile around r 0.25R , where the systematics from uncertainties related to opacity and diffusion become negligible compared to the effect of the fifth force (see also Figure 6). The reference model at standard gravity is the GS98 model (see Table 1 and text). The plots assume for illustration the ideal-gas approximation for the sound-speed profile, where c 2 s ∝ Γ1T /µ, and neglect the small variation of the adiabatic index Γ1 (see also Figure 6). The right panel shows the 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. The almost zero variation of the sound-speed profile around r 0.25R , i.e δc 2 s 0, is due to the fine cancellation between the contribution of the mean molecular weight and temperature. The left panel shows how the picture changes for the fifth-force effect -the induced change in the mean molecular weight under the fifth force, thanks to the change in the hydrogen fraction, in combination with the non-trivial scaling of the temperature profile, leads to the domination of the fifth-force effect compared to opacity/diffusion uncertainties in this solar region.

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 will use our numerical results to attack this issue. We will 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.  Table 1. 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 |Y| = 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 Section 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.71R ) has been explained through analytic relations in Section 2.
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 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 with κ 0 the reference opacity profile based on the OP tabulation, and the chemical mixture GS98 (see also Table 1). 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 will therefore use for the uncertainty on diffusion, 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/diffusion), are shown in the model differences of Figure 6. Here the effect on the thermodynamic properties of the change in composition is illustrated in terms of the mean molecular weight, approximately calculated as 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 fifthforce coupling Y > 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 Y < 0. The increase of density, pressure and temperature before the turn over at r 0.2R in Figure 6, 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 (10) and (18), the structure generally depends on the constant K which is different for each modified model. This directly affects density (see Eq. (10)) and hence pressure, while Eq. (18) shows that the sound speed is independent of K. In contrast,Eq. (20) shows that the fifth force directly affects the sound speed in the convection zone, leading to the behaviour seen in Figure 6.
The peak around r = 0.25R in the fifth-force effect on density, pressure and sound speed shown in Figure 6 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. (7)). Then, where we neglected the insignificant variation of Γ 1 under model differences compared to the other terms (see also the model differences of Figure 6). 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.25R . This is not anymore true for model variations under the fifth force (at unmodified opacity/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 Figure 5. This surprisingly clean signal of the fifth force allows for precision tests of the theory with helioseismic observations. The comparison of indicative computed models under the fifth force against helioseismic data is shown in Figure 10.
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 e.g. Figure 30 of Christensen-Dalsgaard (2021), and it is illustrated in Figure 8. The induced difference in the sound-speed profiles computed with OPAL and OP, respectively, shows a peak at r 0.3R , 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 Section 5 in light of new constraints on Y. 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 we discuss in Section 4.3 (see also Figure 4) the proposed metallicity mixtures A09 shows tension with helioseismology, although the very recent analysis of Magg et al. (2022) suggests that the tension is significantly settled through a revised mixture, close to the GS98 one.

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 Γ 1 , 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 Figure 6. It is evident that fifth-force effect on Γ 1 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 a solar model at standard gravity and GS98 metallicity, but with the so-called Free EoS, instead of the OPAL which we used as our choice of reference. Figure 7 shows a comparison between the choice of EoS at standard gravity and that of the fifth force at fixed EoS (OPAL) for the sound-speed profile. For Y = 10 −2 the fifth-force effect on the sound speed is significantly larger than that of the EoS deep in the radiative zone by a couple of orders of magnitude. However, for Y = 10 −4 the two effects become almost indistinguishable. Clearly, the choice of the EoS for sufficiently small values of Y matters. This is, however, not an issue of concern. Typically, a choice for the EoS is made upfront in stellar modelling. As such, any EoS describing sufficiently well the relevant microphysics should not be considered to contribute to a systematic modelling uncertainty. Furthermore, the excellent agreement between our reference EoS (OPAL) and the alternative choice of Free EoS suggests the consistency of our choice of EoS, that is, 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.

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 Figure 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 2 if we adjust the opacity profile of the A09 model κ through a radius-dependent opacity correction ∆κ according to log κ new (ρ, T, X i ) = log κ(ρ, T, X i ) + ∆ log κ, 2 As argued in those works, this proposal is mostly empirical and to large extent it lacks physical motivation.  Fig. 7: Comparison between the effect of a different EoS choice at standard gravity, and that of the fifth force at fixed EoS on the sound-speed profile. The variation δc 2 s is in the sense (reference model) -(modified model). Our reference model uses the OPAL EoS, whereas the alternative EoS choice at standard gravity is the so-called Free EoS (see Section 4). The models with fifth force in the figure are computed on top of our reference model with OPAL EoS. It is seen that for a fifth-force coupling Y = 10 −2 , the effect of a different EoS on the sound-speed profile remains significantly smaller than that of the fifth force, while for Y = 10 −4 the two effects start becoming indistinguishable. We remind that the fifth-force effect on the sound-speed profile within the radiative zone peaks around r 0.25R (see also Figures 6, 5 and 10). Although the choice of high-level EoS does not constitute a systematic uncertainty in stellar modelling, the above result suggests that for sufficiently small values of Y, the choice of EoS matters. We emphasize that the good agreement between our reference EoS and the alternative one (black, continuous curve) highlights the consistency of our reference EoS in modelling the solar microphysics.
with the ∆ log κ given by Here, δ T log κ ≡ κ GS98 (ρ, T, X) − κ A09 (ρ, T,X) is computed at constant temperature, and a similar relation holds for δ T log ρ. The result of this procedure is shown in Figure 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.
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.

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 convec- , 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 Section 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. The localised opacity variation around r 0.3R could interfere with the signature of the fifth force on the sound-speed profile. 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 Y in Section 5. tive 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.

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 Figure 6. 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 Figure 5). What is more, helioseismic inversions yield highly accurate results in this region, which implies that helioseismically inferred interior profiles show similar level of accuracy.
Focusing on our particular theory setup parametrised by the fifth-force coupling Y, we denote with c 2 s( ) the solar soundspeed profile inferred from inversions of helioseismic data (see also Figure 10), and c 2 s(model) ≡ c 2 s(m) (Y) the profile predicted by the GS98 model with Y = 0 at unmodified opacity/diffusion. To compare the model with the helioseismically inferred sound speed, we will employ a simple likelihood functional defined as  Fig. 9: An illustration of the opacity adjustment to reconcile the two metallicity mixtures, as explained in Section 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 Figure 4, Table 1 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 Section 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. Clearly, the induced opacity variation under the fifth force is significantly smaller and cannot interfere with the metallicity issue in any sense. The horizontal axis of log-temperature corresponds to a radial range between ∼ 0.0002 − 0.73R . The recent results of Magg et al. (2022) suggest the solar metallicity issue is alleviated through the use of a revised metallicity mixture which is very close to the GS98 one.  Rhodes et al. (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 Y = 10 −2 and 10 −4 respectively. Fifth-force models are computed on top of our reference GS98 model at standard gravity (see Table 1). The fifth-force effect in the radiative zone is clearly seen for Y = 10 −2 , as expected from the model differences of Figures 6 and 5. Evidently, models computed with a fifth-force coupling Y 10 −4 start being in agreement with helioseismic predictions for the sound-speed profile.
with the χ 2 constructed with the fractional change of the sound speed as r i denotes the radial points at which the χ 2 is evaluated, between the radial boundaries r 1 , r 2 . We choose not to interpolate the discreet inverted data for the sound-speed profile (c 2 s( ) ), as derived from Rhodes et al. (1997). 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 1 , r 2 approximately the solar region in the radiative zone where our systematics of opacity and diffusion become minimal (see Figures 5 and 10) (r 1 , r 2 ) = (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 Figure 6). For this solar region, we have eight data points for the inverted sound-speed profile available to us. Evaluating (32), we find that its maximum lies at Y −8 · 10 −4 . Integrating it as usual to find the 2σ confidence interval and assuming a flat prior on Y, we find 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.
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 Section 4 (see also Figure 8), one cannot exclude the scenario where the opacity uncertainty acquires a radial dependence which peaks around r 0.3R , 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.3R , which in turn induces a radiallydependent and localised change of the sound speed according to with the σ = 5 · 10 −4 chosen so that the uncertainty fades off sufficiently fast close to the edges at 0.2 and 0.3R 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 Figure 8. We proceed repeating our previous likelihood analysis by adding on top of our theoretically computed sound speed the correction (35). We find that for the "+" and "-" case in (35) the following constraints at 2σ; 1.8·10 −3 Y 3.8·10 −3 , and −4·10 −3 Y −2·10 −3 respectively. Both constraints are in tension with Newtonian gravity (Y = 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 Y. Indeed, Y 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 Section 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 selfconsistent 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.

Convective zone
The solar convection zone offers yet another promising region to test deviation 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 Figure 6). As it was explained in detail in Section 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 Figures 6 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 in-creased 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 Y 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.85R in order to minimise potential biases from surface effects. Repeating our previous analysis, we find a preference for stronger gravity and −10 −3 Y 10 −5 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.

Summary and future work
We employed the Sun as our laboratory to test general extensions of General Relativity. Although we focused on the case of scalar-field extensions of GR aiming to explain dark energy, our analysis easily applies to a broader set of gravity theories. 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, equation of state and metallicity (see e.g. Figure 6).
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 Section 4 and showed in Figure 5, the underlying mechanism for the latter result is due to the nontrivial scaling of 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 Figure 6, 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 Y which improves by approximately 3 orders of magnitude previous constraints from stellar physics (see Section 5). This result relied upon the use of a helioseismically inferred sound-speed profile assuming standard gravity. A full analysis of the fifthforce 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 take other uncertainties of solar modelling consistently 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.