Stellar Oscillations in Modified Gravity
Abstract
Starting from the equations of modified gravity hydrodynamics, we derive the equations of motion governing linear, adiabatic, radial perturbations of stars in scalartensor theories. There are two new features: first, the eigenvalue equation for the period of stellar oscillations is modified such that the eigenfrequencies are always larger than predicted by General Relativity. Second, the General Relativity condition for stellar instability is altered so that the adiabatic index can fall below before unstable modes appear. Stars are more stable in modified gravity theories. Specialising to the case of chameleonlike theories, we investigate these effects numerically using both polytropic LaneEmden stars and models coming from modified gravity stellar structure simulations. We find that the change in the oscillation period of Cepheid star models can be as large as 30% for orderone matter couplings and the change in the inferred distance using the periodluminosity relation can be up to three times larger than if one had only considered the modified equilibrium structure. We discuss the implications of these results for recent and upcoming astrophysical tests and estimate that previous methods can produce new constraints such that the modifications are screened in regions of Newtonian potential of .
I Introduction
Modified theories of gravity (MG) have received a lot of attention in the last decade. Upcoming surveys such as Euclid Amendola et al. (2012) and LSST Ivezic et al. (2008) will test Einstein’s theory of General Relativity (GR) on cosmological scales to incredibly high precision. On the theoretical side, the cosmological constant problem is driving the study of alternate theories of gravity in an attempt to explain why the observed vacuum energy density is 120 orders of magnitude smaller than the quantum field theory prediction (see Clifton et al. (2012) for a review). Cosmology probes the longrange, infrared (IR) behaviour of gravity and so the majority of these theories are IR modifications, often involving new scalardegrees of freedom which either drive the acceleration, screen out the effects of the cosmological constant or include some symmetry such that its value is technically natural. GR is the unique Lorentzinvariant theory of a massless spin2 particle Weinberg (1965) and so any modification of GR necessarily introduces at least one new degree of freedom. These could be decoupled from matter, however, in any generic theory there is no symmetry forbidding this and so any decoupling is unnatural. If a coupling is present then one must face the issues of additional or fifthforces that result whenever a fieldgradient is present. GR has been tested to incredibly high precision over the last century and so one must finetune the strength of these couplings to incredibly small values in order to have a viable theory. Such small couplings usually render any modifications insignificant in all but the strongest gravitational fields. For example, any quintessence model of dark energy requires a mass for the scalar to be of order Copeland et al. (2006) so that the forcerange is of the order of the observable universe. If such a field is then coupled to matter it must pass laboratory tests, which require the force to be submm Will (2004). Therefore, the only way to have a light quintessencescalar coupled to matter is to tune all of the couplings to very small values. In general, there is no symmetry protecting these values and so quantum corrections are expected to induce orderone corrections. The small value of these couplings is therefore tantamount to finetuning.
There is one caveat to the above argument: all of the tests of GR have been performed in our local neighbourhood and so there is nothing preventing a theory which has large fifthforces active over cosmological scales provided that these are somehow absent locally. Such theories are said to posses screening mechanisms and there are several examples in the literature including the chameleon mechanism Khoury and Weltman (2004a, b), the symmetron mechanism Hinterbichler and Khoury (2010), the environmentdependent DamourPolyakov effect Brax et al. (2010), the Vainshtein mechanism Vainshtein (1972) (see also Nicolis et al. (2009); Babichev et al. (2009); Hinterbichler (2012) for applications to modern theories such as Galileons and massive gravity) and the disformal mechanism Koivisto et al. (2012). Most of the theories that include a screening mechanism have the property that the fieldprofile of the new degrees of freedom are sourced by the ambient density. They are constructed such that at high densities either the field’s masses are large enough that the force is shortranged, the effective coupling to matter becomes negligible or the fieldgradients are suppressed with respect to the gradient of the Newtonian potential. The exception is the disformal mechanism where field gradients are sourced by pressure, which is generally negligible for nonrelativistic sources and so no fieldgradient is induced.
The problem is then changed from reconciling fifthforces with local experiments to distinguishing these theories from GR observationally. Theories such as Galileons show interesting effects cosmologically through either the modified background expansion rate, linear perturbation properties or halo clustering effects Barreira et al. (2013a, b, c), however, the same is not true of chameleons and symmetrons, where the cosmological mass of the field is always less than Brax et al. (2012a); Wang et al. (2012), precluding any effects on linear scales. On nonlinear scales, there can still be some galaxy clustering effects (see, for example Hu and Sawicki (2007); Li et al. (2012); Brax et al. (2012b, 2013a) and references therein). Astrophysical tests of these theories can probe regions of parameter space where all cosmological signatures vanish and therefore have the potential to constrain them by several orders of magnitude below what is possible with cosmological tests.
Several astrophysical tests of different modified gravity theories have been investigated, including stellarstructure tests Chang and Hui (2010); Davis et al. (2012); Jain et al. (2012), black hole tests Hui and Nicolis (2012) and galactic dynamics tests Jain and VanderPlas (2011); Lee et al. (2012); Vikram et al. (2013). Here, we will be concerned only with stellar structure. The purpose of this work is to present and investigate a consistent framework for dealing with the equations of modified gravity scalar fields coupled to hydrodynamics at zeroth and firstorder in perturbation theory with an aim to looking at both current and new observational tests. The effects of disformally coupled scalars on the equilibrium structure of nonrelativistic stars are still unknown and a calculation is beyond the scope of this work. We hence focus on theories where the field is coupled to the trace of the energy momentum tensor only and the fifthforce is proportional the the field gradient and the strength of the matter coupling. This includes chameleonlike and Vainshtein screened theories. Starting from the modified equations of hydrodynamics, we derive the equations governing the equilibrium structure of stars (the modified stellarstructure equations) and the evolution of linearised, radial, adiabatic perturbations of both the field and stellar structure. We find three new effects in modified gravity:

