A Implementation of Modified Gravity into MESA

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 scalar-tensor 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 chameleon-like theories, we investigate these effects numerically using both polytropic Lane-Emden 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 order-one matter couplings and the change in the inferred distance using the period-luminosity 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 up-coming 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 long-range, infra-red (IR) behaviour of gravity and so the majority of these theories are IR modifications, often involving new scalar-degrees 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 Lorentz-invariant theory of a massless spin-2 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 fifth-forces that result whenever a field-gradient is present. GR has been tested to incredibly high precision over the last century and so one must fine-tune 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 force-range 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 sub-mm Will (2004). Therefore, the only way to have a light quintessence-scalar 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 order-one corrections. The small value of these couplings is therefore tantamount to fine-tuning.

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 fifth-forces 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 environment-dependent Damour-Polyakov 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 field-profile 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 short-ranged, the effective coupling to matter becomes negligible or the field-gradients 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 non-relativistic sources and so no field-gradient is induced.

The problem is then changed from reconciling fifth-forces 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 non-linear 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 stellar-structure 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 frame-work for dealing with the equations of modified gravity scalar fields coupled to hydrodynamics at zeroth- and first-order 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 non-relativistic 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 fifth-force is proportional the the field gradient and the strength of the matter coupling. This includes chameleon-like and Vainshtein screened theories. Starting from the modified equations of hydrodynamics, we derive the equations governing the equilibrium structure of stars (the modified stellar-structure equations) and the evolution of linearised, radial, adiabatic perturbations of both the field and stellar structure. We find three new effects in modified gravity:

  1. 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.

  2. When the star is unscreened, the star is more stable to these perturbations. There is a well-known result in stellar astrophysics that when the adiabatic index falls below the squared-frequency 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.

  3. 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. Vainshtein-screened 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 chameleon-like theories, although we will present the general derivation and comment briefly on the application to Galileon-like theories.

Chameleon-like 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 under-dense 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 N-body 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 post-main-sequence stars may become partially unscreened leading to strong fifth-forces in their outer layers. The equilibrium structure of stars in chameleon-like 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, Cepheid-variable stars pulsate at a higher frequency and so obey a different period-luminosity (P-L) 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 chameleon-like theories and outline some of the essential features. We then examine the screening mechanism and present the two model-independent 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 zeroth-order, these give the modified equations of stellar structure. At first-order, 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 scalar-radiation 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 Lane-Emden 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 P-L 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 data-sets 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 Scalar-Tensor Gravity

ii.1 The Action

Any theory which contains a scalar gravitationally-coupled to matter fields will generically result in modifications of the Newtonian force-law in the non-relativistic limit and therefore constitute a modified theory of gravity. One well-studied 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 fifth-force (per unit mass)

(4)

is known as the coupling. Theories with screening mechanisms generally have the property that in the unscreened region 1. This form of the action is referred to as the Einstein frame action since the Einstein frame metric, , is used to compute the Ricci scalar and not the Jordan frame metric . One may instead work in the Jordan Frame, in which case the action can be obtained using the Weyl rescaling (1). In what follows we will work exclusively in the Einstein frame since it allows for intuitive physical interpretations of the new features of stellar structure and perturbation theory. A subset of these theories has received a lot of attention recently because certain choices of the scalar potential and conformal factor can lead to local screening of the fifth-force (4) through either the chameleon mechanism Khoury and Weltman (2004a, b), the symmetron mechanism Hinterbichler and Khoury (2010) or the environment-dependent Damour-Polyakov effect Brax et al. (2010). In what follows, we will use a model-independent parametrisation of the theory parameters which is well-suited for astrophysical tests and refer the reader to Jain and Khoury (2010); Khoury (2010); Davis et al. (2012); Brax et al. (2012a); Brax (2012) for the specific details of individual screening mechanisms.

Before going on to discuss the high-density 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 energy-momentum tensor of matter is not covariantly conserved, rather one has , where is the trace of the energy-momentum tensor. The density does not obey the usual continuity equation and one generally works with the conserved density instead. This conserved density is also field-independent. 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 density-dependent effective potential

(6)

It is this density-dependence which is responsible for screening the fifth-force in over-dense environments. Suppose that this potential has a minimum. The field-position of this minimum can differ by many orders of magnitude between the solar-system and the cosmological background (recall ) and so it is possible to engineer effective potentials where either the mass of small oscillations about the high-density minimum is large such that the range of the fifth-force is sub-m 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 fifth-force is undetectable in high-density areas.

One can then ask exactly when an object will be screened. Consider a spherically symmetric over-density with profile and radius , which is a small perturbation about a lower but more spatially-extensive constant density whose characteristic length-scale 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 fifth-force 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, non-relativistic limit we can ignore time-derivatives 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 high-density minimum value so that . When this solution is valid, there are no fifth-forces 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 fifth-forces in this region. In this second case, the field is a small perturbation about the under-dense value so that . Subtracting the equation of motion in both the over and under-dense 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 fifth-force 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 step-function. 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 over-density is ), the fifth-force (4) in the unscreened region is given by the derivative of (9):

(10)

