{'fraction_non_alphanumeric': 0.046861638447066874, 'fraction_numerical': 0.013382931417821105, 'mean_word_length': 4.188796450360511, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stability of spherically symmetric solutions in modified theories of gravity

In recent years, a number of alternative theories of gravity have been proposed as possible resolutions of certain cosmological problems or as toy models for possible but heretofore unobserved effects. However, the implications of such theories for the stability of structures such as stars have not been fully investigated. We use our"generalized variational principle", described in a previous work, to analyze the stability of static spherically symmetric solutions to spherically symmetric perturbations in three such alternative theories: Carroll et al.'s f(R) gravity, Jacobson&Mattingly's"Einstein-aether theory", and Bekenstein's TeVeS. We find that in the presence of matter, f(R) gravity is highly unstable; that the stability conditions for spherically symmetric curved vacuum Einstein-aether backgrounds are the same as those for linearized stability about flat spacetime, with one exceptional case; and that the"kinetic terms"of vacuum TeVeS are indefinite in a curved background, leading to an instability.


I. INTRODUCTION
The idea of using a variational principle to bound the spectrum of an operator is familiar to anyone who has taken an undergraduate course in quantum mechanics: given a Hamiltonian H whose spectrum is bounded below by E 0 , we must have for any normalized state |ψ in our Hilbert space E 0 ≤ ψ|H|ψ . We can thus obtain an upper bound on the ground state energy of the system by plugging in various "test functions" |ψ , allowing these to vary, and finding out how low we can make the expectation value of H.
This technique is not particular to quantum mechanics; rather, it is a statement about the properties of selfadjoint operators on a Hilbert space. In particular, we can make a similar statement in the context of linear field theories. Suppose we have a second-order linear field theory, dependent on some background fields and some dynamical fields ψ α , whose equations of motion can be put into the form for some linear spatial differential operator T . Suppose, further, that we can find an inner product (·, ·) such that T is self-adjoint under this inner product, i.e., for all ψ 1 and ψ 2 . Then the spectrum of T will correspond to the squares of the frequencies of oscillation of the system; and we can obtain information about the lower bound of this spectrum via a variational principle of the form (3) * Electronic address: seifert@uchicago.edu Moreover, if we find that we can make this quantity negative for some test function ψ, then the system is unstable [2]; and we can obtain an upper bound on the timescale of this instability by finding how negative ω 2 0 must be. One might think that this would be a simple way to analyze any linear (or linearized) field theory; in practice, however, it is easier said than done. A serious problem arises when we consider linearizations of covariant field theories. In such theories, the linearized equations of motion are not all of the form (1); instead, we obtain non-deterministic equations (due to gauge freedom) and equations which are not evolution equations but instead relate the initial data for the fields ψ α (since covariant field theories on a spacetime become constrained when decomposed into "space + time".) Further, even if we have equations of the form (1), it is not always evident how to find an inner product under which the time-evolution operator T is self-adjoint (or even if such an inner product exists.) Thus, it would seem that this method is of limited utility in the context of linearized field theories.
In spite of these difficulties, Chandrasekhar successfully derived a variational principle to analyse the stability of spherically symmetric solutions of Einstein gravity with perfect-fluid sources [3,4]. The methods used in these works to eliminate the constraint equations, put the equations in the form (1), and obtain an inner product seemed rather particular to the theory he was examining, and it was far from certain that they would generalize to an arbitrary field theory. In a recent paper [1], we showed that these methods did in a certain sense "have to work out", by describing a straightforward procedure by which the gauge could be fixed, the constraints could be solved, and an inner product could be obtained under which the resulting time-evolution operator T is self-adjoint. We review this procedure in Section II of the paper. We then use this "generalized variational principle" to analyze three alternative theories of gravity which have garnered some attention in recent years: f (R) gravity, the current interest in which was primarily inspired by the work of Carroll et al. [5]; Einstein-aether theory, a toy model of Lorentz-symmetry breaking proposed by Jacobson and Mattingly [6,7]; and TeVeS, a covariant theory of MOND proposed by Bekenstein [8]. We apply our techniques to these theories in Sections III, IV, and V, respectively.
We will use the sign conventions of [9] throughout. Units will be those in which c = G = 1.

A. Symplectic dynamics
We first introduce some necessary concepts and notation. Consider a covariant field theory with an action of the form where L[Ψ] is a scalar depending on some set Ψ of dynamical tensor fields including the spacetime metric. (For convenience, we will describe the gravitational degrees of freedom using the inverse metric g ab rather than the metric g ab itself.) To obtain the equations of motion for the dynamical fields, we take the variation of the four-form Lǫ with respect to the dynamical fields Ψ: where a sum over all fields comprising Ψ is implicit in the first term. Requiring that δS = 0 under this variation then implies that the quantities E Ψ vanish for each dynamical field Ψ. The second term in (5) defines the vector field θ a [Ψ, δΨ]. The three-form θ dual to this vector field (i.e., θ bcd = θ a ǫ abcd ) is the "symplectic potential current." Taking the antisymmetrized second variation of this quantity, we then obtain the symplectic current three-form ω for the theory: In terms of the vector field ω a dual to ω, this can also be written as The symplectic form for the theory is then obtained by integrating the pullback of this three-form over a spacelike three-surface Σ: If we define n a as the future-directed timelike normal to Σ and e to be the induced volume three-form on Σ (i.e., e bcd = n a ǫ abcd ), this can be written in terms of ω a instead: (ω a n a )e. (9) In performing the calculations which follow, it is this second expression for Ω which will be most useful to us.

B. Obtaining a variational principle
In [1], we presented a procedure by which a variational principle for spherically symmetric perturbations of static, spherically symmetric spacetimes could generally be obtained. For our purposes, we will outline this method; further details can be found in the original paper.
The method described in [1] consists of the following steps: 1. Vary the action to obtain the equations of motion (E G ) ab corresponding to the variation of the metric, as well as any other equations of motion E A corresponding to the variations of any matter fields present. This variation will also yield the dual θ a of the symplectic current potential; take the antisymmetrized variation of this quantity (as in (7)) to obtain the symplectic form (9).
2. Fix the gauge for the metric, and choose an appropriate set of spacetime functions to describe the matter fields. Throughout this paper, we will choose our coordinates such that the metric takes the form for some functions Φ and Λ. 1 Our spacetimes will be static at zero-order (i.e., (∂/∂t) a is a Killing vector field at zero-order), but non-static in first-order perturbations. We will therefore have Φ(r, t) = Φ(r)+φ(r, t), where φ is a first-order quantity; similarly, we define the first-order quantity λ such that Λ(r, t) = Λ(r) + λ(r, t). In other words, and Similarly, all matter fields will be static in the background, but possibly time-dependent at first order.
3. Write the linearized equations of motion and the symplectic form in terms of these perturbational fields (metric and matter.) 4. Solve the linearized constraints. One of the main results of [1] was to show that this can be done quite generally for spherically symmetric perturbations off of a static, spherically symmetric background. Specifically, suppose the field content Ψ of our theory consists of the inverse metric g ab and a single tensor field A a1...an b1...bm . (The generalization to multiple tensor fields is straightforward.) Let (E A ) a1...an b1...bm denote the equation of motion associated with A a1...an b1...bm . We define the constraint tensor as  (13) where the summations run over all possible index "slots", from 1 to n and from 1 to m for the first and second summation respectively. It can then be shown that if the background equations of motion hold, and the matter equation of motion also holds to first order, the perturbations of the tensor C ab will satisfy and ∂F ∂r = −r 2 e Λ−Φ δC tt (15) for some quantity F which is linear in the firstorder fields. Moreover, the first-order constraint equations δC tt = δC tr = 0 will be satisfied if and only if F = 0. We can then solve this equation algebraically for one of our perturbational fields, usually the metric perturbation λ. We will refer to the equation F = 0 as the preconstraint equation.

5.
Eliminate the metric perturbation φ from the equations. As φ cannot appear without a radial derivative (due to residual gauge freedom) 2 , we must find an algebraic equation for ∂φ/∂r. The first-order equation δC rr = 0, will, in general, serve this purpose [1]. Use the above relations for λ and ∂φ/∂r to eliminate the metric perturbations completely from the perturbational equations of motion and the symplectic form, leaving a "reduced" set of equations of motion and a "reduced" symplectic form solely in terms of the matter variables.
6. Determine if the reduced equations of motion take the form (1). If so, read off the time-evolution operator T .
7. Determine if the symplectic form, written in terms of the reduced dynamical variables ψ α , is of the form for some three-form W αβ . Then (as we showed in [1]) we must have W αβ = W βα and 8. Determine whether W αβ is positive definite, in the sense that for all ψ α , with equality holding only when ψ α = 0. 3 If this is the case, we can define an inner product on the space of all reduced fields: Equation (17) then shows that T is a symmetric operator under this inner product. Thus, we can write down our variational principle of the form (3).
It is important to note that each step of this procedure is clearly delineated. While the procedure can fail at certain steps (the reduced equations of motion can fail to be of the form (1), for example, or W αβ can fail to be positive definite), there is not any "art" required to apply this procedure to an arbitrary covariant field theory. (In practice, as we shall see, there are certain shortcuts that may arise which we can exploit; however, the "long way" we have described here will still work.) In the following three sections, we will use this formalism to analyze the stability of f (R) gravity, Einstein-aether theory, and TeVeS.

A. Theory
In f (R) gravity, the Ricci scalar R in the Einstein-Hilbert action is replaced by an arbitrary function of R, leaving the rest of the action unchanged; in other words, the Lagrangian four-form L is of the form where A denotes the collection of matter fields, with tensor indices suppressed. Taking the variation of this action with respect to the metric, we obtain the equation of motion (21) where T ab , given by is the matter stress-energy tensor, and f ′ (R) = df /dR. This equation is fourth-order in the metric, and as such is somewhat difficult to deal with. We can reduce this fourth-order equation to two second-order equations using an equivalent scalar-tensor theory [11,12]. 4 This equivalent theory contains two dynamical gravitational variables, the inverse metric g ab and a scalar field α, in addition to the matter fields. The Lagrangian is given by (23) Varying the gravitational part of action with respect to and the vector θ a is our symplectic potential current: The vector θ a mat above is the symplectic potential current resulting from variation of L mat , and θ a Ein is the symplectic potential current for pure Einstein gravity, i.e., The equations of motion are then given by (E G ) ab = 0 and E α = 0. Assuming that f ′′ (α) = 0, this second equation implies that R = α, and substituting this relation into (25) yields the equation of motion obtained in (21).
Hereafter, we will use the form of the equations obtained from the action (23).

B. Obtaining a variational principle
Before we examine the stability of spherically symmetric static solutions in f (R) gravity with perfect fluid matter, we must first consider the question of whether physically realistic solutions exist. In particular, do there exist solutions to the equations of motion which reproduce Newtonian gravity, up to small relativistic corrections? There has recently been a good deal of debate on this subject [14,15,16,17,18]. We will not comment directly on this controversy here except to say that if f (R) gravity (or any theory) does not allow interior solutions with R ≈ 8πGρ to be matched to exterior solutions with R close to zero, then it is difficult to see how such a theory could reproduce Newtonian dynamics in a nearly-flat spacetime. We will therefore give the theory the "benefit of the doubt," and assume that such solutions exist.
An initial attempt to address the stability of such solutions was made by Dolgov and Kawasaki [19]; their perturbation analysis implied that stars would be extremely unstable in Carroll et al.'s f (R) gravity, with a characteristic time scale of approximately 10 −26 seconds. Their results, while suggestive, nevertheless failed to take into account the constrained nature of the theory: the stressenergy tensor, the metric, and the scalar α cannot all be varied independently. The imposition of a constraint can, of course, change whether or not a given system is stable; as a trivial example, consider a particle moving in the potential V (x, y) = x 2 − y 2 , with and without the constraint y = const. 5 In what follows, we will show that Dolgov and Kawasaki's conclusion is, nonetheless, correct: stars in CDTT f (R) gravity do in fact have an ultra-short timescale instability.
To describe the fluid matter, we will use the "Lagrangian coordinate" formalism, as in Section V of [1]. In this formalism, the fluid is described by considering the manifold M of all fluid worldlines in the spacetime, equipped with a volume three-form N. If we introduce three "fluid coordinates" X A on M, with A running from one to three, then the motion of the fluid in our spacetime manifold M is completely described by a map χ : M → M associating with every spacetime event x the fluid worldline X A (x) passing through it. The matter Lagrangian is then given by where ν, the "number density" of the fluid, is given by In turn, N abc , the "number current" of the fluid, is given by We will be purely concerned with spherically symmetric solutions and radial perturbations of the fluid; thus, we will take our Lagrangian coordinates to be of the form in the background, and consider only perturbations δX R ≡ ξ(r, t) at first order (i.e., δX Θ = δX Φ = 0.) We can then easily obtain the background equations of motion; these are we can rewrite the Lagrangian in terms of a rescaled metric g ab = Ω −2 (σ)g ab = e σ/2 √ 3π g ab , putting it in the form [11,12] where σ is related to the scalar σ by f ′ (α) = e σ/2 √ 3π and the exact form of the potential V (σ) is determined by the function f (R). The decoupling of the scalar perturbations, however, requires V (σ 0 ) = 0 and Ω ′ (σ 0 )/Ω = 0, where σ 0 is the background value of σ. The former condition will not hold for a generic f (R), and the latter condition will not hold for any f (R) for which f ′′ (R) = 0. and where ̺ ′ = d̺/dν. These equations are not all independent; in particular, the Bianchi identity implies that (34c) is automatically satisfied if the other four equations are satisfied. Note that under the substitutions ̺ → ρ, ̺ ′ ν − ̺ → P , the "matter terms" in these equations take on their familiar forms for a perfect fluid. We next obtain the symplectic form for the theory. The contribution θ a grav to the symplectic potential current from the gravitational portion of the action is given by (27); to obtain the symplectic form, we then take the antisymmetrized second variation of θ, as in (7). Performing this variation, we find that Ein is defined, analogously to θ a Ein , to be the symplectic current associated with pure Einstein gravity. This is equal to [21] The contribution to the symplectic current from the matter terms in the Lagrangian coordinate formalism was calculated in [1]; we simply cite the result for the t-component of ω a in a static background here: where the antisymmetrization in the second term is over the tensor indices only (not the fluid-space indices). Writing out ω t in terms of our perturbational variables, we have (Note that the t-component of the symplectic current about a spherically symmetric static solution vanishes in pure Einstein gravity, i.e., ω t Ein = 0.) The first main step towards obtaining a variational principle is to solve the linearized constraints. To do this, we must calculate the constraint tensor C ab , as defined in (13). Since all the fields in our theory other than the metric are scalars with respect to the spacetime metric, the resulting expression is particularly simple: where (E G ) ab , the metric equation of motion, is given by (25). We can then find algebraic equations for λ and ∂φ/∂r. The algebraic equation for λ will be given by the solution of the equation F = 0, where F is given by (14) and (15). Similarly, an algebraic equation for ∂φ/∂r can be found by solving the equation The situation is thus very similar to the example of pure Einstein gravity coupled to perfect fluid matter described in [1]; only the precise form of the perturbational equations of motion F = 0 and (δE G ) rr = 0 are different. In the f (R) gravity case, these equations are and where we have introduced the new quantity for notational convenience. Equation (42) then implies that our algebraic equation for λ is We could plug this result into (43) to obtain an algebraic equation for ∂φ/∂r; however, the resulting expression is somewhat complicated. In fact, we can derive a simpler expression for ∂φ/∂r. To do so, we note that In the case of f (R) gravity, this equation is given by Using the trace of this equation to eliminate the f ′ (α) term, along with the equation R = α, we can show that To zero-order, the θθ-component of this equation is and to first order, it is This equation, with λ given by (45), can then be solved algebraically for ∂φ/∂r in terms of the matter variables and their derivatives. Our next step is to obtain the reduced matter equations of motion, i.e., eliminate the metric degrees of freedom λ and ∂φ/∂r from the remaining first-order equations of motion. These remaining equations of motion are the "scalar" equation R − α = 0, which at first order is and the matter equation of motion, which as in the case of pure Einstein gravity is We can then eliminate the gravitational equations of motion using (45) and (49), leaving equations solely in terms of the matter variables b and ξ. Equivalently, we can write our "reduced" equations in terms of b and a new variable ζ, defined by (The utility of this new variable will become evident when we obtain the reduced symplectic form.) Performing these algebraic manipulations, the resulting reduced evolution equation for ζ takes the form (53) where the coefficients A i are dependent on the background fields. Similarly, the evolution equation for b takes the form We see from the form of the above equations that f (R) gravity has cleared another hurdle required to have a valid variational principle: the reduced equations do in fact take the form (1), containing only second derivatives in time and only up to second radial derivatives. We can also apply these constraint equations to the symplectic form. For the symplectic form, we only need to eliminate λ using the first-order constraint equation (45). Performing an integration by parts to eliminate the mixed derivatives that result, and applying the background equations of motion, we obtain our reduced symplectic form: As noted above (16), this reduced symplectic form defines a three-form W αβ . For a valid variational principle to exist, this W αβ must be positive definite in the sense of (18). We can see that the W αβ defined by (55) is positive definite if and only if f ′ (α) > 0 and ̺ ′ ν > 0 in our background solutions. The latter condition will hold for any matter satisfying the null energy condition, since ̺ ′ ν = ρ + P ; however, the former condition must be checked in order to determine whether a valid variational principle exists. 6 For the particular f (R) chosen by Carroll et al. [5], we have f ′ (R) = 1 + µ 4 /R 2 > 0, and so W αβ is always positive definite in the required sense for this choice of f (R). All that remains is to actually write down the variational principle for f (R) gravity. As noted above (3), the variational principle will take the form where T is the time-evolution operator. In our case, ψ denotes a two-component vector, whose components are the functions ζ and b. The inner product in which T is self-adjoint can be read off from (55), and thus the denominator of (56) will be (57) The numerator of (56) is, of course, rather more complicated. After multiple integrations by parts and applications of the background equations of motion, we can put this quantity in the form where the C i coefficients are functions of the background fields. The exact forms of the C i are given in Appendix A.

C. Discussion
All of our results thus far have been independent of the choice of f (R) (assuming, of course, that f ′′ (R) = 0.) In the case of Carroll et al.'s f (R) gravity, we can now use this variational principle to show that this theory is highly unstable for Newtonian solutions. For the choice f (R) = R − µ 4 /R, we have f ′ (R) = 1 + µ 4 /R 2 and f ′′ (R) = −2µ 4 /R 3 . Moreover, for a quasi-Newtonian stellar interior, we will have α = R ≈ ρ, where ρ is the matter density in the star. In particular, this implies that since the above term will dominate over all the others in C 4 . (Note that for a star with the density of the Sun, and the choice of µ ≈ 10 −27 m −1 made by Carroll et al., ρ/µ 2 ≈ 10 10 .) This implies that for a test function ψ with ζ set to zero, we will have (note that f ′ (α) = 1 + µ 4 /ρ 2 ≈ 1 in the stellar interior.) As a representative mass distribution, we take the Newtonian mass profile of an n = 1 polytrope: where R is the radius of the star and ρ 0 is its central density. We take ρ 0 to be of a typical stellar density, ρ 0 ≈ 10 −24 metres −2 . Numerically integrating (59) with a test function of the form we find that which corresponds to an instability timescale of τ ≈ 4 × 10 −27 seconds. This timescale is of the same magnitude as the instability found in [19]. We see, then, that for the f (R) originally chosen by Carroll et al., the theory is extremely unstable in the presence of matter. Moreover, a similar argument will obtain for any choice of f (R) for which quasi-Newtonian solutions exist and for which f ′′ (R) is sufficiently small and negative at stellar-density curvature scales. We can always pick a test function b(r) lying entirely inside the stellar interior such that ∂b/∂r is of order b/R. Thus, if |f ′ (ρ)/f ′′ (ρ)| ≫ 1/R 2 for a typical stellar density ρ, the b 2 term in (59) will dominate the (∂b/∂r) 2 term, and the resulting lower bound on ω 2 0 will then be of order f ′ (α)/f ′′ (α). If the choice of f (α) results in this quantity being negative, quasi-Newtonian stellar solutions will be unstable in the corresponding theory. We can thus rule out any theory (e.g., [5,22]) for which f ′′ (ρ) is sufficiently small and negative.
This result seems to be closely related, if not identical to, the "high-curvature" instabilities found in cosmological solutions by Song, Hu, and Sawicki [23]. Particularly telling is the fact that the instability time scale found in their work is proportional to |f ′′ (R)/f ′ (R)|, the same time-scale found in the present work.

A. Theory
Einstein-aether theory [6,7] was first formulated as a toy model of a gravitational theory in which Lorentz symmetry is dynamically broken. This theory contains, along with the metric g ab , a vector field u a which is constrained (via a Lagrange multiplier Q) to be unit and timelike. The Lagrangian four-form for this theory is where A denotes any matter fields present in the theory, and The c i constants determine the strength of the vector field's coupling to gravity, as well as its dynamics. 7 (Note 7 Note that our definitions of the coefficients c i differ from those in [6,7,24,25], as do the respective metric signature conventions. The net result is that to translate between our results and those of the above papers, one must flip the signs of all four coefficients. that in the case c 3 = −c 1 > 0 and c 2 = c 4 = 0, we have the conventional kinetic term for a Maxwell field; this special case was examined in [26], prior to Jacobson and Mattingly's work.) In the present work, we will work in the "vacuum theory", i.e., in the absence of matter fields A.
Performing the variation of the Lagrangian four-form, we find that and In the above, we have introduced the notationu a = u b ∇ b u a and J a b = K ac bd ∇ c u d ; θ a Ein is defined by (28), as above. Equation (67c) is the constraint that u a is unit and timelike, while (67a) and (67b) are the equations of motion for g ab and u a , respectively. If desired, we can eliminate the Lagrange multiplier Q from these equations by contracting (67b) with u a , resulting in the equation We now take the variation of the symplectic potential current to obtain the symplectic current. The resulting expression can be written in the form where ω a Ein is the usual symplectic current for pure Einstein gravity (given by (36)), and ω a vec is obtained by taking the antisymmetrized variation of the last two terms in (68). The exact form of this tensor expression is rather complicated, and can be found in Appendix B.