When the star is unscreened, the period of oscillation can be reduced by an amount. For Cepheid stars, which are good probes of MG Jain et al. (2012), we find that the change in the period can result in differences from the GR inferred distance measurement of up to three times what has been previously been calculated using an approximation and therefore that current bounds can be improved using the same data sets.

When the star is unscreened, the star is more stable to these perturbations. There is a wellknown result in stellar astrophysics that when the adiabatic index falls below the squaredfrequency of the fundamental mode is negative and so there is an unstable mode. We will show that in MG the index can fall below without any instability. In MG, there is no universal bound on the critical index for the appearance of the instability; its precise value depends on the structure and composition of the star as well as how unscreened it is.

The perturbations of the star can source scalar radiation and vice versa.
The first two effects require the star to be unscreened while the third does not. Vainshteinscreened theories are too efficient to leave any star unscreened and so it is only the third effect that may be used to probe this theory. In this work, we will only investigate the first two effects and so for this reason, we will only deal explicitly with chameleonlike theories, although we will present the general derivation and comment briefly on the application to Galileonlike theories.
Chameleonlike theories screen according to the local Newtonian potential; objects where this is lower than a universal model parameter (to be defined below) are unscreened whereas they are screened when the converse is true. One is then led to look for MG effects in the most underdense regions of the universe. In terms of the Newtonian potential, these are dwarf galaxies. Dwarf galaxies located in clusters are screened by the Newtonian potential of their neighbours whilst isolated dwarf galaxies may show MG effects. Recently, MG Nbody codes have been combined with SDSS data to produce a screening map of the nearby universe Cabre et al. (2012). This has allowed dwarf galaxy tests of MG using current data to be performed for the first time Jain et al. (2012); Vikram et al. (2013).
In unscreened dwarf galaxies located in cosmological voids, main and postmainsequence stars may become partially unscreened leading to strong fifthforces in their outer layers. The equilibrium structure of stars in chameleonlike theories has been well studied Chang and Hui (2010); Davis et al. (2012). These stars must provide a stronger outward pressure gradient in order to combat these additional forces than they would had they been GR stars. The extra energy per unit time required to source this gradient comes from an increased rate of nuclear burning in the core. A star under the influence of MG is therefore more luminous, more compact and has a higher effective temperature than its GR doppelgänger. A modified version of the publicly available code MESA, which can simulate the structure and evolution of stars in these theories was presented in Davis et al. (2012) and has become a powerful tool in the quantitative prediction of astrophysical effects. At the level of perturbations, Cepheidvariable stars pulsate at a higher frequency and so obey a different periodluminosity (PL) relation to the GR one. Using a simple approximation to find the new period, Jain et al. (2012) have used MESA models to place the strongest constraints in the literature to date. The MG frequencies can be numerically computed by solving the new equations coming from modified gravity hydrodynamics and one aim of this work is to compare the new MG predictions with this approximation. The others are to investigate the stability of stars in these theories as well as the implications for current and possible future astrophysical tests.
This paper is then organised as follows: in section II, we present the action for chameleonlike theories and outline some of the essential features. We then examine the screening mechanism and present the two modelindependent parameters that can be used to fully specify the effects of MG on the structure and evolution of stars. We then give a brief introduction to astrophysical tests of modified gravity and outline the current constraints. In section III we derive the equations of modified gravity hydrodynamics. At zerothorder, these give the modified equations of stellar structure. At firstorder, we derive the equations of motion describing the coupled linear, radial, adiabatic perturbations of both the field and radial displacements. In section IV we specialise to the case where the scalarradiation is negligible and elucidate the properties of the new equation of motion for stellar perturbations. In section V we show how the criterion for the star to be stable to these perturbations is altered in and argue that stars are more stable in MG. In section VI, we present numerical results. First, we investigate LaneEmden models of convective stars and calculate the new oscillation periods and critical value of the adiabatic index. We then turn our attention to MESA models. We calculate the change in the period and inferred distance from the PL relation for the same Cepheid models used in Jain et al. (2012) in order to asses the accuracy of their approximation and find that when the full effects of hydrodynamic perturbations are included, the differences from GR can be up to three times larger than the approximation predicts. Using this, we estimate the new constraints that can be found using the same datasets and procedure as Jain et al. (2012). These correspond to constraints where the modifications are screened in regions of Newtonian potential . Finally, we discuss our results in light of recent and future astrophysical tests and conclude in section VII.
Ii Screened ScalarTensor Gravity
ii.1 The Action
Any theory which contains a scalar gravitationallycoupled to matter fields will generically result in modifications of the Newtonian forcelaw in the nonrelativistic limit and therefore constitute a modified theory of gravity. One wellstudied class of model has this coupling arise because the metric in the matter Lagrangian differs from the metric appearing in the Ricci scalar by a conformal factor such that
(1) 
Typically, so that the theories do not differ too greatly between frames. These theories are described by the general action
(2)  
(3) 
and give rise to a Newtonian fifthforce (per unit mass)
(4) 
is known as the coupling. Theories with screening mechanisms generally
have the property that in the
unscreened region
Before going on to discuss the highdensity screening mechanism, it is important to address the notion of a conserved matter density. The conformal coupling of the field to the visible sector has the result that the energymomentum tensor of matter is not covariantly conserved, rather one has , where is the trace of the energymomentum tensor. The density does not obey the usual continuity equation and one generally works with the conserved density instead. This conserved density is also fieldindependent. From here on we shall refer to as the conserved matter density corresponding, for example, to the stellar mass per unit volume.
ii.2 The Screening Mechanism
The coupling to matter has the effect that the scalar potential appearing in the action (2) is not the same scalar potential which governs the field’s dynamics. Instead, the equation of motion is
(5) 
corresponding to a densitydependent effective potential
(6) 
It is this densitydependence which is responsible for screening the fifthforce in overdense environments. Suppose that this potential has a minimum. The fieldposition of this minimum can differ by many orders of magnitude between the solarsystem and the cosmological background (recall ) and so it is possible to engineer effective potentials where either the mass of small oscillations about the highdensity minimum is large such that the range of the fifthforce is subm Khoury and Weltman (2004b) or the coupling and the magnitude of the force is negligible Hinterbichler and Khoury (2010); Brax et al. (2010). In this way, the fifthforce is undetectable in highdensity areas.
One can then ask exactly when an object will be screened. Consider a spherically symmetric overdensity with profile and radius , which is a small perturbation about a lower but more spatiallyextensive constant density whose characteristic lengthscale is . For example, a galaxy in the cosmological background or a star inside a galaxy. If is large enough such that the field can reach its minimum at the higher density and remain there over most of the object’s interior then the fifthforce will be negligible and the object is screened. If the converse is true and the field deviates only from the value in the background by a small perturbation then the full effects of the modification of GR will manifest. In the static, nonrelativistic limit we can ignore timederivatives and so the equation of motion (5) is:
(7) 
There are then two cases to consider. First, if the field can reach its minimum we have and equation (7) is unsourced. In this case the field is just constant and equal to the highdensity minimum value so that . When this solution is valid, there are no fifthforces but in general there is some radius , known as the screening radius, where the field begins to move away from this value and we expect large fifthforces in this region. In this second case, the field is a small perturbation about the underdense value so that . Subtracting the equation of motion in both the over and underdense regions and linearising we have
(8) 
where , is the mass of the field in the background, and we have used the fact that . In all theories of interest, we have so that the fifthforce gives rise to novel features on large scales and so we have ignored the first term, it is negligible compared with the Laplacian. Now we know that is related to the Newtonian potential via the Poisson equation, and so we may substitute this into (8) and integrate twice to find the field profile:
(9) 
where is the Heaviside stepfunction. We refer to the region as the screened region and the region as unscreened. Using the definition of the Newtonian potential , where is the mass enclosed within a sphere of radius (the total mass of the overdensity is ), the fifthforce (4) in the unscreened region is given by the derivative of (9):
(10) 
where . We will assume that there
are no fifthforces when and that the force is given by equation
(10) when from here on
All that remains is to find the screening radius . Using the fact that when as does , equation (9) gives us an implicit expression for the screening radius:
(11) 
It will be useful later to recast this as an integral equation:
(12) 
The parameter is known as the selfscreening parameter and is of paramount importance to the screening properties of these theories. Since and , there are no solutions when
(13) 
In this case, the object is fully screened and . When the object will be at least partially unscreened.
The screening and fifthforce properties in any region are fully specified
by and , however
these are very environmentdependent and are nonlinearly related to the field
values in different regions of the universe. The exception to this is unscreened
objects, where the field is only a small perturbation around the value in the
background and so these values are roughly constant. Astrophysical tests of MG
Davis et al. (2012); Jain and VanderPlas (2011); Cabre et al. (2012); Jain et al. (2012); Vikram et al. (2013) are
performed by comparing the properties of screened and unscreened galaxies.
Whether or not a galaxy is screened or not depends on its Newtonian potential
relative to the cosmological value of ; the same is true of denser
objects in unscreened galaxies. In each case, the strength of the fifthforce in
any unscreened region is proportional to the cosmological value of and so astrophysical tests of MG probe the cosmological values of
and . The theory is then parametrised in a
modelindependent manner by the cosmological values of these parameters, which
we denote by , which measures the strength of the fifthforce in
unscreened regions and , which determines how small the Newtonian
potential must be in order for the object to be unscreened
ii.3 Astrophysical Tests of Modified Gravity
We will not attempt any comparison with data here but for completeness and motivation we will describe the basic premise behind observational tests using astrophysical effects below.
Current Constraints
Only unscreened objects will show deviations from GR and so we must look for such objects in the universe if we wish to observe these deviations. There are two classes of screened objects: selfscreened, where the selfNewtonian potential of the object is larger than and environmentscreened, where the object would be unscreened in isolation but is still screened due to the Newtonian potential of its neighbours (see Hui et al. (2009) for a nice discussion of equivalence principle violation effects and their relation to environmental screening). The cosmological bound coming from cluster constraints is Schmidt et al. (2009). The Newtonian potential of spiral galaxies is and so this would leave most galaxies, including our own milky way selfunscreened. Clearly if this were the case then our own galaxy would show deviations from GR, which are not observed and so we must conclude that our own galaxy is screened, either because or because we are screened by the Newtonian potential of the local group. While this has been debated in the past, recently, an independent constraint from comparing water maser distances to NGC 4258 with tip of the red giant branch distances has placed the bound Jain et al. (2012). Finally, the strongest constraints on come from comparing tip of the red giant branch distances with Cepheid distances and place the bound Jain et al. (2012).
Dwarf Galaxy Tests
With this bound, there are very few unscreened objects in the universe. There
are no cosmological signatures, either at the level of the background or linear
perturbations (although see Brax et al. (2013b, c) for an
unusual model with a very strong coupling to matter) and spiral galaxies
are screened. The only selfunscreened objects in the universe are dwarf
galaxies, whose Newtonian potentials are and large
postmainsequence stars
Observational signatures using dwarf galaxies generally fall into two classes: modifications of the galaxy dynamics Jain and VanderPlas (2011) and those resulting from effects on stellar structure Chang and Hui (2010); Davis et al. (2012). Both have been carried out Jain et al. (2012); Vikram et al. (2013) using SDSS data and, currently, the former is not competitive with the latter. In this work we are interested in stellar structure effects and so we shall only describe these here.
The most practical method for using stellar structure tests is to compare different distance indicators to the same galaxy. The distance to a galaxy using stellar phenomena is either calibrated on Milky Way (or local group galaxies in some cases) objects or derived from the GR theoretical model. Either way, if the stellar property used to calculate this distance is altered in MG then the measured distance to an unscreened galaxy will differ from the GR distance. If, on the other hand, the indicator is insensitive to MG then any observation will give precisely the GR distance. If one then compared two distances to the same galaxy, one using a screened distance indicator and the other using an unscreened one, any discrepancy probes MG. By using many different galaxies and comparing to a screened sample, a statistically significant constraint can be achieved. This is what was done in Jain et al. (2012).
Cepheid Distances in Modified Gravity
In practice, the best unscreened distance indicators are Cepheid variable stars. They have Newtonian potentials of and there is a large amount of data available.
Cepheid stars are 510 stars which have evolved off the mainsequence and onto the red giant branch. During this phase, their HertzsprungRussell (HR) diagrams show a socalled blue loop, where the temperature rapidly increases and decreases at approximately constant luminosity. Whilst traversing the loop, the tracks cross the instability strip, where the star is unstable to Cepheid pulsations. Along this narrow strip in the  plane, the star is unstable due to a layer of partially ionised helium coinciding with the region where the motion becomes nonadiabatic. Upon compression, the energy does not go into raising the pressure, which would result in an increased outward force giving a stabilising effect but rather into further ionising the helium. This damming up of the energy acts like an engine and drives large changes, up to 10%, in the stellar radius. Along this strip, there is a wellknown period luminosity relation:
(14) 
where is the bolometric magnitude and , , and are empirically determined constants. Empirically, Freedman and Madore (2010). Now in MG, and we will see this derived below, the period of a Cepheid at fixed luminosity, effective temperature and mass change by some amount so that this empiricallyfitted relation predicts a deviation from the GR distance
(15) 
Previously Jain et al. (2012), was calculated using the approximation , where is calculated using some appropriate average. In this work, we will derive the equations governing its value and show how it can be calculated numerically using some simple examples and compare our results with this approximation.
Iii Modified Gravity Hydrodynamics
In this section we shall start from the modified equations of hydrodynamics and proceed to derive the equations governing both the equilibrium structure and linearised radial perturbations of nonrelativistic stars.
We will describe bulk quantities such as the pressure and density in the Eulerian picture, where these quantities are to be considered as fields which give the value of said quantities at any point in space as a function of time. In contrast, we will describe the position of individual fluid elements (and, when needed, the pressure perturbations) in the Lagrangian picture, where the motion of individual fluid elements are followed as a function of time. In this case, the Lagrangian position of a fluid element satisfies the momentum equation:
(16) 
where is the pressure, is the density and is the external force density. In MG, the fluid moves under its own Newtonian gravity and the fifthforce (4) due to the scalar field so that
(17) 
This is the only hydrodynamical equation that is altered relative to GR;
changing gravity only changes the motion of the fluid elements and does not
directly
(18) 
which may be integrated once to give
(19) 
Since mass is a locallyconserved quantity we also have the continuity equation:
(20) 
where is the velocity of the fluid element. In general, one must also consider the energy generation and radiative transfer equations but these are only important if one wishes to study the effects of perturbations coupled to stellar atmospheres, which is irrelevant in the context of modified gravity. We will include their effects when simulating the equilibrium stellar configuration numerically in order to produce the correct stellar properties, however, we will not include them in our perturbation analysis. Instead, we will work in the socalled adiabatic approximation, where the density and pressure evolve according to
(21) 
The quantity
(22) 
is the first adiabatic index or equation of state. It is of paramount importance to the study of stellar pulsation and stability and we will return to discuss it later on. The simplest gases are typically characterised by a constant equation of state so that , where () for nonrelativistic (relativistic) gases. In general, one must couple the hydrodynamic equations to the full radiative transfer system and extract the equation of state from the result.
iii.1 Equilibrium Structure
We are interested in small perturbations about the equilibrium structure of stars and so we must first find the modified zerothorder equations. These were presented heuristically in Davis et al. (2012), however, here we shall derive them from this more formal setup. The equilibrium stellar configuration is both static and spherically symmetric and can be found by setting timederivatives to zero and in the hydrodynamic equations so that this now represents the Eulerian coordinate. We will denote all equilibrium quantities with a subscriptzero except for , which is defined at the background level only. It is important to note that is not a property of the star but is the cosmological value of found by evaluating (11) using the cosmological values of and . In what follows, is the equilibrium fieldprofile throughout the star and not the cosmological value. With no timedependence, (20) is trivially satisfied and . This simple form of the density profile allows us to find the mass enclosed in any given radius:
(23) 
The momentum equation (17) then reduces to the modified hydrostatic equilibrium equation
(24) 
Physically, this equation describes the pressure profile the star must assume in order for the star to support itself against gravitational collapse. The second term is the fifthforce due to the scalar field; stars in modified gravity need to provide larger pressure gradients in order to combat this extra inward component Chang and Hui (2010); Davis et al. (2012); Jain et al. (2012). These equations are then supplemented by the radiative transfer equation
(25) 
which describes how the temperature varies due to the flux of energy with luminosity away from regions of energy generation governed by
(26) 
Here, is the opacity, the crosssection for radiation absorption per unit mass and is the energy generation rate per unit mass from processes such as nuclear burning.
Taken by themselves, these equations do not close and one must specify the equations of state relating and , which are themselves determined by further equations involving energy transfer and nuclear burning networks. In the subsequent sections, we will solve these equations analytically by assuming different equations of state. When computing the equilibrium structure using MESA, lookup tables are used to find the opacity and pressure at a given density and temperature while energy generation rates are found using nuclear burning networks for individual elements. We refer to reader to the MESA instrumentation papers Paxton et al. (2011, 2013) for the full details.
iii.2 Linear Perturbations
The dynamics of stellar oscillations are governed by the hydrodynamics of small perturbations about the equilibrium configuration and so we shall linearise the equations (17), (18), (20) about some assumed background profile ignoring secondorder terms in the perturbations.
The fundamental object of interest is the linearised perturbation to the
Lagrangian position of each fluid element , which
describes the oscillation of the fluid from equilibrium at each radius. So far,
our treatment has been completely general and we have worked with the full
threedimensional problem. Continuing in this fashion would result in a
complete treatment of both radial and nonradial modes of oscillation. The aim
of this work is to provide a consistent framework with which to predict the
oscillation properties of partially unscreened stars residing in unscreened
dwarf galaxies. The extragalactic nature of these stars ensures that only their
fundamental radial mode (and possibly the firstovertone) are
observable and so from here on we will specify to the case of radial
oscillations so that is a purely radial vector
(27) 
and the goal of the present section is to derive is equation of motion.
We shall work with Eulerian perturbations of the background profile, which we will distinguish from Lagrangian perturbations by the use of a tilde so that
(28)  
(29)  
(30)  
(31) 
As remarked above, the Lagrangian perturbations may provide more physical insight on occasion and the two are related, for example, via
(32)  
(33) 
It will be useful later to have the Lagrangian pressure perturbation in terms of our systemvariables. This is given by
(34) 
where the first equality holds only in the adiabatic approximation.
We begin by perturbing equation (20) to find
(35) 
This may be used in the perturbed form of (17) to
find
(36) 
The perturbed Poisson equation is
(37) 
which may be integrated once using (35) to find
(38) 
Now
(39) 
and using (21) this is
(40) 
We wish to eliminate and and so we substitute (35) into (40) to find
(41) 
We then have
(42) 
This may then be used with (37) and (36) to find:
(43) 
This is the equation of motion governing the evolution of . In GR, and we obtain a single equation. In MG, however, we need a separate equation for , which is found by perturbing the equation of motion (5):
(44) 
where is the mass of the unperturbed field at zero density and equation (35) has been used. We are interested in stationarywave solutions and so we expand
(45)  
(46) 
to yield two coupled equations
(47)  
(48) 
Equations (47) and (48) constitute the main result of this section. One could combine them into a single equation, however, it is more instructive to treat the system as two coupled equations. In GR, we have and (48) reduces to
(49) 
which describes linear, adiabatic, radial waves moving in the stellar interior. It is known at the linear adiabatic wave equation (LAWE) and its eigenfrequencies give the frequency of stellar oscillations about equilibrium. We will hence refer to (48) as the modified linear adiabatic wave equation (MLAWE). Its properties will be the subject of the next section.
Using the profile (10), we have
(50) 
which we shall use in all analytic and numerical computations from here on.
iii.3 Application to Other Modified Gravity Theories
Before proceeding, it is worth noting that although this work focuses on chameleonlike theories, the hydrodynamic analysis performed so far is more general than this. Any scalartensor theory whose Einstein frame matter coupling can be written in the effective form
(51) 
where is the trace of the matter energymomentum tensor, will give rise to an additional force
(52) 
provided that is the canonically normalised field
Iv Properties of the Modified Linear Adiabatic Wave Equation
The MLAWE describes the behaviour of stellar oscillations under modified gravity. There are two major differences with respect to the GR equation. Firstly, there are two additional terms in the homogeneous part, proportional to the first and second derivatives of the background field. When these are negligible and the homogeneous part behaves as it would in GR, however, these are comparable to the other terms when and encode the effect of modified gravity on wave propagation in the region exterior to the screening radius. Physically, the term proportional to acts as a varying enhancement of Newton’s constant given by equation (10) and the term proportional to can schematically be viewed as and so it encodes the effect of a varying Newtonian force in the outer regions.
The second effect is a driving term proportional to . This is clearly the effect of the fifthforce due to perturbations in the field. This was modeled by Silvestri (2011) as an inhomogeneous forcing term at a single frequency. Here it appears as coupling between the field and stellar perturbations: the stellar perturbations source the scalar field perturbations and vica versa. Physically, corresponds to scalar radiation (or rather the flux at infinity). There is evidence from previous studies Silvestri (2011); Upadhye and Steffen (2013) that this is negligible in chameleonlike systems and so from here on we will neglect the dynamics of the field perturbations and treat only the homogeneous part of the MLAWE (48). There may, nevertheless, be interesting features associated with the coupled oscillatorlike problem such as beating and scalar radiation and this is left for future work.
iv.1 Boundary Conditions
The MLAWE requires two boundary conditions in order to fully specify the
solution given a specific value of
(53) 
The surface boundary condition depends on the stellar atmosphere model (see Cox (1980) for a discussion) but the lowest modes, where the period of oscillation is far longer than the inertial response time of the atmosphere can be described by solutions with so that . This gives the surface condition Cox (1980)
(54) 
where the Lagrangian pressure perturbation is given by (34). Note that the additional terms in the MLAWE vanish at the stellar centre and radius if we take so that these conditions are identical to GR.
iv.2 SturmLiouville Nature of the Equation
The MLAWE can be written in SturmLiouville form
(55) 
where the weight function and the operator can be written
(56)  
(57) 
The problem of finding the pulsation frequencies is then one of finding the eigenvalues of these equations that correspond to eigenfunctions satisfying the boundary conditions (53) and (54). In practice, it is not possible to solve these equations analytically for physically realistic stars and numerical methods must be used. We will do just this in section VI. Despite the need for numerics, a lot of the new modified gravity features can be discerned and elucidated using wellknown SturmLiouville techniques and so we shall investigate these first.
iv.3 Scaling Behaviour of the Eigenfrequencies
Using the dimensionless quantities:
(58)  
(59)  
(60) 
the MLAWE (48) can be cast into dimensionless form:
(61)  
(62) 
where
(63) 
is the dimensionless eigenfrequency and the term proportional to is only present when . In GR, and one can solve this given some equilibrium stellar model to find . Since this must be a dimensionless number one has . In MG, there are two sources of change to this value at fixed mass: the change due to the different equilibrium structure described by (24) and the change due to the additional term in the MLAWE. At the level of the background, we expect that stars of fixed mass have smaller radii and larger values of (where by we mean some appropriate average over the entire star) so that the frequencies are higher in MG and at the level of perturbations, one can replace in the GR equation by the effective frequency so that and we therefore expect the MG eigenfrequency to be larger still.
The behaviour and importance of each effect requires a full numerical treatment, however, one can gain some insight by considering scaling relations when the star is fully unscreened so that . The equations of stellar structure are selfsimilar so that we can replace each quantity, such as the pressure, by some characteristic quantity, in this example the central pressure , multiplied by some dimensionless function in order to see how these different characteristic quantities scale when others are varied. Let us assume an equation of state of the form
(64) 
where is a constant and differs from since the system need not be adiabatic. In this case, equations (23) and (24) give
(65)  
(66) 
which can be combined to find the scaling of the radius for a fullyunscreened star in MG:
(67) 
at fixed mass. Ignoring the MG perturbations, one would then expect
(68) 
We shall confirm this limit numerically for some simple models later in section VI.1. In the fully unscreened limit, we would then expect the eigenfrequencies to scale approximately like
(69) 
so that they are always larger than the GR prediction, at least when .
V Stellar Stability
Given the SturmLiouville nature of the problem, we can find an upper bound on the fundamental frequency using the variational principle. Given an arbitrary trial function , one can construct the functional
(70) 
which has the property that . Ignoring MG for now and taking the simplest case where is constant, the fundamental eigenfrequencies of the LAWE (III.2) satisfy
(71) 
When the right hand side is negative we have and the
eigenfunctions have growing modes. This is the wellknown result in stellar
astrophysics that stars where the adiabatic index falls below are unstable
to linear perturbations and cannot exist
In modified gravity, we have
(72) 
and so this stability condition is altered in stars which are at least
partially unscreened. This is not surprising given the form of equation
(48). The term proportional to the derivative of
behaves like a position dependent mass for , which is negative when
so that one would expect growing modes. The additional terms in
(72) are due to the two new terms in (48), which are
of precisely the same varying mass form with the opposite sign. A negative mass,
which would signify an instability, coming from the GR term can then be
compensated by the new terms in MG, restoring stability
Physically, the adiabatic index is a measure of how the pressure responds to a compression of the star. Given a compression from one radius to a smaller radius , a larger adiabatic index will result in more outward pressure. If this increase in pressure is faster than the increase in the gravitational force, the star can resist the compression and is hence stable. Below the critical value of , the converse is true and the star is unstable. We have already seen above that the new terms contributing to the stability correspond to a varying value of and its derivative in the outer layers. Since modified gravity enhances the gravity, one may näively expect that its effect is to destabilise stars, however, we will argue below that this is not the case. The MLAWE describes acoustic waves propagating in the star. If the gravity and its gradient is larger in the outer layers then it is more difficult for these waves to propagate and hence modes which would usually have been unstable are stabilised.
Once again, we must disentangle the effects of the modified equilibrium structure and the perturbations on the critical value of . Consider first the modified equilibrium structure only. In this case, the stability condition is given by the GR expression, equation (71), however the pressure and density profiles will be different. Clearly, the critical value for the instability is still since this is the only value which makes the integral vanish but this does not necessarily mean the stability is altered away from this value. Scaling the pressure and density using the dimensionless quantities defined in (58), (59) and (60), we have
(73) 
where is a dimensionless function which depends on the composition of
the star
Let us now turn our attention to the effects of perturbations. Using equation (50), we have
(74) 
The additional term is clearly positive and so one may lower the value of below and still find positive eigenfrequencies, confirming our earlier intuition that the effect of MG is to stabilise stars compared with GR. Unlike GR, there is no universal critical index in MG and a numerical treatment is necessary to extract the value for a star of given composition and screening radius. We will investigate the stability of some simple semianalytic models in section VI.1.2 but before doing so, one can gain some insight into the full effects of MG on the stability by using the same scaling relations as section IV.3. In particular, we can set the numerator in (74) to zero to find the modified critical index:
(75) 
where is a dimensionless number which encodes the effects of the structure and composition of the star. In theories, Brax et al. (2008) and so we expect the critical value of to change by assuming that . We will verify numerically that this is indeed the case for a simple model in section VI.1.2.
Vi Numerical Results
We will now proceed to solve the MLAWE for various different stars. We will do this for two different stellar models: LaneEmden models, which describe the equilibrium configuration of spheres of gas collapsing under their own gravity, and MESA models, which are fully consistent and physical models that include effects such as nuclear burning networks, metallicity, convection, energy generation and timeevolution. The first models are simple compared to the second but they have the advantage that the nongravitational physics (e.g. nuclear burning) is absent, which will allow us to gain a lot of physical intuition about the new modified gravity features. They are simple semianalytic models and this allows us to first investigate the MLAWE using a controlled system with known scaling properties and limits without the complications arising from things like radiative transfer. This also allows us to test that the code is working correctly since we can compare our results with both the GR case, which has been calculated previously, and the fully unscreened case, which can be predicted analytically given the GR one. Their perturbations can also be described using an arbitrary value of , independent of their composition and so we will use them to study the modifications to stellar stability. These models are not realistic enough to compare with observational data and the power of MESA lies in that it can produce realistic models of stars such as mainsequence and Cepheid stars, which will allow us to predict the effects of MG on realistic stars in unscreened galaxies. MESA predictions have already been used to obtain the strongest constraints on chameleonlike models Jain et al. (2012) and combining these models with the modified gravity oscillation theory has the potential to provide new constraints.
The details of the implementation of MG into MESA is given in Appendix A (see also Davis et al. (2012); Jain et al. (2012) for some previous uses of the same code) and, where necessary, details of the numerical procedure used to solve the MLAWE are given in Appendix B. The shooting method has been used to solve the MLAWE in all instances.
vi.1 LaneEmden Models
LaneEmden models describe spheres of gas that support themselves against gravitational collapse by producing an outward pressure given by the solution of the (in our case modified) hydrostatic equilibrium equation. The key assumption is that the radiation entropy per unit mass is constant, which decouples the effects of nuclear burning (26) and radiative transfer (25) from the stellar structure equations derived in section (III.1). In this case, only the massconservation equation (23) and the hydrostatic equilibrium equation (24) are relevant. In GR, these equations are selfsimilar, which allows one to scale out the physical quantities in favour of dimensionless ones but this property is lost if one assumes a modified gravity profile of the form (10) and so, in order to retain a comparison with the GR case, we will solve the modified hydrostatic equilibrium equation
(76) 
which physically corresponds to setting in the region exterior to the screening radius. LaneEmden models describe polytropic gases defined by
(77) 
where is constant for a star of given mass but varies between different
stars. , is known as the polytropic index. The case was
studied in Davis et al. (2012); here we shall generalise their procedure to
arbitrary values. We can eliminate in terms of the LaneEmden coordinate
(78) 
where and are the central pressure and density respectively. Next, we rewrite the pressure and density in terms of a dimensionless function such that and . Substituting this into the mass conservation and modified hydrostatic equilibrium equations (equations (23) and (24) respectively) we can combine to two to find the modified LaneEmden Equation:
(79) 
where is the LaneEmden screening radius such that . The boundary conditions for this equation are () and (this follows from the fact that in the hydrostatic equilibrium equation). The radius of the star in LaneEmden coordinates, is then found from the condition (the pressure is zero at the stellar radius), from which the physical radius, may be found. For convenience, we introduce the quantities
(80)  
(81) 
In GR, and there is a unique solution for any given value of . We will denote the GR values of and using and respectively. In MG, there is a twodimensional space of solutions at fixed given by specific values of and , each with different values of , , and .
By writing equation (12) in LaneEmden variables, we find an implicit relation for and hence the screening radius:
(82) 
where