where . We will assume that there are no fifth-forces when and that the force is given by equation (10) when from here on2. When there is no fifth-force and the object is screened whereas when the object is full unscreened and we have , where is the Newtonian force. The parameter therefore determines the strength of the fifth-force.

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 self-screening 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 fifth-force properties in any region are fully specified by and , however these are very environment-dependent and are non-linearly 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 fifth-force 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 model-independent manner by the cosmological values of these parameters, which we denote by , which measures the strength of the fifth-force in unscreened regions and , which determines how small the Newtonian potential must be in order for the object to be unscreened3. We will work with this parametrisation in what follows. theories are chameleon theories with Brax et al. (2008) and is equivalent to the derivative of the function today, .

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: self-screened, where the self-Newtonian potential of the object is larger than and environment-screened, 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 self-unscreened. 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 self-unscreened objects in the universe are dwarf galaxies, whose Newtonian potentials are and large post-main-sequence stars 4. Isolated dwarf galaxies located in voids should therefore show deviations from GR compared to their counterparts located in clusters. The basic premise then is to compare the observational properties altered by MG in a sample of screened and unscreened dwarf galaxies. Any statistically significant discrepancy is a sign of MG whereas an agreement within experimental errors constrains and . One then needs some method of discerning whether or not an observed galaxy is screened. Recently, the results of MG N-body simulations have been used in conjunction with Sloan Digital Sky Survey (SDSS) data to provide a screening map of the universe for different values of Cabre et al. (2012). Using this map, it has already been possible to place new constraints on these theories Jain et al. (2012); Vikram et al. (2013).

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 5-10 stars which have evolved off the main-sequence and onto the red giant branch. During this phase, their Hertzsprung-Russell (HR) diagrams show a so-called 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 non-adiabatic. 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 well-known 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 empirically-fitted 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 non-relativistic 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 fifth-force (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 directly5 alter other processes such as mass conservation, energy generation and radiation flux. The quantity () is the mass enclosed inside a radius from the centre, and is given via the Poisson equation

(18)

which may be integrated once to give

(19)

Since mass is a locally-conserved 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 so-called 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 non-relativistic (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 zeroth-order equations. These were presented heuristically in Davis et al. (2012), however, here we shall derive them from this more formal set-up. The equilibrium stellar configuration is both static and spherically symmetric and can be found by setting time-derivatives to zero and in the hydrodynamic equations so that this now represents the Eulerian coordinate. We will denote all equilibrium quantities with a subscript-zero 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 field-profile throughout the star and not the cosmological value. With no time-dependence, (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 fifth-force 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 cross-section 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, look-up 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 second-order 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 three-dimensional problem. Continuing in this fashion would result in a complete treatment of both radial and non-radial 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 extra-galactic nature of these stars ensures that only their fundamental radial mode (and possibly the first-overtone) are observable and so from here on we will specify to the case of radial oscillations so that is a purely radial vector6. A more convenient quantity to work with is the relative displacement, given by

(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 system-variables. 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 find7

(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 stationary-wave 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 chameleon-like theories, the hydrodynamic analysis performed so far is more general than this. Any scalar-tensor theory whose Einstein frame matter coupling can be written in the effective form

(51)

where is the trace of the matter energy-momentum tensor, will give rise to an additional force

(52)

provided that is the canonically normalised field 8. The only assumption we have made so far is but this was not necessary, only appropriate, and if one has a scalar-tensor theory where this is not the case then one would find additional terms in (48). In terms of the theories studied so far we have . Another class of theories that has attracted recent interest are those which screen via the Vainshtein mechanism Vainshtein (1972) such as Galileons Nicolis et al. (2009) and massive gravity de Rham et al. (2011); Hinterbichler (2012), which are described by where is a constant describing the strength of the fifth-force. In this case, one has so that . The MLAWE (48) therefore applies equally to Vainshtein screened theories provided that one provides the counterpart to the equation of motion for the field perturbation (47). Vainshtein screening is too efficient to show any modified gravity effects on stellar structure at the background level and so one would expect the terms proportional to derivatives of to be negligible and hence that the oscillation frequencies are identical to GR, however, this is not necessarily true at the level of perturbations and scalar radiation may provide an observational test. We will not investigate this possibility here but postpone it for future work.

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 fifth-force 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 chameleon-like 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 oscillator-like 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 9. Firstly, our system is spherically symmetric and so we must impose at . The MLAWE then requires

(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 Sturm-Liouville Nature of the Equation

The MLAWE can be written in Sturm-Liouville 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 well-known Sturm-Liouville 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 self-similar 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 fully-unscreened 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 Sturm-Liouville 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 well-known result in stellar astrophysics that stars where the adiabatic index falls below are unstable to linear perturbations and cannot exist10.

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 stability11.

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 star12. In MG, the radius of the star will be smaller than its GR counterpart and the effective value of is larger. Hence, when the maximum possible frequency is greater than in GR whereas when the maximum frequency is more negative. If a star is unstable in GR then MG enhances the instability, moving the frequency further away from zero. This can also be seen from the scaling relation (68). If then is even more negative. At the background level, the effects of MG are to destabilise stars that are already unstable, without altering the stability condition.

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 semi-analytic 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: Lane-Emden 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 time-evolution. The first models are simple compared to the second but they have the advantage that the non-gravitational 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 semi-analytic 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 main-sequence 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 chameleon-like 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 Lane-Emden Models

Lane-Emden 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 mass-conservation equation (23) and the hydrostatic equilibrium equation (24) are relevant. In GR, these equations are self-similar, 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. Lane-Emden 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 Lane-Emden coordinate 13 where is defined by

(78)

where and are the central pressure and density respectively. Next, we re-write 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 Lane-Emden Equation:

(79)

where is the Lane-Emden 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 Lane-Emden 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 two-dimensional space of solutions at fixed given by specific values of and , each with different values of , , and .

By writing equation (12) in Lane-Emden variables, we find an implicit relation for and hence the screening radius:

(82)

where