B. Obtaining a variational principle
As before, we are primarily concerned with the tcomponent of ω a , and its form in terms of the perturbational variables. We will use the usual spherical gauge (10) for g ab . We will further assume that u a ∝ t a in the background solution (the so-called "static aether" assumption), and that to first order we have i.e., the first-order perturbation to u t is −φe −Φ while the first-order perturbation to u r is υ. (Note that under these conventions, g ab u a u b = −1 to first order, as required.) The background equations of motion can then be calculated to be and where we have defined c 14 = c 1 +c 4 and used the equation to eliminate Q. The final equation of motion, (E u ) a = 0, is satisfied trivially if the aether is static. In terms of these variables, the t-component of ω a can be calculated to be where we have defined c 123 = c 1 +c 2 +c 3 . In what follows, we will assume that c 14 = 0, and, except where otherwise noted, that c 123 = 0 as well. Our next step, as usual, is to solve the constraints. However, in this theory we have the added complication of the presence of a vector field. This means, in particular, that the tensor C ab (as defined in (13)) is not merely proportional to (E G ) ab , as in the previous section, but is instead Writing out δC rt in terms of our perturbational quantities, we find that it is indeed a total time derivative, While we could follow the methods outlined in Section II to reduce these equations to the basic form (1), it is actually simpler to pursue a different path. If we solve (78) for ∂φ/∂r, rather than λ as usual, and plug the resulting expressions into (79) and the reduced symplectic form can be written as The form W αβ thus defined may be positive or negative definite, depending on the signs of c 123 and c 14 (recall that we are assuming that c 123 and c 14 are nonvanishing); however, it is never indefinite. Thus, the symplectic form is in fact of the required form (16), and the equations of motion can be put in the form (1). We can therefore write down our variational principle of the form (3); the denominator is and the numerator is We can now examine the properties of this variational principle to determine the stability of Einstein-aether theory. The simplest case to examine is that of flat space. Let us denote the coefficient of ψ 2 in (90) by Z(r), i.e., In the case of flat spacetime, Z(r) vanishes, and we are left with (92) Since both integrands are strictly positive, we conclude that flat space is stable to spherically symmetric perturbations in Einstein-aether theory if and only if This is the same condition found in [7] for the stability of the "trace" wave mode (modulo the aforementioned sign differences). The other stability conditions found in [7] those corresponding to the "transverse traceless metric" and "transverse aether" wave modes -are incompatible with spherical symmetry, and thus our analysis cannot reproduce these conditions. More broadly, we can also analyze the stability of the general spherically symmetric solutions described in [24]. These (exact) solutions are described in terms of a parameter Y : where the constants Y ± are given by and r min is a constant of integration related to the mass M via We can then obtain Z(Y ) by writing out (91) in terms of these functions of Y (noting that, for example, ∂Λ/∂r = (∂Λ/∂Y )/(∂r/∂Y )), and then plot Z and r parametrically. The resulting function is shown in Figure 1, for c 123 = ±c 14 and c 14 = −0.1. The first thing we see is that in the asymptotic region, the sign of Z(r) is determined by the sign of c 123 . To quantify this, we can obtain an asymptotic expansion of Z(r) as r → ∞, noting that r → ∞ as Y → 0. Doing this, we find that to leading order in M/r, Thus, the coefficient of the ψ 2 term in (90) is, to leading order in the asymptotic region, the same sign as that of the (∂ψ/∂r) 2 term. We can therefore conclude that for a test function in the "asymptotic region" of one of these solutions (i.e., r ≫ M ), the conditions on the c i 's are the same as those in flat spacetime. 8 For a realistic spacetime, of course, there will be some ball of matter in the central region, and the region where the vacuum solution holds will be precisely the region where r ≫ M ; thus, we can conclude that for a normal star in Einsteinaether theory, the exterior is stable if and only if (93) holds. (Henceforth, we will assume that the c i coefficients have been chosen with this constraint in mind, unless otherwise specified.) We also note that as r approaches r min , Z(r) diverges negatively. One might ask whether a test function in this region (in a spherically symmetric spacetime surrounding a compact object, say) could lead to an instability. We investigated this question numerically using our variational principle; however, our results for (ψ, T ψ) were positive for all test functions ψ that we tried. Roughly speaking, the derivative terms in (89) always won out over the effects of the negative Z(r). 9 This is, of course, far from a definitive proof of the positivity of (90); and it should be emphasized that the above analysis has not considered the effects of matter on the stability of such solutions. Nevertheless, the above results are at least indicative that the spherically symmetric vacuum solutions Einstein-aether theory do not possess any serious stability problems. 9 We have been assuming a "static aether" throughout this work; however, there also exist solutions in which the aether vector u a is not globally aligned with the time-translation vector field [25]. Thus, even supposing that such "extremely compact stars" are phenomenologically realistic, and that the exterior solutions for such objects are unstable, these "non-static aether" solutions might still be physically viable.
D. The case of c123 = 0 The above analysis assumed that c 123 was nonvanishing; however, the action originally considered by Jacobson and Mattingly [6] used a "Maxwellian" action, i.e., the kinetic terms for the aether field were of the form F ab F ab , corresponding to c 3 = −c 1 and c 2 = 0. This case is therefore of some interest.
In the preceding analysis, the perturbational equations of motion (79) and (80) were obtained without any assumptions concerning c 123 , as was the preconstraint equation (78). We can therefore simply set c 123 to zero in these equations and use (78) to solve for ∂φ/∂r, as we did in the general case. After application of the background equations of motion, there result the equations and where we have defined Z 0 in terms of the background fields: We can then combine these two equations to make a single second-order equation for υ: Since there are no radial derivatives to contend with, the solutions to this equation are all of the form (104) where f (r) is an arbitrary function of r.
Stability of these solutions will thus depend on the sign of the quantity ( 2 c14 + 1)Z 0 (r). We can plot Z 0 (r) parametrically, as was done in the general case; the resulting function, shown in Figure 2, falls off rapidly in the asymptotic region. We can also perform an asymptotic expansion similar to that done in the general case to find the large-r behaviour of Z 0 (r); the result is that We conclude that in the c 123 = 0 case, the spherically symmetric static solutions of Einstein-aether theory are unstable unless i.e., unless either c 2 > 1 and c 14 > −2, or c 2 < 1 and c 14 < −2. Note that this excludes the limit of "small c i " often considered in works such as [7]. In the "Maxwellian" case, this instability is likely related to the non-boundedness of the Hamiltonian [27]; however, our analysis also suggests that non-standard kinetic terms can in fact stabilize the theory (at least in the timelike aether case), contrary to the arguments made in that work.
Finally, we note that even if (106) holds, the perturbational solutions of the form (104) may still be problematic, as the radial gradients of υ will grow linearly with time. To quantify what magnitude of gradients would be acceptable, we note that in the c 123 = 0 case, there is only one dynamical degree of freedom; thus, υ uniquely determines λ and ∂φ/∂r. Setting c 123 = 0 in (78) and (79) allows us to solve these equations for λ and ∂φ/∂r; in particular, which simplifies in a nearly-flat spacetime to be Applying this to our solution (103), we find that for sufficiently large t We can estimate typical scales for ∂φ/∂r by looking at other sources, such as planets. The scale of perturbations due to the Earth's gravitational field (at its surface) is given by ∂φ/∂r ≈ M ⊕ /r 2 ⊕ , where M ⊕ and r ⊕ are the mass and radius of the Earth, respectively. Thus, the perturbation to ∂ 2 φ/∂r 2 near the Earth's surface will be of the order where R is the radius of Earth's orbit and M ⊙ is the mass of the Sun. This enhancement to ∂ 2 φ/∂r 2 will lead to an observable change in the tidal effects due to the Sun's gravity; demanding that these remain small relative to the normal Newtonian tidal effects, i.e., ∂ 2 φ/∂r 2 ≪ ∂ 2 Φ/∂r 2 ≈ 2M ⊙ /R 3 , we then have for t ≈ 5 × 10 9 years (the approximate age of the Earth.) We see, therefore, that the requirement that perturbational tidal effects remain small severely constrains our choices of c 14 and c 2 .

A. Theory
TeVeS (short for "Tensor-Vector-Scalar") is a modified gravity theory proposed by Bekenstein [8] in an attempt to create a fully covariant theory of Milgrom's Modified Newtonian Dynamics (MOND). The fields present in this theory consist of the metric; a vector field u a which (as in Einstein-aether theory) is constrained by a Lagrange multiplier Q to be unit and timelike; and two scalar fields, α and σ. The Lagrangian four-form is where L g is the usual Einstein-Hilbert action, R; (113) L s is the "scalar part" of the action, with k and ℓ positive constants of the theory (with "length dimensions" zero and one, respectively), and F (x) a free function; L v is the vector part of the action, with K a positive dimensionless constant and F ab = ∇ a u b − ∇ b u a ; and L m the matter action, non-minimally coupled to the metric: where L mat [A, g ab ] would be the minimally coupled matter Lagrangian for the matter fields A. Note that if we ignore the scalar and matter portions of the Lagrangian, this Lagrangian is the same as that for Einstein-aether theory (64), with c 1 = −c 3 = −K/16π and c 2 = c 4 = 0.
In the present work, we will work exclusively with the "vacuum" (L m = 0) theory. Taking the variation of (112) to obtain the equations of motion and the symplectic current, we find that and θ a = θ a Ein − σ 2 (∇ a α − u aα )δα In the above, we have definedα ≡ u a ∇ a α. The remaining equation E Q = 0 is identical to that in Einstein-aether theory, (67c). For a static solution, in our usual gauge, the background equations of motion become and where we have defined y(x) = −xF (x) − 1 2 x 2 F ′ (x) (as in [8]) and used the equation for Q, to simplify. (As in Einstein-aether theory, the background equation (E u ) a = 0 is satisfied trivially by a static aether.) The symplectic form for the theory can now be calculated from (119). The result will be essentially the same as that in Einstein-aether theory (with the appropriate values of the c i 's), with added terms stemming from the variations of L s : where ω a Ein is given by (36), ω a vec is given by (B1) with c 1 = −c 3 = −K/16π and c 2 = c 4 = 0, and

B. Applying the variational formalism
The next step is to write out the t-component of ω a in terms of the perturbational variables. We will take our metric perturbations to have the usual form, and our vector perturbation to be of the same form as was used for Einstein-aether theory (71). For the two scalar fields, we define δα ≡ β and δσ ≡ τ . Calculating ω t s in terms of these perturbational variables, and using the results of (75), we find that We now turn to the question of the constraints. The tensor C ab is again given by C ab = (E G ) ab + u a (E u ) b , and thus to first order we have δC rt = (E G ) rt + δu r (E u ) t (since u r = 0 in the background.) Calculating δC rt in terms of our perturbational variables, we find that the preconstraint equation is In principle, we could now use this equation, together with the equation to obtain equations for λ and ∂φ/∂r in terms of the "matter" variables υ, β, and τ . However, it is simpler to pursue a similar tactic to the one we used in the reduction of Einstein-aether theory. To wit, the r-component of (δE u ) a is We can use (125) to simplify this equation; the result is We then define the new variable χ as The evolution equation for β, meanwhile, is given by the equation δE α = 0; in terms of the perturbational fields, it is We can thus follow the following procedure to write the evolution equations solely in terms of β and χ: first, we use the preconstraint equation (125) to eliminate the combination ∂φ/∂r + e 2Λ−Φ ∂υ/∂t in favour of λ and β; next, we use (131) to eliminate τ and its spatial derivatives; and finally, we use the definition of χ (129) to eliminate λ and its derivatives in favour of β and χ. Applying our procedure to (126) results in an equation depending on ∂υ/∂t as well as β, ∂β/∂r, and χ; combining this equation with (128) then yields a second-order evolution equation for χ of the form where the coefficients V i and W are dependent on the background fields. Similarly, applying this procedure to (130), and using (126) to eliminate a resulting ∂υ/∂t term, we obtain a second-order equation for β: where, again, the U i 's and V i 's depend on the background fields. Note that (132) has no dependence on the spatial derivatives of χ, only upon χ itself.
We also need to reduce the symplectic form and express it in terms of β and χ. Applying the preconstraint equation (125) to (124), and using the definition of χ (129) along with (128), we find that where We can also apply the background equations of motion, along with the Equation (121), to this quantity to obtain For a variational principle to exist for this theory, the form W αβ defined by (134) must be positive definite. In our case, this means that the coefficients of both (∂χ 1 /∂t)χ 2 and (∂β 1 /∂t)β 2 in the integrand of (134) must be positive for the background solution about which we are perturbing. While the coefficient for the latter term is obviously always positive, the situation for the former coefficient (namely, H) is not so clear. To address this issue, we need to know the properties of the spherically symmetric static background solutions of TeVeS. These solutions (with a "static aether") are described in [28]. In our gauge, they are most simply described in which is negative for any positive choice of k and K. Therefore, it is not possible to straightforwardly apply the variational principle to TeVeS, since the quadratic form defined by W αβ is indefinite and thus cannot used as an inner product. 11

C. WKB analysis
While we cannot derive a variational principle with TeVeS, the fact that we have been able to "reduce" the equations of motion to an unconstrained form still allows us to analyse the stability of its spherically symmetric solutions. Let us consider a WKB ansatz, of the form with κ very large compared to the scale of variation of the background functions U i , V i , and W. We will further choose the the functions f i (r) and ω(r) are chosen to be "slowly varying" relative to the scale defined by κ, i.e.
Under this assumption, we will then have Now let us apply the time-evolution operator T implicitly defined by (133) and (132) to our ansatz (143). We see that for sufficiently large κ, the highest-derivative terms will dominate the lower-derivative terms. Thus, to a good approximation we will have (146) Then, in the limit of large κ, our ansatz (143) will be an approximate eigenvector of T if there exist a f β (r) and f χ (r) such that In other words, in the limit of large κ, the problem of finding modes of T is a simple two-dimensional eigenvalue problem where the eigenvalues are functions of r. In this limit, the eigenvalues of this matrix are (to leading order in κ) ω 2 (r) ≈ κ 2 U 1 (r), −W(r) + V 1 (r)V 3 (r) U 1 (r) .
It is not difficult to verify that for sufficiently large κ, our assumptions for the ansatz (144) are satisfied. It remains to write out the functions U 1 (r), V 1 (r), V 3 (r), and W(r) in terms of the background functions. These can be shown to be: We can see that this first eigenvalue (153) is always positive as long as y ′ (x) > 0; indeed, the choice of y(x) made in [8] does satisfy this inequality. The second eigenvalue (154), however, is just Since H(r) is always negative, we conclude that this second mode is unstable for 0 < K < 2. Further, for a spherically symmetric solution outside a Newtonian star, the approximation (142) is valid; thus, to lowest nonvanishing order in r −1 , we have

D. Discussion
The above result for ω 2 (r) implies that for stability of Newtonian solutions, we must have K > 2, ruling out the range of parameters originally considered in [8]. To estimate the time scale of the instability when K < 2, we first note that unless k and K are much larger than the ratio of the star's radius R to its Schwarzchild radius m g , the parameters m g , m s , and z g are all approximately equal, differing by terms of O(m g /R) [8]. In particular, if k and K are of the same order and K ≈ 2, we find that the timescale of this instability will be on the order of 10 6 seconds-approximately two weeks-for points near the surface of the Sun.
Note also that the effective gravitational potential for nonrelativistic motion is where Ξ = e −2αc (1 − K/2) −1 , Φ N is the Newtonian gravitational potential, and α c is the asymptotic value of α, determined by the cosmological boundary conditions [29]. Expanding the asymptotic behaviour of Φ N and α about infinity, we find that where in the weak-field limit, we have and m s ≈ e −2αc m g . The coefficient in parentheses in (159) can be thought of as renormalizing the gravitational constant G by some overall factor; the first term in parentheses comes from the usual perturbations to the "physical metric" componentg tt , while the second term comes from variation of the scalar field α. However, if K > 2, this would make the metric contribution to the "effective gravitational constant" negative. This could be remedied by choosing k > 8π/(K − 2). However, it is unclear whether the theory, in this parameter regime, would still be experimentally viable; for example, Big Bang nucleosynthesis constraints require that k < 0.75, and the CMB power spectrum also constrains the viable regions of parameter space [30]. 12 We also note that the eigenvector corresponding to the unstable mode of T will satisfy (to leading order) or, in the limit of large κ, f β ≈ 0. In other words, this unstable mode should manifest itself in growth of the radial component of the vector field u a rather than growth of the scalar α (cf. (128)). An instability of the vector field was also found in [32] by Dodelson and Liguori. However, it seems unlikely that this is the same instability for three reasons. First, the instability found in [32] was found in a cosmological context, not a Newtonian-gravity context; since the cosmological solutions of TeVeS are in a very real sense separate from the Newtonian solutions (existing on two different branches of the function y(x)), it is difficult to draw a direct correspondence between the stability properties of these two types of solutions. Second, the instability found in [32] manifests itself only in the limit of a matter-dominated Universe; the instability we have found exists in vacuo. Third, Dodelson and Liguori's instability requires K to be sufficiently small relative to k, while our instability is present for all k > 0 and 0 < K < 2.
Finally, it is perhaps notable that the vector field u a in TeVeS has "Maxwellian" kinetic terms, which we found in Section IV to be unstable in the context of Einsteinaether theory. It is possible, then, that this instability could be cured via a more general kinetic term for the vector field.