Spherically symmetric solutions in torsion bigravity

TL;DR

This paper explores spherically symmetric solutions in a novel torsion bigravity theory, revealing differences from general relativity and bimetric gravity, especially in high-compactness star models, with implications for gravitational physics.

Contribution

It introduces torsion bigravity as a geometric framework with unique properties, including a specific spectrum of excitations and distinct perturbative behavior compared to bimetric gravity.

Findings

Spherically symmetric solutions have the same degrees of freedom as ghost-free bimetric gravity.

No inverse powers of the spin-2 mass parameter appear at low orders in perturbation.

High-compactness star models differ significantly from GR predictions.

Abstract

We study spherically symmetric solutions in a four-parameter Einstein-Cartan-type class of theories. These theories include torsion, as well as the metric, as dynamical fields, and contain only two physical excitations (around flat spacetime): a massless spin-2 excitation and a massive spin-2 one (of mass $ m_2 \equiv \kappa$). They offer a geometric framework (which we propose to call "torsion bigravity") for a modification of Einstein's theory that has the same spectrum as bimetric gravity models. We find that the spherically symmetric solutions of torsion bigravity theories exhibit several remarkable features: (i) they have the same number of degrees of freedom as their analogs in ghost-free bimetric gravity theories ( i.e. one less than in ghost-full bimetric gravity theories); (ii) in the limit of small mass for the spin-2 field ($ \kappa \to 0$), no inverse powers of $\kappa$ arise at the first two orders of perturbation theory (contrary to what happens in bimetric gravity where $1/\kappa^2$ factors arise at linear order, and $1/\kappa^4$ ones at quadratic order). We numerically construct a high-compactness (asymptotically flat) star model in torsion bigravity and show that its geometrical and physical properties are significantly different from those of a general relativistic star having the same observable Keplerian mass.

Figures

Equations

$$S_{\mathrm{total}}=S_{\mathrm{TBG}}[{e^i}_{\mu },A_{ij\mu }]+S_{\mathrm{matter}}.$$

$$S_{\mathrm{TBG}}[{e^i}_{\mu },A_{ij\mu }]=\int d^{4}x\sqrt{g}{L}_{\mathrm{TBG}}[e,\partial e,{\partial }^{2}e,A,\partial A],$$

$$L_{\mathrm{TBG}}$$

$$c_R+c_F$$

$$\frac{c_F}{c_R}$$

$$c_{F^2}$$

$${K^i}_{j\mu }\equiv {A^i}_{j\mu }-{{\omega }^i}_{j\mu }(e).$$

$$K_{ijk}$$

$$T_{i[jk]}$$

$${T^{\lambda }}_{\mu \nu }\equiv {{\Gamma }^{\lambda }}_{\mu \nu }-{{\Gamma }^{\lambda }}_{\nu \mu }\equiv 2{{\Gamma }^{\lambda }}_{[\mu \nu ]}$$

$$\int d^{4}x{\Lambda }^{\lambda (\mu \nu )}{Q}_{\lambda ,(\mu \nu )}\equiv \int d^{4}x{\Lambda }^{\lambda \mu \nu }{\nabla }_{\lambda }^{\Gamma }{g}_{\mu \nu }$$

$$e^{2\Phi }$$

$$e^{2\Lambda }$$

$$d{s}^2=-e^{2\Phi }d{t}^2+e^{2\Lambda }d{r}^2+r^2\left(d{\theta }^2+{\sin }^2\theta d{\varphi }^2\right).$$

$$e^{\hat{0}}=e^{\Phi }dt,e^{\hat{1}}=e^{\Lambda }dr,e^{\hat{2}}=rd\theta ,e^{\hat{3}}=r\sin \theta d\varphi .$$

$$A_{\hat{a}\hat{0}\hat{0}}=\overline{V}(r)n_a,$$

$$A_{[\hat{a}\hat{b}]\hat{c}}=\overline{W}(r)\left(n_a{\delta }_{bc}-n_b{\delta }_{ac}\right),$$

$$V(r)$$

$$W(r)$$

$${A^{\hat{1}}}_{\hat{0}}$$

$${A^{\hat{1}}}_{\hat{2}}$$

$${A^{\hat{1}}}_{\hat{3}}$$

$${A^{\hat{2}}}_{\hat{3}}$$

$${{\omega }^{\hat{1}}}_{\hat{0}}$$

$${{\omega }^{\hat{1}}}_{\hat{2}}$$

$${{\omega }^{\hat{1}}}_{\hat{3}}$$

$${{\omega }^{\hat{2}}}_{\hat{3}}$$

$${K^{\hat{1}}}_{\hat{0}\hat{0}}$$

$${K^{\hat{1}}}_{\hat{2}\hat{2}}$$

$$T_{\hat{0}[\hat{1}\hat{0}]}$$

$$T_{\hat{2}[\hat{1}\hat{2}]}$$

$$G_{ij}(A)\equiv F_{ij}(A)-\frac{1}{2}{\eta }_{ij}F(A)$$

$$G_{\hat{t}\hat{t}}$$

$$G_{\hat{r}\hat{r}}$$

$$G_{\hat{\theta }\hat{\theta }}$$

$$\hat{t}=\hat{0},\hat{r}=\hat{1},\hat{\theta }=\hat{2},\hat{\varphi }=\hat{3}.$$

$$S=S_{\mathrm{field}}+S_{\mathrm{m}}.$$

$$\delta S_{\mathrm{m}}=\int \delta (\sqrt{g}{L}_m)d^{4}x=\frac{1}{2}\int \sqrt{g}{T}^{\mu \nu }\delta g_{\mu \nu }{d}^{4}x$$

Full text

I Introduction

Einstein’s theory of gravitation, i.e. General Relativity (GR), has, so far, been found to be in excellent accord with all gravitational observations and experiments. In particular, its foundational stone, the weak equivalence principle, has been recently confirmed at the $10^{-14}$ level Touboul:2017grn , while gravitational-wave observations have confirmed several basic dynamical predictions of GR LIGOScientific:2019fpa ; Abbott:2018lct . [See, e.g., chapter 20 in Ref. Tanabashi:2018oca for a review of the experimental tests of GR.]

However, since the discovery of GR more than a century ago, the quest for possible extensions of GR has been going on. We shall not discuss here the various motivations underlying the study of modified theories of gravity (see Ref. Capozziello:2011et for a review). Let us only mention that, from a pragmatic point of view, it is useful to have alternative theories of gravity to conceive and interpret tests of gravity Berti:2015itd .

Here, we study a class of geometric theories of gravitation that generalize the Einstein-Cartan theory. The original idea of Cartan Cartan:1923zea ; Cartan:1924yea ; Cartan1925 was to extend GR by considering the metric and the (affine) connection as a priori independent fields (first-order formalism), and by allowing the connection111Cartan worked within a vielbein formalism, in which the affine connection is naturally restricted to be metric-preserving; see below. to have nonzero torsion. Cartan added the idea that torsion might be sourced by some sort of intrinsic spin density along the matter worldlines222“En admettant la possibilité d’ éléments de matière doués de moments cinétiques non infiniment petits par rapport à leur quantité de mouvement.”; bottom of p. 328 in Cartan:1923zea . Later, Weyl pointed out that it is natural, in such a first-order formalism, to consider that fermions (Dirac spinors) directly couple to the connection, so that the torsion ${T^{i}}_{[jk]}=-{T^{i}}_{[kj]}$ is sourced by the microscopic (quantum) spin density of fermions $\sim\frac{1}{2}\bar{\psi}\gamma^{i}\gamma_{[j}\gamma_{k]}\psi$. [As explained in detail below, the latin indices $i,j,k,\cdots$ denote frame indices.] He also showed that if one follows Einstein and Cartan in using as gravitational action the first-order form of the scalar curvature, the torsion is algebraically determined by its source and that the first-order action is equivalent to a second-order (purely metric) action containing additional “contact terms” quadratic in the torsion source, and therefore quartic in fermions. The ideas of Cartan and Weyl were further developed by Sciama Sciama1962 , Kibble Kibble:1961ba and many others (see Hehl:1976kj for a review of later work on this approach based on gauging the Poincaré group).

A new twist in the story started after the discovery of supergravity Freedman:1976xh , and especially of its first-order formulation Deser:1976eh . Indeed, the first-order formulation of supergravity is similar to the Einstein-Cartan-Weyl approach, with a gravitational term linear in the scalar curvature, and a nonzero torsion ${T^{i}}_{[jk]}$ algebraically determined in terms of its gravitino source $\sim\bar{\psi}_{j}\gamma^{i}\psi_{k}$, leading, after replacement in the action, to contact terms quartic in the gravitino. However, quantum loops generate an effective action containing terms at least quadratic in the curvature. When considered in a purely metric, second-order formulation, terms quadratic in the curvature lead to higher-order field equations, which raise difficulties Stelle:1976gc , in the form of “ghosts” (negative-energy modes), even at the classical level Stelle:1977ry . This raised the issue of finding ghost-free theories of gravity with an action containing terms quadratic in curvature, but treated in a first-order formulation. Indeed, such a formulation leads to second-order-only field equations for the metric and the connection Neville:1978bk ; Neville:1979rb , so that the torsion now propagates away from its source. The most general solution to finding such ghost-free and tachyon-free (around Minkowski spacetime) theories with propagating torsion was obtained in parallel work by Sezgin and van Nieuwenhuizen Sezgin:1979zf ; Sezgin:1981xs , and by Hayashi and Shirafuji Hayashi:1979wj ; Hayashi:1980av ; Hayashi:1980ir ; Hayashi:1980qp . It was found that there are twelve six-parameter families of ghost-free and tachyon-free theories with propagating torsion Hayashi:1980qp ; Sezgin:1981xs . These theories always contain an Einsteinlike massless spin-2 field, together with some (generically) massive excitations coming from the torsion sector. The possible spin-parity labels of the excitations propagated by the torsion sector are: $2^{+}$, $2^{-}$, $1^{+}$, $1^{-}$, $0^{+}$, $0^{-}$. Only certain combinations of these spin-parities can be present in the various six-parameter families of ghost-free and tachyon-free theories with propagating torsion (see Table I in Hayashi:1980qp or Table I in Sezgin:1981xs ).

One of these classes of theories (with torsion propagating both massive $2^{+}$ and massive $0^{-}$ excitations) has recently been studied with the hope that the massive spin-2 field it contains will define a new, more geometric, solution to having a healthy and cosmologically relevant infrared modification of gravity Nair:2008yh ; Nikiforova:2009qr . We recall that the physics of an ordinary, massive333especially with a very small mass, say of cosmological scale. Fierz-Pauli-type Fierz:1939zz ; Pauli:1939xp ; Fierz:1939ix spin-2 field raises many subtle issues going by the names of: vanDam-Veltman-Zakharov discontinuity vanDam:1970vg ; Zakharov:1970cc , Vainshtein (conjectured) mechanism Vainshtein:1972sx , and Boulware-Deser ghost Boulware:1973my . A breakthrough in the problem of defining a class of consistent, ghost-free nonlinear theories of a massive spin-2 field was achieved in Ref. deRham:2010kj . This then allowed the construction of a class of consistent, ghost-free nonlinear theories of bimetric gravity Hassan:2011zd .

The aim of the present paper is to study the four-parameter subclass of the propagating-torsion models of Refs. Sezgin:1979zf ; Sezgin:1981xs ; Hayashi:1979wj ; Hayashi:1980av ; Hayashi:1980ir ; Hayashi:1980qp ; Nair:2008yh ; Nikiforova:2009qr that it similar to the bimetric gravity models of Hassan:2011zd in the sense that it contains only two types of excitations: an Einsteinlike massless spin-2 excitation, and a positive-parity massive spin-2 one. To emphasize this similarity we shall often refer to the models we study as defining a theory of torsion bigravity. We think that the geometric origin of the massive spin-2 additional field (contained among the torsion components, rather than through a second metric) makes such a torsion bigravity model an attractive alternative to the usually considered bimetric gravity models. In particular, the fact that massive gravity is described in these models by a different Young tableau than the more familiar (symmetric tensor) models completely changes the various issues linked to nonlinear effects, and renders the study of their physical properties a priori interesting. Some of the results of previous work on such models Nair:2008yh ; Nikiforova:2009qr ; Nikiforova:2016ngy ; Nikiforova:2017saf ; Nikiforova:2017xww ; Nikiforova:2018pdk has shown them to be remarkably healthy and robust around various backgrounds (though Ref. Nikiforova:2018pdk found the presence of gradient instabilities around the self-accelerating torsionfull cosmological solution found in Nikiforova:2016ngy ; but these instabilities might be due to the endemic stability problems of self-accelerating cosmological universes rather than to the theory itself). Anyway, let us emphasize here that the existence of the self-accelerating solution of Ref. Nikiforova:2016ngy necessarily relied on the presence in the spectrum of both $2^{+}$ and $0^{-}$ excitations. In the present work we focus on the minimal model containing only the $2^{+}$ excitation (besides the Einstein massless graviton). This minimal torsion bigravity model has not yet received any specific attention in the literature beyond its linearized approximation (which follows from the general linearized-limit results of Refs. Hayashi:1980ir ; Zhang:1982jn ; Nikiforova:2009qr ).

Let us note in passing, for the cognoscenti, that we are talking here about positive-parity spin-2 excitations contained in the torsion field ${T^{i}}_{[jk]}$, and not of the “dual gravity”, negative-parity spin-2 excitation contained in the irreducible SO$(3,1)$ Young tableau $T_{[ab]c}$ (satisfying $T_{[ab]c}+T_{[bc]a}+T_{[ca]b}=0$) introduced by Curtright Curtright:1980yk ; Curtright:1980yj . Among the propagating torsion models of Refs. Hayashi:1980qp ; Sezgin:1981xs some give rise to massive $2^{-}$ excitations and some to massive $2^{+}$ ones, but the two parities cannot be simultaneously present in ghost-free models.

As we started this Introduction by recalling that the source of torsion is the microscopic (quantum) spin of elementary fermions, the reader might worry that this would prevent the existence of phenomenologically relevant, macroscopic torsion fields in ordinary, non spin-polarized systems, such as stars, planets, or even neutron stars444We leave to future work a study of the amount of spin-polarization in a strongly magnetized neutron star.. However, as was already noticed in Refs. Hayashi:1980ir and Nikiforova:2009qr , and as will be clear in the present work, the mere presence of a usual, Einsteinlike energy-momentum tensor $T^{\mu\nu}$ suffices to generate macroscopic torsion fields. In the following, we shall then, for simplicity, set the torsion source to zero and only consider the effect of the energy-momentum source $T^{\mu\nu}$.

II Formalism and action of torsion bigravity

Here, we essentially follow the notation of Refs. Hayashi:1979wj ; Hayashi:1980av ; Hayashi:1980ir ; Hayashi:1980qp (which we also used in our previous paper Nikiforova:2018pdk ). Latin indices $i,j,k,\ldots=0,1,2,3$ (moved by the Minkowski metric $\eta_{ij},\eta^{ij}$) are used to denote Lorentz-frame indices referring to a vierbein ${e_{i}}^{\mu}$ (with inverse ${e^{i}}_{\mu}$) , while Greek indices $\mu,\nu,\ldots=0,1,2,3$ (moved by the metric $g_{\mu\nu}\equiv\eta_{ij}{e^{i}}_{\mu}{e^{j}}_{\nu}$) are used to denote spacetime indices linked to a coordinate system $x^{\mu}$. When there is a risk of confusion we add a hat, e.g. ${e^{\hat{i}}}_{\mu}$, on the frame indices. The signature is mostly plus.

The (first-order) action is expressed in terms of two basic independent fields: (i) the (inverse) vierbein $e^{i}_{\mu}$; and (ii) a general SO$(3,1)$ connection ${A^{i}}_{j\mu}$, which is metric-preserving (i.e. $A_{ij\mu}=-A_{ji\mu}$, where $A_{ij\mu}\equiv\eta_{is}{A^{s}}_{j\mu}$). The most general ghost-free and tachyon-free (around Minkowski spacetime) action containing only a massless spin-2 excitation and a (positive-parity) massive spin-2 one has four parameters555See Appendix B for a discussion, and the link with our previous notation. and can be written as:

[TABLE]

The torsion bigravity part, $S_{\rm TBG}$, of the action reads

[TABLE]

where $\ \sqrt{g}\equiv\sqrt{-\det g_{\mu\nu}}\equiv\det{e^{i}}_{\mu}$, and

[TABLE]

Here, we use the letter $R$ to denote the various curvatures defined by the Riemannian structure (curvature tensor ${R^{i}}_{jkl}\equiv{R^{i}}_{j\mu\nu}{e_{k}}^{\mu}{e_{l}}^{\nu}$, Ricci tensor $R_{ij}={R^{k}}_{ikj}$ and curvature scalar $R=\eta^{ij}R_{ij}$), and the letter $F$ to denote the corresponding Yang-Mills-type curvatures defined by the SO$(3,1)$ connection ${A^{i}}_{j\mu}$ ((curvature tensor ${F^{i}}_{jkl}\equiv{F^{i}}_{j\mu\nu}(A){e_{k}}^{\mu}{e_{l}}^{\nu}$, Ricci tensor $F_{ij}(A)={F^{k}}_{ikj}$ and curvature scalar $F(A)=\eta^{ij}F_{ij}$). Note that, because of the projections on the frame, the frame components of the $F$-type curvature depend algebraically on the vierbein ${e^{i}}_{\mu}$, besides depending on ${A^{i}}_{j\mu}$ and its first derivatives. See Appendix A for more details on the definition of these objects, and for the relation with the notation used in our previous paper Nikiforova:2018pdk . [An explicit form of the general field equations can also be found in the latter reference.]

The torsion bigravity lagrangian (II) a priori depends on four parameters: $c_{R}$, $c_{F}$, $c_{F^{2}}$, and $c_{34}$. Actually, the last one, $c_{34}$, will not enter in the discussion of spherically symmetric solutions. This leaves us with three relevant parameters. The analysis of Refs. Hayashi:1980qp ; Sezgin:1981xs has shown that the absence (around a Minkowski background) of pathologies (ghosts or tachyons) require the three parameters $c_{R}$, $c_{F}$, $c_{F^{2}}$ to be positive. Actually, they are related to the gravitationlike coupling constants $G_{0}$ (linked to massless spin-2 exchange) and $G_{m}$ (linked to massive spin-2 exchange), and to the mass666Here, the “mass”, $\kappa$, of the massive spin-2 field refers to the inverse of its (reduced) Compton wavelength, i.e. the parameter entering the exponential decay $\propto{\rm e}^{-\kappa r}$ of a static torsion field. $\kappa\equiv m_{2}$ of the massive spin-2 excitation, by the relations

[TABLE]

Here, we have introduced (following Sezgin:1979zf ) the notation ${\lambda}$ for the sum $c_{R}+c_{F}$ of the two curvature coefficients. It is indeed this sum which measures (at least in the weak field limit) the usual Einsteinian gravitational coupling constant $1/(16\pi G_{0})$. We have also introduced the notation $\eta$ for the dimensionless ratio $c_{F}/c_{R}$, which measures (within a factor $\frac{4}{3}$ linked to the difference between the massless, $S_{0}^{\mu\nu}=T^{\mu\nu}-\frac{1}{2}T\eta^{\mu\nu}$, and massive, $S_{m}^{\mu\nu}=T^{\mu\nu}-\frac{1}{3}T\eta^{\mu\nu}$, spin-2 matter couplings777In the Newtonian limit, we have, indeed, $S_{0}^{00}=\frac{1}{2}T^{00}$ while $S_{m}^{00}=\frac{2}{3}T^{00}=\frac{4}{3}S_{0}^{00}$.) the ratio of couplings to matter. It is tempting to conjecture that, for general solutions, the ultraminimal class of theories defined by the three parameters $G_{0}$, $G_{m}$ and $\kappa=m_{2}$, taking $c_{34}=0$, will have the best possible nonlinear behavior.

The difference between the affine connection ${A^{i}}_{j\mu}$ and the torsionless Levi-Civita connection ${\omega^{i}}_{j\mu}(e)$ defined by the vierbein $e^{i}_{\mu}$ is called the contorsion tensor

[TABLE]

The frame components ${K^{i}}_{jk}\equiv{e_{k}}^{\mu}{K^{i}}_{j\mu}$ of the contorsion tensor are related to the frame components ${T^{i}}_{[jk]}=-{T^{i}}_{[kj]}$ of the torsion tensor by the relations

[TABLE]

[Note that ${T}_{i[jk]}=-{T}_{i[kj]}$ while ${K}_{ijk}=-{K}_{jik}$.] The field equations are linear in the second-order derivatives of $e^{i}_{\mu}$ and ${A^{i}}_{j\mu}$ when using these quantities as basic fields in the action. One should avoid to use the vierbein and the torsion as basic fields because this introduces, in view of the link (5) which involves first derivatives of the vierbein, third derivatives of the vierbein in the field equations. One should rather consider the torsion as a field that is a posteriori derived from the basic fields.

Let us emphasize that the first-order formalism used in the Einstein-Cartan(-Weyl-Sciama-Kibble) theory considered here (which is often called “Poincaré gauge theory”) is fundamentally different from the often considered Palatini-type (“metric-affine”) first-order formalism. In both formalisms one independently varies the metric and the connection, and one a priori allows for the presence of torsion, i.e. for a nonsymmetric part of the connection:

[TABLE]

However, in the Palatini approach (which is usually performed in a coordinate frame) one independently varies all the components of a (symmetric metric) $g_{\mu\nu}$ and of a (non-symmetric) connection ${\Gamma^{{\lambda}}}_{\mu\nu}$. This yields $10$ equations obtained by varying $g_{\mu\nu}$ together with $4^{3}=64$ equations obtained by varying the connection ${\Gamma^{{\lambda}}}_{\mu\nu}$. By contrast, in the Cartan-type approach used here, one gets $16$ equations by varying $e^{i}_{\mu}$ and only $24$ equations by varying $A_{ij\mu}=-A_{ji\mu}$. Because of the (chosen) Local-Lorentz invariance of the action the $16$ vierbein equations are submitted to $6$ Noether identities (linked to infinitesimal local Lorentz rotations $\omega_{[ij]}$; see, e.g., Hehl:1976kj ; Hayashi:1979wj ; Nikiforova:2017saf ) and are therefore essentially equivalent to $10$ field equations obtained by varying $g_{\mu\nu}$. By contrast, the $64$ connection equations of the Palatini approach are stronger than the $24$ equations obtained by varying $A_{[ij]\mu}$. For instance, if the connection does not directly couple to matter, it has been shown Afonso:2017bxr ; BeltranJimenez:2017doy that a general Palatini action of the $\sqrt{g}f({\cal R}_{(\mu\nu)})$ type (where ${\cal R}_{(\mu\nu)}$ denotes the symmetric part of the Ricci tensor defined by the nonsymmetric connection ${\Gamma^{{\lambda}}}_{\mu\nu}$) yields algebraic equations for the connection that determine it (modulo an additional “projective” term $\delta^{\lambda}_{\mu}A_{\nu}$) to be the torsionless Levi-Civita connection of the auxiliary gothic metric $\sqrt{q}q^{\mu\nu}\equiv\delta\left[\sqrt{g}f({\cal R}_{(\mu\nu)})\right]/\delta{\cal R}_{(\mu\nu)}$. As the projective term drops out of the action (because it does not contribute to ${\cal R}_{(\mu\nu)}$ and is assumed not to couple directly to matter) one ends up with a theory of gravity where the metric $q_{\mu\nu}$ is an Einstein-frame metric having the usual Einstein-Hilbert dynamics, but where the matter is coupled to the different metric $g_{\mu\nu}$, with some nonlinear relation between these two metrics and the matter stress-energy tensor $T_{\mu\nu}$. In these theories, there are no dynamical effects linked to a propagating torsion. On the other hand, in the generalized Cartan-type theories considered here, the torsion field is a dynamical field, which is generated by the matter stress-energy tensor $T_{\mu\nu}$ even in absence of direct coupling of the connection to matter, which propagates away from the material sources, and which has physical effects via its coupling to the physical metric $g_{\mu\nu}$.

From the technical point of view, the crucial difference between the Cartan-type and Palatini-type approaches is that the SO$(3,1)$ connection $A_{[ij]\mu}$ is algebraically constrained to be metric-preserving. This means that, in order to derive the Cartan-type field equations within a coordinate-based Palatini approach one needs to add to the action density a Lagrange multiplier term , say

[TABLE]

where $Q_{{\lambda},(\mu\nu)}\equiv\nabla_{{\lambda}}^{\Gamma}g_{\mu\nu}$ denotes the covariant derivative of the metric with respect to the general ( a priori nonsymmetric) affine connection ${\Gamma^{{\lambda}}}_{\mu\nu}$. Note that the presence of this term in the action then contributes to the $64$ equations obtained by varying the connection by additional terms involving the $40$ unknown Lagrange multipliers $\Lambda^{{\lambda}(\mu\nu)}$.

III Static spherically symmetric metrics and connections

In the present paper, we investigate static spherically symmetric solutions of torsion bigravity. We assume from the beginning that the solutions are: (i) time-reversal invariant; (ii) SO(3) invariant; and (iii) parity invariant. Under these assumptions, we can use a Schwarzschildlike radial coordinate, and denote

[TABLE]

so that the metric reads:

[TABLE]

We then correspondingly define the co-frame $e^{\hat{i}}={e^{\hat{i}}}_{\mu}dx^{\mu}$ as

[TABLE]

The structure of a general (possibly torsionfull) connection under the just stated assumptions (i)–(iii) has been determined by Rauch and Nieh Rauch:1981tva . This structure is clear when using Cartesianlike coordinates $x^{0},x^{a}$, (with $a=1,2,3$) and a corresponding Cartesianlike co-frame $e^{\hat{0}},e^{\hat{a}}$. Time-reversal invariance implies that the only nonvanishing components of a general connection must form a vector $A_{\hat{a}\hat{0}\hat{0}}=-A_{\hat{0}\hat{a}\hat{0}}$ and a three-index tensor $A_{\hat{a}\hat{b}\hat{c}}=-A_{\hat{b}\hat{a}\hat{c}}$. Then spherical symmetry implies that the vector $A_{\hat{a}\hat{0}\hat{0}}$ must be in the radial direction $n^{a}$, say

[TABLE]

with some radial function ${\overline{V}}(r)$, while spherical symmetry, and parity invariance (which forbids the presence of the Levi-Civita tensor $\epsilon_{\hat{a}\hat{b}\hat{c}}$) imply that the three-index tensor $A_{[\hat{a}\hat{b}]\hat{c}}$ must be of the form

[TABLE]

with a second radial function ${\overline{W}}(r)$. Therefore the most general affine connection (under the assumptions (i)–(iii)) involves two a priori unknown radial functions. When re-expressing these results in terms of the polar-type frame (12), one finds that the two unknown radial functions parametrizing a general affine connection can be chosen as being

[TABLE]

Note that $V$ and $W$ are components along our basic orthonormal frame (12).

Then the nonvanishing components of the connection one-form are found to be

[TABLE]

Note that the last component (in the $\theta,\varphi$ 2-plane) is independent of the unknown functions $V,W$, but only depends on the use of a polar-type frame, with a Schwarzschildlike radial coordinate.

The nonzero components of the torsionless Levi-Civita connection one-form, ${\omega^{i}}_{j\mu}(e)$, defined by the metric (11), are found to be (using Eq. (243))

[TABLE]

Note that the last component is (as necessary) the same as for the general affine connection $A$, and that the nonzero components of the contorsion tensor are then found to be (modulo the antisymmetry with respect to the first two spatial indices in the second equation)

[TABLE]

Because of the restricted number of nonzero components, the nonzero components of the torsion tensor $T_{i[jk]}$ (which is antisymmetric with respect to the last two indices, are the same (modulo some permutation of indices) as those of the contorsion tensor $K_{ijk}=K_{[ij]k}$ (which is antisymmetric with respect to the first two indices), e.g.

[TABLE]

Using (III) we can construct the Einstein tensor of the $A$ connection:

[TABLE]

This tensor happens to be symmetric, $G_{ij}(A)=G_{ji}(A)$ under our (static, spherically symmetric) assumptions. Its nonzero components read,

[TABLE]

For additional clarity, we used here a more explicit notation for the frame indices:

[TABLE]

IV Torsion bigravity action

Using the previous formulas we can now write down the action, and derive from it the field equations. [We have checked that varying the spherically-symmetric-reduced action does yield field equations that are equivalent to the spherically-symmetric-reduced field equations, as derived directly from the general field equations in Ref. Rauch:1981tva .] We recall that the structure of the action is

[TABLE]

The variation of the matter action $S_{\rm{m}}$ with respect to the metric reads

[TABLE]

while we assume here that its variation with respect to the SO(3,1) $A$ connection (linked to the local, quantum, spin density) vanishes.

The field action is the sum of various contributions:

[TABLE]

Here (neglecting to write the “double-zero” term $\propto F_{[ij]}^{2}$)

[TABLE]

and

[TABLE]

where

[TABLE]

For notational simplicity, we shall often omit below to include in the action the trivial (field-independent) volume factor $dt\,(\sin\theta\,d\theta\,d\phi)$, so as to work with a radial action $S^{\prime}=\int drw(r)L$.

The usual Einstein-Hilbert term is explicitly computed as being

[TABLE]

which can be rewritten in the form

[TABLE]

where

[TABLE]

Note that in this form, the first term is linear in the first derivatives of the metric variables (actually linear in $\Lambda^{\prime}$). The affine-connection analog of the Einstein-Hilbert term is obtained by inserting Eqs. (III) in

[TABLE]

where we used the fact that

[TABLE]

To streamline the structure of the terms depending on the derivatives of $V$ and $W$, it is useful to introduce a shorthand notation for the kind of covariant derivatives of $V$ and $W$ entering Eqs. (III), namely

[TABLE]

We also introduce a shorthand notation for the term involving the square of $W$, namely

[TABLE]

With this notation, we have

[TABLE]

Concerning the contribution quadratic in $F_{ij}(A)$, it is easy to see that

[TABLE]

so that $L_{F^{2}}$ can be directly expressed in terms of $G_{ij}(A)$ as

[TABLE]

Inserting the expressions (III) for the components of $G_{ij}$, and using the shorthand notation introduced above, leads to

[TABLE]

At this stage, the various contributions to the action take the form

[TABLE]

A remarkable fact about this action is that the only term containing the square of derivatives is the contribution $\propto(\nabla V+\nabla W)^{2}$ in $L_{F^{2}}$. It is then convenient to add a so-called ”double-zero” term to the action, so as to end up with an equivalent action which is only linear in derivatives. [In the present case, this is also equivalent to making a Legendre transform.]

To explain the idea behind this transformation, let us first consider a toy model with the Lagrangian

[TABLE]

We can eliminate the square of the derivative of $q$ by adding the following double-zero term to the Lagrangian, involving a new, independent variable $\pi$:

[TABLE]

Indeed, the equation of motion of $\pi$ obtained by varying $L^{\rm{old}}_{\rm{toy}}+\Delta L$ is

[TABLE]

Then the modification of the equation of motion of $q$ coming from varying $\Delta L(\pi,\dot{q},q)$ will involve (because of the quadratic nature of $\Delta L(\pi,\dot{q},q)$) a factor $\left[\pi-(\dot{q}+A(q))\right]$, which vanishes when $\pi$ is on-shell. This shows that the action

[TABLE]

leads to equivalent equations of motion. But the latter action is first-order in derivatives. Indeed:

[TABLE]

On the last line we recognize the result of making a Legendre transformation from $\dot{q}$ to $2\pi=\delta L^{\rm{old}}_{\rm{toy}}/\delta\dot{q}$.

In our case, we choose to introduce as new variable the only combination of covariant derivatives of $V$ and $W$ that enters quadratically in the action, namely

[TABLE]

We then add to the original action the double-zero term

[TABLE]

which yields

[TABLE]

The field action then looks as follows,

[TABLE]

We use the macroscopic energy-momentum tensor,

[TABLE]

i.e., using $u^{0}=e^{-\Phi}$,

[TABLE]

so that the variation of the matter action reads

[TABLE]

V Torsion bigravity field equations

Let us now write down the equations obtained from varying the action $S^{\prime}_{\rm{field}}+S^{\prime}_{m}=\int drw(r)L$ (considered in its first-order form, with $\pi$ as an independent variable) with respect to the five field variables $x^{a}=(\Phi,\Lambda,V,W,\pi)$, $a=1,2,3,4,5$. Note that, introducing $x^{0}\equiv r$, as a fictitious sixth timelike variable, the latter first-order action has the structure

[TABLE]

Here, we denoted ${\dot{x}}^{a}=dx^{a}/dx^{0}$, and $x^{\mu}=(x^{0},x^{a})$, with $\mu=0,1,2,3,4,5$. The six components $A_{\mu}(x^{\nu})$ of the one-form $A_{\mu}(x)dx^{\mu}$ depend on the six variables $x^{\nu}$. The one-form $A_{\mu}(x)dx^{\mu}$ is just the usual Hamilton-Cartan one-form $p_{a}dq^{a}-Hdt$ of a first-order action, but we find useful to view it as the Maxwelllike action for a massless charged particle of worldline $x^{\mu}$ interacting with an external electromagneticlike potential $A_{\mu}(x^{\nu})$.

Let us write separately the contributions coming from varying the various pieces of the action $S^{\prime}=\int drw(r)L$ with respect to the five field variables $x^{a}=(\Phi,\Lambda,V,W,\pi)$.

[TABLE]

Here we have introduced (after variation) the shorthand notation $F$ for the radial derivative of $\Phi$:

[TABLE]

We use as basic equations for the five field variables $x^{a}=(\Phi,\Lambda,V,W,\pi)$ the five first-order equations

[TABLE]

with (denoting $c\equiv c_{F}$, so that ${\lambda}-c\equiv c_{R}$)

[TABLE]

where each term $\delta(w\,L)/\delta x^{a}$ is obtained by summing the corresponding terms among Eqs. (55)-(67). [The factors $\kappa^{2}(c-{\lambda})=-\kappa^{2}c_{R}$ have been included to eliminate the denominator implicitly present in $c_{F^{2}}=\frac{\eta{\lambda}}{\kappa^{2}}=\frac{c{\lambda}}{({\lambda}-c)\kappa^{2}}$.]

The five (geometric) field equations above must be supplemented (when considering the interior of a star) by the usual (universal) matter equation following from the (radial) conservation law $\nabla^{g}_{\mu}T^{\mu\nu}=0$ for a spherically-symmetric configuration with macroscopic energy-momentum tensor, Eqs. (IV), namely

[TABLE]

with

[TABLE]

VI Reduction of the field equations to a ghost-freelike system of three first order equations

Let us recall that the basic aim of the present work is to study the geometric torsion bigravity model as an alternative to the usually considered bimetric gravity models. The latter models are defined by considering two independent dynamical metric tensors, say $g_{\mu\nu}$ and $f_{\mu\nu}$, having separate Einstein-Hilbert actions, and being coupled to each other (besides some matter coupling) via some generalized Fierz-Pauli potential ${\cal V}(f,g)$. These models are generalizations of the massive gravity models where the metric $f_{\mu\nu}$ is non dynamical, and frozen into some given background value (e.g. a Minkowski background $f_{\mu\nu}=\eta_{\mu\nu}$). For many years, it was thought that massive gravity models (and, consequently, their bimetric generalizations) were plagued by the necessary presence of an additional, ghostlike, degree of freedom Boulware:1973my ; Deffayet:2005ys ; Creminelli:2005qk . The latter Boulware-Deser ghost enters only at the nonlinear level (because, at the linear level, the Fierz-Pauli potential Fierz:1939zz ; Pauli:1939xp ; Fierz:1939ix ensures the presence of only five, healthy degrees of freedom in the massive-gravity sector).

It was emphasized by Babichev, Deffayet and Ziour Babichev:2009us that the presence of the Boulware-Deser ghost in generic massive gravity models888We start by considering generic massive-gravity (and bimetric) models containing a Boulware-Deser ghost to contrast them with the properties of ghost-free massive-gravity (and bimetric) models. is already apparent when considering (co-diagonal) spherically symmetric solutions. More precisely, a generic massive gravity model has (when using a Schwarzschild radial coordinate $r$ for the physical metric $g_{\mu\nu}$) three variables: $\Phi(r)$, $\Lambda(r)$ (defined as in Eq. (11) above) together with a third “gauge” variable $\mu(r)$ relating the Schwarzschildlike radius $r$ to the “flat” radial variable $r_{f}$ defined by the background metric $f_{\mu\nu}$, namely $r_{f}=re^{-\mu(r)/2}$. The crucial point (which can also be seen in the explicit field equations of Ref. Damour:2002gp ) is that the massive-gravity field equations are first order in $\Phi(r)$, and $\Lambda(r)$, but second order in $\mu(r)$. This means that the total differential order of the massive-gravity $\Phi(r)$, $\Lambda(r)$, $\mu(r)$ system is four. Equivalently, the general999Here, “general” means that we do not impose boundary conditions at infinity. exterior spherically-symmetric solution of a generic massive-gravity model contains four arbitrary integration constants. One of them will be an additional constant $c_{0}$ in $\Phi(r)$, which is physically irrelevant because it can be gauged away by renormalizing the time variable: $t\to t^{\prime}=e^{-c_{0}/2}t$. We conclude that the general exterior spherically-symmetric solution of a generic (ghostfull) massive-gravity model contains three physically relevant arbitrary integration constants. This is one more constant than for the general exterior spherically-symmetric solution of the Fierz-Pauli linearized massive gravity model. Indeed, the latter general linearized solution for $h_{\mu\nu}=g_{\mu\nu}-f_{\mu\nu}$ is (see Boulware:1973my )

[TABLE]

where $\kappa$ denotes the mass of the massive graviton, and where $Y_{\kappa}(r)$ is the general exterior spherically-symmetric solution of the Yukawa equation

[TABLE]

which contains two integration constants, $c_{+},c_{-}$, namely

[TABLE]

We recall in passing that the trace of $h_{\mu\nu}$ is locally related to the matter density via

[TABLE]

Summarizing: the presence of a sixth field degree of freedom in a generic (ghostfull) massive-gravity model is visible when considering the general exterior spherically-symmetric solution: indeed, the latter solution, generically involves three physically relevant integration constants, which is one more than the two physically relevant integration constants $c_{+}$, $c_{-}$ entering the corresponding linearized solution of the five-degree-of-freedom Fierz-Pauli model. In addition, we recall that the linearized massive-gravity solution necessarily involves a $\frac{1}{\kappa^{2}}$ factor in some of its components, and that this feature is the origin of the appearance of a Vainshtein radius below which one cannot trust the usual weak-field perturbation expansion of massive gravity Vainshtein:1972sx ; Damour:2002gp .

Let us emphasize that the ability of the spherically-symmetric limit to detect the presence of the Boulware-Deser ghost is somewhat obscured if one focusses, from the beginning, on exponentially decaying solutions, rather than on general exterior solutions. [See, in this respect, Refs. Damour:2002gp , Babichev:2009us and Babichev:2010jd .]

When extending a massive-gravity model into a corresponding bimetric gravity one, we must add to the count of the physically relevant integration constants entering a general exterior solution the Schwarzschildlike mass $m$ parametrizing the physics of the massless- spin-2 sector. We therefore conclude that the general exterior solution of a ghostfull bimetric gravity model will involve four physically relevant integration constants, while the general exterior solution of a ghost-free bimetric gravity model will involve only three physically relevant integration constants (corresponding to $m$, $c_{+}$, $c_{-}$ parametrizing the corresponding linearized system). [We recall that we discounted here the physically irrelevant additional constant entering $\Phi(r)$.] The fact that the general exterior solution of ghost-free bimetric gravity models (using the restricted class of potential ${\cal V}(f,g)$ discovered in deRham:2010kj ) indeed involves only three physically relevant integration constants has been explicitly shown by Volkov Volkov:2012wp . Indeed, he showed how to reduce the (co-diagonal) field equations to a system of three first-order differential equations, for the three variables $N$, $Y$, and $U$, see Eqs. (5.7) in Volkov:2012wp . [The variable $\Phi=\ln Q$ is then obtained by a quadrature: $\Phi=\int dr{\cal F}_{5}+c_{0}$, see Eq. (5.3c) in Volkov:2012wp .]

We are now going to show that the torsion bigravity model is similar to the ghost-free bimetric gravity models in that its general exterior spherically-symmetric solution only involves three physically relevant integration constants. [We will see later that these three integration constants do correspond to the constants $m$, $C_{+}$, $C_{-}$ parametrizing the corresponding linearized torsion bigravity system.] This will be shown by reducing the system of five first-order field equations $E_{1}$–$E_{5}$ written in the previous section to a system of three first-order differential equations (together with a quadrature for $\Phi(r)$). In view of the fact, recalled above, that the presence of the Boulware-Deser ghost was visible in spherically-symmetric solutions of generic ghostfull bimetric gravity models, we consider this property of torsion bigravity as a suggestion (though not a proof) that it might be ghost-free in a general (time-dependent and non-spherically-symmetric) situation.

As our reduction process is algebraically involved, we will not display all the technical details, but only explain the algorithm by which we could explicitly derive a reduced system of three first-order equations for three unknowns. Explicit calculations are anyway better done by using algebraic manipulation programmes, starting from the explicit basic field equations written in the previous section.

Before explaining the explicit reduction process we used, let us briefly indicate how the reduction issue could be formulated in terms of the Hamilton-Cartan action (54). The variational equations of motion coming from the first-order action (54) are

[TABLE]

where $F_{\mu\nu}(x)=\partial A_{\nu}(x)/\partial x^{\mu}-\partial A_{\mu}(x)/\partial x^{\nu}$ are the components (with $\mu=0,1,2,3,4,5$) of the two-form $F=dA$, and where we recall that $x^{0}$ simply denotes the radial variable $r$, which plays the role of time in our action. Because of the antisymmetry of $F_{\mu\nu}$, there are only five independent equations among the equations ${\cal E}_{\mu}$, Eq. (81) (say ${\cal E}_{a}$, for $a=1,2,3,4,5$). A necessary condition for the variational equations (81) to have nontrivial solutions in the phase-space “velocity” $v^{\mu}=\frac{dx^{\nu}}{dx^{0}}$ is that the determinant of the six-by-six matrix $F_{\mu\nu}$ be vanishing. As $F_{\mu\nu}$ is antisymmetric and even, its determinant is the square of its Pfaffian

[TABLE]

This shows that a necessary condition following following from the five equations ${\cal E}_{a}=0$ (which are equivalent to the equations $E_{a}=0$ of the previous section) is the constraint

[TABLE]

The latter constraint is purely algebraic in the five variables $x^{a}=(\Phi,\Lambda,V,W,\pi)$ (and depends on $x^{0}=r$). In turn, the (primary) constraint (83) implies as secondary constraint an equation linear in the velocities $v^{a}=\frac{dx^{a}}{dx^{0}}$ (i.e. the radial derivatives of $(\Phi,\Lambda,V,W,\pi)$), namely

[TABLE]

This argument indicates that the basic system of five equations $E_{1}$–$E_{5}$ of the previous section implies (at least) one algebraic constraint, (83), together with the extra differential condition (84). To check what is the precise import of these constraints on the number of free data determining the general exterior solution of our system we need to explicitly write down and study these constraints, as we will do next (starting directly from the explicit form (70)–(74) of our five basic equations $E_{1}$–$E_{5}$).

When doing so, it it convenient to start by noticing that the gauge symmetry $t\to t^{\prime}=e^{-c_{0}/2}t$, which corresponds to changing $\Phi(r)$ into $\Phi(r)+c_{0}$ shows that our basic five field equations can be entirely expressed in terms of $F(r)\equiv\Phi^{\prime}$, without any explicit appearance of the undifferentiated variable $\Phi(r)$. Actually, the various $e^{-\Phi}$ factors in our definitions (70)–(74) were designed to realize this disappearance of $\Phi(r)$. In other words, we can consider the system $E_{1}$–$E_{5}$ as being algebraic in $F$, and differential (of first order) only in the four variables $\Lambda$, $V$, $W$, and $\pi$.

It is also useful to work with a slightly modified set of variables. In the following we shall replace the set of variables $F$, $\Lambda$, $V$, $W$, and $\pi$ by the new set $F$, $L$, $V$, $\overline{Y}$, and $\pi$ where

[TABLE]

where we used also the intermediate notation

[TABLE]

The usefulness of this change of variables is that it allows one to easily show that two combinations of our five basic equations $E_{1}$–$E_{5}$, (70)–(74), yield two algebraic constraints in the five variables $F$, $L$, $V$, $\overline{Y}$, $\pi$.

On the one hand, the equation $E_{1}$ turns out to be algebraic in $F$, $L$, $V$, $\overline{Y}$, $\pi$ (without involving any derivative):

[TABLE]

Moreover $E_{1}$ is linear in $F$ and quadratic in $L$. As indicated, $E_{1}$ also involves the pressure $P(r)$ as a matter source. The constraint $E_{1}=0$ will be used to algebraically eliminate $F$ by expressing it in terms of the other variables.

On the other hand, the only derivative entering the two equations $E_{3}$ and $E_{5}$ is $\overline{Y}^{\prime}$. This implies that a linear combination of $E_{3}$ and $E_{5}$ yields an algebraic constraint. More precisely the new expression

[TABLE]

is an algebraic expression in our (redefined) variables, namely

[TABLE]

which is is linear in both $F$ and $L$.

The reduction process we use is then the following. First, we solve the algebraic constraint $E_{1}(F,L,V,\overline{Y},\pi)=0$ (which is linear in $F$) with respect to $F$ to get

[TABLE]

Then, we replace $F\to F_{\rm sol}[L,V,\overline{Y},\pi;P]$ in the other algebraic constraint $E_{35}$, Eq. (90), to get a reduced algebraic constraint involving only the four geometric variables $L,V,\overline{Y},\pi$, say

[TABLE]

The so-obtained algebraic constraint $E^{\rm red}_{35}(L,V,\overline{Y},\pi)=0$ turns out to be quadratic in $L$. There is a unique root of this quadratic equation in $L$ 101010This is the smallest root, i.e. the root with a negative coefficient in front of the discriminant when writing the equation with a positive coefficient for $L^{2}$., say

[TABLE]

which is such that it has the physically desirable feature of asymptotically behaving like its Schwarzschild counterpart

[TABLE]

when the arguments $V,\overline{Y},\pi$ asymptotically decay at infinity in a Schwarzschildlike manner. This requirement follows from the physical requirement that the contorsion tensor (being entirely generated, at the linear level, via a massive-spin-2 excitation; see below) must decay $\propto e^{-\kappa r}$ so that $V$ and $W$, and the corresponding $\overline{Y},\pi$, must asymptotically decay as their Schwarzschild counterparts, i.e. as the corresponding frame components of the Levi-Civita connection, see (III).

Then, by substituting $L\to L^{-}_{\rm sol}[V,\overline{Y},\pi;P]$, from Eq. (93), into the expression (95), we get an explicit expression for $F$ in terms of the three geometric variables $V,\overline{Y},\pi$, say

[TABLE]

The final stage of our reduction process consists in replacing $F\to F^{\rm red}_{\rm sol}[V,\overline{Y},\pi]$ and $L\to L^{-}_{\rm sol}[V,\overline{Y},\pi]$ into the remaining equations $E_{2}$, $E_{4}$, and $E_{5}$ to get three first-order equations for the three unknowns $V$, $\overline{Y}$, and $\pi$ (involving also $P$ and $\rho$ as source terms), say

[TABLE]

By construction, these three equations are linear in all the radial derivatives. When replacing the radial derivative of the pressure which appears in $E_{2}^{\rm red}$ by the matter equation (98) discussed next, one can solve the three equations (VI) for the three derivatives $V^{\prime},\overline{Y}^{\prime},\pi^{\prime}$ so as to get an explicit first-order radial-evolution system, say

[TABLE]

When considering the solution inside a star one must augment this system by the reduction of the equation (76) constraining the radial evolution of the pressure, namely

[TABLE]

and by giving an equation of state relating $\rho$ to $P$, say $\rho=\rho(P)$.

After integrating the system (VI), (98), for the four variables $V,\overline{Y},\pi,P$, one can compute the values of the variables $F$, $L$ (or $\Lambda$), and $W$ by using Eqs. (95), (93), and (85). Finally, the value of the gravitational potential $\Phi(r)$ is obtained by a quadrature:

[TABLE]

where we fixed the arbitrary additional constant in $\Phi$ by the requirement that $\Phi(r)\to 0$ at radial infinity.

VII Linearized approximation

Let us study the linearized approximation to our five basic field equations $E_{1}$–$E_{5}$,(70)–(74). We are, in particular, interested in understanding how the linearized solutions behave in the small-mass limit $\kappa\rightarrow 0$. In the next section we will then consider the second-order (post-linear) solutions. We will see that, both at the linear level, and at the postlinear level, the limit $\kappa\rightarrow 0$ of torsion bigravity is much better behaved than in massive gravity and bimetric gravity. Some aspects of the linearized approximation of dynamical-torsion models have already been considered in Refs. Hayashi:1980ir ; Nikiforova:2009qr , and in Ref. Zhang:1982jn for spherically symmetric solution, but our treatment will be more extensive and detailed.

In absence of material source (i.e. when $\rho\to 0$ and $P\to 0$), the torsion bigravity field equations admit the solution $\Phi=0$, $F\equiv\Phi^{\prime}=0$, $\Lambda=0$, $V=0$, $W=-\frac{1}{r}$, $\overline{Y}\equiv V+W+\frac{1}{r}=0$. We denote with a subscript 1 a first-order deviation from this trivial solution, i.e. $F_{1}$, $\Lambda_{1}$, $V_{1}$ and $\overline{Y}_{1}$. The explicit form of the linearized approximation of the field equations looks as follows,

[TABLE]

The hat added on the ${\hat{E}}^{\rm{lin}}_{n}$’s indicate that these equations differ by a factor from the linearization of the corresponding equations $E_{n}$, as defined in Eqs. (70)–(74) above. In keeping with what was already the case at the nonlinear level, the first equation ${\hat{E}}^{\rm{lin}}_{1}$ is algebraic in the variables $\Lambda_{1}$, $F_{1}$, $V_{1}$, $\overline{Y_{1}}$ and thus can be used to express $F_{1}$ in terms of other three, $F_{1}=F_{1}\left(\Lambda_{1},V_{1},\overline{Y_{1}}\right)$. Furthermore, the equation ${\hat{E}}^{\rm{lin}}_{5}$ is algebraic in $\pi_{1}$ and can be used to express $\pi_{1}$ in terms of the variables $V_{1}$, $\overline{Y_{1}}$, and $\overline{Y_{1}}^{\prime}$. Henceforth, we solve ${\hat{E}}^{\rm{lin}}_{1}=0$ for $F_{1}$, and ${\hat{E}}^{\rm{lin}}_{5}=0$ for $\pi_{1}$ so as to eliminate

[TABLE]

from the system.

It is then easily seen that the replacement of Eq. (106) in Eq. (102) eliminates the derivative $\overline{Y_{1}}^{\prime}$ and yields an equation which is algebraic in $\Lambda_{1}$, $V_{1}$, $\overline{Y_{1}}$. We can then use the latter algebraic equation (which is equivalent to the combination ${\hat{E}}^{\rm{lin}}_{35}\equiv{\hat{E}}^{\rm{lin}}_{3}-{\lambda}\,{\hat{E}}^{\rm{lin}}_{5}$) to express $\Lambda_{1}$ in terms of $V_{1}$ and $\overline{Y_{1}}$, say

[TABLE]

After inserting all the replacements Eqs. (105), (106), (107), one ends up with two remaining equations to solve: Eq. (103), which is second order in $\overline{Y_{1}}$ and first-order in $V_{1}$, and Eq. (101), which is first order in $\overline{Y_{1}}$, and $V_{1}$. The explicit form of the latter two equations is streamlined by introducing the new variables,

[TABLE]

We find that these variables must satisfy the following equations

[TABLE]

Given a solution of these two linear equations, the full linearized solution is given by the inverse of Eqs. (108), (109), i.e.

[TABLE]

as well as by

[TABLE]

The total differential order of the system Eqs. (110), (LABEL:eqVmk) is three, i.e. the same order as we found above for the full, nonlinear system.

One should note the remarkable fact that these linearized-approximation equations never involve the inverse of the squared mass $\kappa^{2}$ of the massive spin-2 excitation. This is in sharp contrast with the corresponding linearized massive-gravity, or bimetric gravity, equations which always involve an inverse power of $\kappa^{2}$, see, e.g., Eqs. (VI), (80). We will see below that the absence of inverse powers of $\kappa^{2}$ persists at the postlinear order.

Let us recall the structure of the solutions of equations of type (110) and (LABEL:eqVmk), with general source terms on the right-hand sides,

[TABLE]

These equations have unique solutions that are regular at the origin and decaying at infinity. They are given by the following formulas,

[TABLE]

In the second equation, the Green’s function $G_{\kappa}(r,\hat{r})$, satisfying the equation

[TABLE]

is constructed as

[TABLE]

where $\theta(x)$ denotes Heaviside’s step function, while

[TABLE]

are two appropriate homogeneous solutions, incorporating the boundary conditions. Namely, $X_{>}(r)$ decays at infinity, while $X_{<}(r)$ is regular at $r=0$. In addition,

[TABLE]

is the appropriate (constant) Wronskian of the two solutions.

Note that $V_{m0}$ and $V_{mk}$ are “pure” variables corresponding to the massless and massive linear excitations respectively. We then see on Eqs. (VII), (113), (114), how each metric or connection variable is some combination of these two pure variables.

Let us explicitly display the above linearized solution in the simple case where the source is a constant density star, say

[TABLE]

But, first, let us note that the source terms in the linearized equations (110), (LABEL:eqVmk) have different perturbative orders of magnitude. Indeed, we can consider that the primary source of all the variables is the matter density $\rho(r)$, and that it defines the formal expansion parameter $\varepsilon$ of our weak-field expansion: $\rho\equiv\varepsilon\rho_{1}$. Here, $\varepsilon$ is a bookkeeping device, which will be set to one at the end. The linearized variables $\Phi_{1}$, $F_{1}$, $V_{1}$, etc. are first-order in $\varepsilon$. E.g. $\Phi=\varepsilon\Phi_{1}+O(\varepsilon^{2})$, $F=\varepsilon F_{1}+O(\varepsilon^{2})$ (where $F_{1}=\Phi_{1}^{\prime}$), etc. On the other hand, the pressure-gradient equation (76) has the structure

[TABLE]

The boundary condition that $P(r)$ vanishes at the surface of the star then shows that the pressure $P$ is actually of second order in $\varepsilon$: $P=P_{2}\varepsilon^{2}+O(\varepsilon^{3})$, with

[TABLE]

To determine $P_{2}$ we must first determine the value of $F_{1}$ generated by $\rho\equiv\varepsilon\rho_{1}$. We shall take into account, in the next section, the second-order effects induced by the source terms involving the pressure $P=P_{2}\varepsilon^{2}+O(\varepsilon^{3})$ in the linearized equations (110), (LABEL:eqVmk), (113), (114). In the present section, we can define the pure linearized fields $F_{1}$, $V_{1}$, etc. by neglecting all the pressure-related source terms in the field equations (110), (LABEL:eqVmk), (113), (114), and by using the constant density ansatz (123). This leads to the following explicit solutions of the system (110), (LABEL:eqVmk):

[TABLE]

[TABLE]

Here, we recall that ${\lambda}\equiv 1/(16\pi G_{0})$, while $R_{s}$ denotes the radius of the star, and we have defined

[TABLE]

The “form factor” $\mathcal{F}(z)$, entering the magnitude $C_{mk}$ of $V_{mk}$ outside the star, has been defined so that $\mathcal{F}(z)\to 1$ when its argument $z=\kappa R_{s}\to 0$.

There are apparent factors $\propto 1/\kappa^{3}$ entering the inner solution for $V_{mk}$. However, these factors (which come from the Wronskian $W=\kappa^{3}$ in the Green’s function) are cancelled by $O(\kappa^{3})$ terms in the numerators. Indeed, the Green’s function itself is seen to have a finite limit as $\kappa\to 0$, because

[TABLE]

This ensures that the linearized solution has a finite limit when $\kappa\rightarrow 0$ (at a fixed value of $r$). In the limit $\kappa\to 0$ (keeping fixed both $R_{s}$ and $r$) one has indeed the following limit for $V_{mk}$,

[TABLE]

Let us also give the expression for $F_{1}$,

[TABLE]

where

[TABLE]

and

[TABLE]

The full (interior and exterior) solutions for the other variables are easily derived from the expressions given above. Let us only write down here the exterior ($r\geq R_{s}$) solutions for all the variables. [We recall in passing that all variables have zero background values, except for $W=-\frac{1}{r}+\varepsilon W_{1}+O(\varepsilon^{2})$.]

[TABLE]

where

[TABLE]

It is important to display also the linearized values of the two independent components of the contorsion (and torsion), as defined in Eq. (III)

[TABLE]

They read

[TABLE]

[TABLE]

Note that the (con-)torsion components are exponentially decaying. [This remains true at all orders of perturbation theory.] By contrast, the geometric variables $\Phi_{1},F_{1},\Lambda_{1},V_{1},W_{1}$ contain an additive mixture of massless (power-law decaying) and massive (exponentially decaying) spin-2 excitations.

VIII Second order perturbations

Let us consider the solutions of torsion bigravity at the second order in the source $\rho=\varepsilon\rho_{1}$ (for the case of a constant density star: $\rho(r)=e_{0}$). Each variable (except $\rho=e_{0}$ itself which is left unexpanded) is now written as

[TABLE]

At second order, we define the second-order values of the functions $V_{m0}$ and $V_{mk}$, by (conventionally) using the same formulas as at first order, i.e.

[TABLE]

We can use the inverse of these equations (see Eqs. (VII)) to express $V_{2}$ and $\overline{Y}_{2}$ in terms of $V_{m0(2)}$ and $V_{mk(2)}$.

When expanding to second order our basic field equations $E_{1}$–$E_{5}$,(70)–(74), we first get algebraic equations for $F_{2}$ and $\Lambda_{2}$ of the form

[TABLE]

where $N_{2}^{F}$ and $N_{2}^{\Lambda}$ are additional second-order contributions which are either quadratic in the first-order variables $F_{1},\Lambda_{1},V_{1},\overline{Y}_{1}$ (and, eventually, their derivatives), or linear in the pressure $P_{2}$.

We also get differential equations for $V_{m0(2)}$ and $V_{mk(2)}$ of the form

[TABLE]

where the second-order source terms $S_{m0(2)}$ and $S_{mk(2)}$ consist of terms bilinear in $V_{1}$, $\overline{Y}_{1}$, $F_{1}$, $\Lambda_{1}$, together with additional contributions linear in the pressure $P_{2}$ (remembering that $P$ is second-order, see section VII). We recall that $P_{2}$ is obtained by solving the matter equation

[TABLE]

with the condition that $P_{2}$ vanishes at the radius of the star $r=R_{s}$.

The second-order solution is then explicitly obtained by using our general Green’s function representation

[TABLE]

We found that it was possible to explicitly compute all the integrals generated by inserting the first-order solution in the source terms $S_{m0(2)}$, $S_{mk(2)}$ entering the latter second-order expressions. The final expressions involve, besides elementary functions, some exponential-integral functions ${\rm Ei}(-x)$ with various arguments proportional to $z=\kappa r$ or $z_{s}=\kappa R_{s}$. We recall that, with $x>0$,

[TABLE]

It would take too much space to display here in full detail the second-order solutions (both in the interior and in the exterior of the star) for all our variables. We will only display here the function of most physical importance at the second order, namely the variable $F_{2}$, which is the radial derivative of the second-order gravitational potential $\Phi=\varepsilon\Phi_{1}+\varepsilon^{2}\Phi_{2}+O(\varepsilon^{3})$. As we shall explicitly discuss below, this is indeed the only variable whose second-order value is needed to discuss the usual first post-Newtonian approximation. In addition, it is enough to know its value outside the star to discuss its phenomenological implications as a modification of the usual Schwarzschild metric outside a spherical mass distribution.

The full, second-order exterior solution $F_{2}$ has a rather complicated structure, which can, however, be explicitly displayed as follows:

[TABLE]

where $z\equiv\kappa r$, $z_{s}\equiv\kappa R_{s}$, and where the dependence on the source characteristics of the various coefficients can be expressed in terms of two form factors: the previously defined form factor $\mathcal{F}(z_{s})$, (131), and a new one denoted $\mathcal{E}(z_{s})$ and defined as

[TABLE]

With this notation, the various terms in Eq. (154) are:

[TABLE]

In order to better understand the structure of $F_{2}$, let us study it under the two limits: (i) $r\to\infty$ at fixed $\kappa$ (so that $z=\kappa r\to\infty$); and (ii) $\kappa\to 0$ at fixed $r>R_{s}$ (so that $z=\kappa r\to 0$ and $z_{s}=\kappa R_{s}\to 0$). The first limit studies the asymptotic structure of the solution at spatial infinity, while the second one would be the relevant one if (as is often done in massive-gravity studies) one would consider a Compton wavelength $\kappa^{-1}$ for the massive gravity excitation of cosmological size.

VIII.1 Limit $r\to\infty$ at fixed $\kappa$

Let us start by recalling that the first-order approximation to the exterior solution for $F=F_{1}+F_{2}+\cdots$ reads, according to (135), as follows,

[TABLE]

$F_{1}$ is the sum of a usual Newtonlike (and Schwarzschildlike) power-law contribution $m_{1}/r^{2}$, and of a decaying Yukawa contribution $\propto\partial_{r}\left(e^{-\kappa r}/r\right)=-e^{-\kappa r}(1+\kappa r)/r^{2}$. Let us now consider the spatial asymptotics $r\rightarrow\infty$ of the second-order exterior solution $F_{2}$. To this end, we must take into account the asymptotic behavior of the exponential integral ${\rm Ei}(-z)$ (when $z\to+\infty$)

[TABLE]

Using the latter asymptotic behavior, one concludes that $F_{2}$, (154), contains four types of terms with different behavior at infinity:

[TABLE]

As a consequence the leading terms in the limit $r\rightarrow\infty$ of $F_{1}+F_{2}$ read,

[TABLE]

where $m_{2}\equiv m_{2}(z_{s})$ is given by Eq. (156), while $C_{1}^{F}(z_{s})$, $C_{2}^{F}(z_{s})$ and $C_{LN}$ are given by Eqs. (137), (162),(157).

We see that if we define the total mass parameter $m$ of the star in torsion bigravity (in the Schwarzschild sense of $m=GM$, i.e. a length scale associated with the mass) as the coefficient of $1/r^{2}$ in $F(r)$, as $r\to\infty$ (i.e. $\Phi(r)\approx-m/r$ in this limit), we have

[TABLE]

Here, we set the bookkeeping parameter $\varepsilon$ back to 1 in the first two terms, but kept it in the error term as a reminder that there are higher-order contributions that are at least cubic in the matter-density source $\rho$.

Before looking at the value of $m_{2}$ in various limits, let us note that the term $\frac{2m_{1}^{2}}{r^{3}}$ is the second-order term in the $m/r$ expansion of a Schwarzschild solution, say $F_{S}(r)$, of mass $m_{1}$, indeed

[TABLE]

More generally, one can show by considering the structure of perturbation theory in torsion bigravity that, to all orders of perturbation theory, the asymptotic spatial behavior of the solution will be such that the two independent (con-)torsion components (III) are exponentially decaying (modulo power-law and logarithmic factors),

[TABLE]

As a consequence, the variables $\Phi,F,\Lambda,V,W$ will asymptotically approach (modulo exponentially small corrections) some Schwarzschildlike geometric data (for some mass parameter $m$)

[TABLE]

Let us look more closely at the value of the asymptotic mass $m=m_{1}+m_{2}+O(\varepsilon^{3})$, and in particular at its second-order contribution $m_{2}(z_{s})$. We recall that

[TABLE]

where $G_{0}=1/(16\pi\lambda)$ is the (conventionally defined) massless spin-2 gravitational constant, and where

[TABLE]

is the (conventionally defined) bare mass-energy of the constant-density star. We recall in this respect that in GR, the total Schwarzschild mass of a constant-density star is actually, simply given by the Newtonlike expression

[TABLE]

where $G_{N}$ denotes Newton’s gravitational constant. If we identify the torsion bigravity massless spin-2 gravitational constant $G_{0}=1/(16\pi\lambda)$ with Newton’s constant, $G_{N}$, we see that our first-order mass parameter $m_{1}$ (with units of length) is equal to the (full) general relativistic mass parameter $m_{\rm GR}$.

On the other hand, the second-order contribution to the torsion-bigravity mass reads

[TABLE]

where the form factor $\mathcal{E}(z_{s})$ (with $z_{s}=\kappa R_{s}$) was defined in Eq. (155).

We recall that the dimensionless parameter

[TABLE]

is a measure of the ratio between the coupling constant $G_{m}$ of the massive graviton and the coupling constant $G_{0}$ of the massless one. Therefore the ratio between $m_{2}$ and $m_{1}$ can be written as

[TABLE]

This expression is compatible with the idea that in the limit where $G_{m}/G_{0}\to 0$ (at fixed $\kappa$) the torsion degrees of freedom decouple from the matter so that torsion bigravity reduces to GR with $G_{N}=G_{0}$, and the total mass parameter $m=m_{1}+m_{2}+\cdots$ reduces to its general relativistic value (177).

It is interesting to discuss the physical consequences of the form factor $\mathcal{E}(z_{s})=\mathcal{E}(\kappa R_{s})$ entering $m_{2}$. It is easily checked that the form factor $\mathcal{E}(z_{s})=\mathcal{E}(\kappa R_{s})$ has the following properties: (i) in spite of the prefactor $z_{s}^{-5}$ in its definition, $\mathcal{E}(z_{s})$ is regular when $z_{s}\to 0$, and has the finite limit

[TABLE]

(ii) $\mathcal{E}(z_{s})$ is negative in the interval $0\leq z_{s}<z_{*}$, and positive for $z_{s}>z_{*}$, where $z_{*}\approx 1.6969326$; and (iii) $\mathcal{E}(z_{s})$ tends to zero like $+1/z_{s}^{2}$ when $z_{s}\to\infty$.

As a consequence of this behavior of the form factor $\mathcal{E}(z_{s})$ we have the following limiting value for $m_{2}$ as $\kappa\to 0$ (i.e. $z_{s}\to 0$)

[TABLE]

[We will discuss the small $\kappa$ limit in more details in the next section.] The negative value of $m_{2}$ in this limit is probably due to the fact that the massive-gravitational binding energy $-\frac{3}{5}G_{m}M_{\rm bare}^{2}/R_{s}$ (due to exchange of massive spin-2 excitations, in the small mass limit) dominates over other forms of binding energy (e.g. pressure-related energy).

Another limit is the limit of very heavy massive spin-2 excitation ($\kappa\to\infty$), i.e. of a very short-range modification of gravity, $\kappa^{-1}\ll R_{s}$. In this case the second-order correction to the mass parameter mass is found (as expected) to go to zero,

[TABLE]

VIII.2 Limit $\kappa\rightarrow 0$ with fixed $r>R_{s}$

Let us now study in more detail the limit where $\kappa$ becomes very small, i.e. where the Compton wavelength $1/\kappa$ is much larger than all the other scales of the problem (and notably $R_{s}$), being, e.g., of cosmological magnitude. This is the situation which is usually considered for massive gravity and bimetric gravity. As is well known since the work of Vainshtein Vainshtein:1972sx , the perturbation expansion of massive gravity (and bimetric gravity) involves negative powers of $\kappa^{2}$, which render the perturbative expansion invalid for radii $r$ smaller than some Vainshtein radius $R_{V}$ given, in generic (ghostfull) massive-gravity theories, by the formula

[TABLE]

More precisely, at the second-order approximation in $G$, the perturbative solution of the field equations of generic massive-gravity (and bimetric gravity) theories contain terms that fractionally modify the linear approximation, say $\Phi_{1}\sim m/r$ by terms of the type (see, e.g., Deffayet:2001uk )

[TABLE]

The latter expansion is performed in the domain $R_{s}<r\ll\kappa^{-1}$, in the limit where $\kappa^{-1}$ is much larger $R_{s}$. [In this domain, and in this limit, one does not see the Yukawa exponential decay $\propto e^{-\kappa r}$.]

By contrast, we found the rather remarkable fact that, when considering the same limit, no terms involving inverse powers of $\kappa$ enter the perturbative expansion of torsion bigravity (in the domain $R_{s}<r\ll\kappa^{-1}$) up to the second order included.

For instance, the second-order contribution to $F$, considered in this limit, takes the following form

[TABLE]

where the $\kappa$-dependent piece tends to zero a $\kappa\to 0$. We have shown that, similarly, all the other field functions in second order perturbation theory, i.e. $V_{2}$, $\overline{Y}_{2}$, $\Lambda_{2}$ have finite limits (i.e. contain no denominators $\propto 1/\kappa^{2}$) as $\kappa\to 0$. Such a result was a priori not all guaranteed because the field equations of torsion bigravity do contain denominators $\propto 1/\kappa^{2}$. Indeed, such denominators come from the fact that the coefficient $c_{F^{2}}$ of the $F_{ij}^{2}$ terms in the action is proportional to $1/\kappa^{2}$, see Eq. (II).

The absence of $O(1/\kappa^{2})$ terms at second order is due to a special cancellation. Let us explain it. We recall that the second-order variables $F_{2}$ and $\Lambda_{2}$ are expressed in terms of the second-order potentials $V_{m0(2)}$ and $V_{mk(2)}$ via the equations

[TABLE]

Here, the additional (nonlinear) terms $N_{2}^{F},N_{2}^{\Lambda}$ (which are bilinear in $V_{1}$, $\overline{Y}_{1}$, $F_{1}$, $\Lambda_{1}$ and their derivatives) do contain some $1/\kappa^{2}$ factors, but all these factors have a special structure: each monomial containing a factor $\kappa^{-2}$ simultaneously contains at least one power of $\overline{Y}_{1}$ or of one of its derivatives. Similarly, the potentials $V_{m0(2)}$ and $V_{mk(2)}$ satisfy the differential equations (VIII) where the source functions $S_{m0\,(2)}$ and $S_{mk\,(2)}$ consist of terms bilinear in $V_{1}$, $\overline{Y}_{1}$, $F_{1}$, $\Lambda_{1}$ and their derivatives. Again the latter bilinear expressions $S_{m0\,(2)}$, $S_{mk\,(2)}$ do contain some $1/\kappa^{2}$ factors, but the latter a priori dangerous (when $\kappa\to 0$) terms have the same special structure as $N_{2}^{F},N_{2}^{\Lambda}$. Each factor $\kappa^{-2}$ multiplies a monomial which is at least linear in $\overline{Y}_{1}$ or one of its derivatives.

In turn, the reason why the terms $\propto\kappa^{-2}\overline{Y}_{1}$, or $\propto\kappa^{-2}\overline{Y}_{1}^{\prime}$, $\cdots$, turn out to be innocuous in the limit $\kappa\to 0$ is that the variable $\overline{Y}_{1}$ happens to be of order $O\left(\kappa^{2}\right)$ as $\kappa\to 0$, so that $\kappa^{-2}\overline{Y}_{1}$ has a finite limit as $\kappa\to 0$. Indeed, from the definition (108) one gets that

[TABLE]

Then, using Eq. (116) and the derivative of Eq. (115), one can see that $\overline{Y}_{1}$ satisfies the following differential equation,

[TABLE]

At the linear level, the source terms $S_{m0(1)}$, $S_{mk(1)}$, read, according to Eqs. (110), (LABEL:eqVmk),

[TABLE]

so that the combination of source terms entering the equation for $\overline{Y}_{1}$ cancells:

[TABLE]

Finally, $\overline{Y}_{1}$ satisfies an equation whose right-hand-side is explicitly $O(\kappa^{2})$, namely

[TABLE]

This explains why $\overline{Y}_{1}$ is of order $O(\kappa^{2})$, thereby ensuring the absence of denominators $1/\kappa^{2}$ in the second-order solution.

It is not a priori clear whether this (or a similar) cancellation mechanism will continue to work at the third order of perturbation theory. [The specific property (190) does not seem to persist for $S_{m0(2)}$ and $S_{mk(2)}$.] We note that one cannot apply the same reasonings to the next (third) order of perturbations because the property (190) is not true for $S_{m0\,(2)}$ and $S_{mk\,(2)}$. This means that it is a priori possible that the perturbation theory will involve $1/\kappa^{2}$ factors in the third order. We leave the investigation of this subject to future work, and comment below on what would be the consequences of the presence of $1/\kappa^{2}$ factors at the third order of perturbation theory. For the time being, we shall continue studying the consequences of our results at the second order of perturbation theory.

IX Numerically constructing exact star solutions

In GR, it is possible to write down analytically the exact solution for the metric generated by a constant-density perfect fluid Schwarzschild:1916ae . Though the exterior Schwarzschild solution Schwarzschild:1916uq is an exact exterior solution of torsion bigravity (with zero contorsion), this is not true for the interior Schwarzschild solution. Indeed, as we saw in our perturbation theory analysis, the presence of a nonzero $T_{\mu\nu}$ in space, necessarily generates some nonzero contorsion field, i.e. a difference between the affine connection $A^{i}_{j\mu}$ and the Levi-Civita connection ${\omega^{i}}_{j\mu}$. And indeed, one can check that the interior Schwarzschild solution (with zero contorsion) does not satisfy the field equations of torsion bigravity.

As the analytic construction of an exact analytical solution of the complicated system of torsion bigravity spherically-symmetric field equations discussed in Sec. VI seems difficult, we have appealed to numerical methods to confirm the global existence of regular solutions of torsion bigravity satisfying the boundary conditions imposed in our perturbation theory. Let us recall that these boundary conditions are: (1) geometric regularity of all our variables at the origin $r\to 0$, and (2) decay of all our variables at spatial infinity $r\to\infty$.

We recall that the system of equations to be satisfied (in presence of matter) consists either of: (i) the original six field equations comprising $E_{1}$–$E_{5}$, together with the matter equation $E_{m}$ (knowing that this system is constrained by two other equations that must be satisfied); or (ii) a reduced system made of the three radial-evolution equations (VI), plus the radial-evolution equation (98) for the pressure $P(r)$. In our numerical simulations, we have used the reduced first-order system of four ordinary differential equations (ODE’s) defined by Eqs. (VI), (98), for the four variables $V,\overline{Y},\pi,P$. This system is completed by giving an equation of state for the matter. In our simulations we use the simple condition of constant density: $\rho(r)=e_{0}$. After integrating this system, the values of the variables $F$, $L$ (or $\Lambda$), and $W$ was obtained by using Eqs. (95), (93), and (85).

As we have seen in Sec. VII, in perturbation theory the boundary conditions (1) and (2) (together with the choice of the radius $R_{s}$ of the star) uniquely determine (at each order of perturbation theory) a torsion bigravity solution. The main motivation for constructing numerical solutions was to prove that this uniqueness property holds in the full nonlinear theory. To do this we need to study what the conditions of regularity at the origin impose as constraints on the initial conditions (at $r\to 0$) of our four variables $V(r),\overline{Y}(r),\pi(r),P(r)$. First, the geometric character (scalar, vector, tensor, etc) of our variables show that, near the origin, they must admit general Taylor expansions of the following restricted type:

[TABLE]

together with

[TABLE]

[$\Lambda(0)=0$ is necessary to have a locally flat metric at the origin.] By inserting these expansions into the equations of our system, we get, at each order in $r$ some relations between the various expansion coefficients. The crucial point is that, if we consider the central value $P_{0}=P(r=0)$ of the pressure as a given quantity (that will determine the radius, given the constant density $e_{0}$), the equations of our system give enough relations to determine all the other expansion coefficients $v_{n},y_{n},\pi_{n},P_{n}$ in terms of only one of them. We have chosen $v_{1}$ as unique free inital datum. For instance, at the lowest order in the $r$ expansion, one finds that $y_{1}$, $\pi_{0}$, $f_{1}$ and $\Lambda_{2}$ are determined by $v_{1}$ and $P_{0}$ (and $e_{0}$) by the following formulas

[TABLE]

Similar formulas also determine the next order coefficients in the $r$ expansion: $v_{3}$, $y_{3}$, $\pi_{2}$, $P_{2}$, etc.

In other words, a single “shooting parameter” at the origin, namely $v_{1}$, uniquely determines (after having chosen $P_{0}$) the solution of torsion bigravity. When integrating the system, the value $R_{s}(P_{0},v_{1})$ of the star radius will be obtained as the (first) radius where $P(r)$ (vanishes).

For $r>R_{s}(P_{0},v_{1})$ one sets $\rho(r)=0$ and $P(r)=0$ and continues integrating the three field equations (VI) to get the exterior solution for the three variables $V,\overline{Y},\pi$. For a generic value of $v_{1}$, the so-constructed exterior solution for $V,\overline{Y},\pi$ (and the associated values of $F$, $W$ and $\Lambda$) will not decay at infinity, but will contain some growing exponential pieces $\propto e^{+\kappa r}$. We have seen in Sec. VII that the general exterior solution contains three parameters: one parameter, say $m$, (Schwarschild-type total mass) parametrizing all the power-law behaviour of the solution (asymptotically described by a Schwarzschild metric and connection); together with two parameters, say $C_{+}$ and $C_{-}$ respectively parametrizing the exponentially growing, and decaying, Yukawa-type contributions to the solution. At the linear level, each variable contains different coefficients $C_{+}$ and $C_{-}$, e.g.,

[TABLE]

but all the exponential-mode coefficients are related between themselves by the field equations, so that only two of them are independent.

In order to satisfy the decaying boundary condition at spatial infinity, we finally have a one-parameter shooting problem, namely it is enough to impose that (given some value of $P_{0}$) the coefficient $C_{+}(v_{1})$ of one variable vanishes. To numerically extract from numerical data an estimate of the (common, underlying) $C_{+}(v_{1})$ coefficient, we worked with the variable $V_{mk}(r)\equiv 2V(r)-\overline{Y}(r)$ which does not contain a mass-type, power-law contribution. In practical terms, this meant tuning the value of $v_{1}$ at $r=0$ until reducing essentially to zero the value of, say,

[TABLE]

taken at some large value of $r_{0}$ (such that $e^{+\kappa r_{0}}\gg 1$, so that the exponentially decaying contribution to $V_{mk}(r_{0})$ is fractionally negligible). [In practice, we used $\kappa r_{0}=10$ corresponding to $e^{+\kappa r_{0}}\approx 2\times 10^{4}$]. The tuning of $v_{1}$ is obtained by a simple dichotomy procedure, i.e. alternating the signs of $C_{+}^{\rm eff}(r_{0};v_{1})$ by changing the value of $v_{1}$ until $C_{+}^{\rm eff}(r_{0};v_{1})$ is smaller than what is permitted by the numerical accuracy of our simulation.

We implemented this simple, one-parameter shooting strategy for several star models, of various radii and compactnesses. Let us only indicate here our results for one such star model. Without loss of generality, we used units where $\kappa=1$ and ${\lambda}=1$. The first condition says that we measure lengths in units of $\kappa^{-1}$, while the second one defines the (independent) unit for the Newtonian constant such that $16\pi G_{0}=1$. Here, we shall exhibit a specific star model having the following physical characteristics. First we set the dimensionless torsion bigravity parameter $\eta$ to the value $\eta=1$, i.e. $c_{F}=c_{R}$ (both being equal to $\frac{1}{2}$ in our units where ${\lambda}=c_{F}+c_{R}=1$). The other physical choices concern: (a) the radius of the star in units of $\kappa^{-1}$, i.e. the dimensionless quantity $z_{s}=\kappa R_{s}$, and (b) the value of the star compactness111111We normalize the definition of the compactness so that it is equal to 1 for a black hole in GR. See below the exact definition of $C_{s}$., $C_{s}\simeq 2G_{0}M_{s}/R_{s}$, with $M_{s}\equiv\frac{4\pi}{3}e_{0}R_{s}^{3}$. The two quantities $z_{s}$ and $C_{s}$ are dimensionless, and physically depend on the two independent values of $e_{0}$ and $P_{0}$. We have chosen (in our units) the specific values

[TABLE]

These values were chosen by using, as guideline, our perturbation-theory expressions, with the aim of getting a star model having $\kappa R_{s}\sim 1$, and a sufficiently high compactness $C_{s}\sim 0.3$ (comparable to the expected compactness of a neutron star in GR).

Anyway, after doing the choices (IX), we found that we needed to tune $v_{1}$ to the value

[TABLE]

to get a sufficiently small value of $C_{+}^{\rm eff}(r_{0})$, i.e. a solution exhibiting numerical decay up to $r\sim 10/\kappa$. As said above, we obtained $v_{1}^{\rm tuned}$ by dichotomy, using as first guesses the analytical estimates of $v_{1}$ obtained either directly from linearized perturbation theory, namely

[TABLE]

or, alternatively, by combining the relation between $v_{1}$ and the value $f_{1}$ of $F^{\prime}|_{r\to 0}$ with the analytical estimate for $f_{1}$ deduced from our linear linear solution (135), i.e.

[TABLE]

The numerical solution was found to have a star radius equal (in our units where $\kappa=1$)

[TABLE]

The value of the star radius was numerically determined by looking at the point where the pressure $P(r)$ vanishes.

We display in Fig. 1 the numerical values (both inside and outside the star) of four variables encapsulating the essential geometrical properties of our solution, namely: $F(r)$, $\Lambda(r)$, and the two independent (con)torsion components, namely ${K^{\hat{1}}}_{\hat{0}\hat{0}}$ and ${K^{\hat{1}}}_{\hat{2}\hat{2}}$, as defined in Eq. (III). While $F(r)$ and $\Lambda(r)$ decay for large $r$ in a power-law fashion ($F(r)\propto 1/r^{2}$ and $\Lambda(r)\propto 1/r$), the torsion components decay exponentially. Note that the order of magnitude of the torsion inside the star is comparable to the value of $F$. As ${K^{\hat{1}}}_{\hat{0}\hat{0}}=V-e^{-\Lambda}F$ (from (III)), we see that the matter density of the star generates a torsion which is of roughly the same magnitude as the component ${\omega^{\hat{1}}}_{\hat{0}\hat{0}}=e^{-\Lambda}F$ of the Levi-Civita connection. [From Eqs. (VII), this remains true even when $\eta\to 0$.]

In order to measure the deviation from GR implied by our numerical star model, we have extracted several observable, gauge-invariant characteristics of our solution. First, we extracted an estimate of the total Keplerian-Schwarzschildian mass parameter $m_{S}$ (as measured faraway) by fitting (in the interval $6<r<10$) the numerical value of $r^{2}F(r)$ to its analytically predicted asymptotic expansion $\sim m_{S}(1+2m_{S}/r)+C_{-}^{F}e^{-\kappa r}(1+\kappa r)+C_{+}^{F}e^{+\kappa r}(1-\kappa r)$. This gave us

[TABLE]

where the digit in parenthesis is a rough measure of the uncertainty (in the last digit) on the numerical determination of $m_{S}$. Note that this is only slightly smaller than what would be the value of the total mass in Einstein’s theory, namely $m_{\rm GR}=e_{0}R_{s}^{3}/12\approx 0.101111$. We have verified that such a value is compatible with our second-order-corrected mass value, $m_{1}+m_{2}$, with $m_{2}$ given by Eq. (180). [It happens that the form factor $\mathcal{E}(z_{s})$, though still negative, is quite small, thereby explaining why one does not see the expected larger self-gravity binding effect due to a high compactness $\sim 0.3$.]

The formally defined compactness $2m_{S}/R_{s}$ would then be $2m_{S}/R_{s}\approx 0.272$. However, such a formal definition (directly copied on GR expressions) does not correspond to any observable characteristics of a star in torsion bigravity. We therefore extracted other (in principle) observable features and numbers from our solution.

We have seen above that if one probes our bigravity field at, say, distances $r\gtrsim 5/\kappa$, the geometry will look like a GR metric of mass $m_{S}$. On the other hand, the exact torsion bigravity metric functions $F=\Phi^{\prime}$ and $\Lambda$ start significantly differing from their GR counterparts $F_{S}(r)\equiv F_{\rm GR}(r,m_{S})$ and $\Lambda_{S}(r)\equiv\Lambda_{\rm GR}(r,m_{S})$ as $r$ gets smaller and comparable to $1/\kappa$. This is illustrated in Figs. 2 and 3. These figures show that, near the star, the torsion bigravity solution differs by $\gtrsim 100\%$ from its GR counterparts.

Let us observationally define the compactness of a star by the surface value of

[TABLE]

In our torsion bigravity model, we can compute the surface value $\Phi(R_{s})$ of $\Phi$ (relative to a zero value at infinity: $\Phi(\infty)=0$) by integrating $F(r)$, namely

[TABLE]

A numerical evaluation of this integral gave us

[TABLE]

and therefore

[TABLE]

Note that this is significantly larger than the corresponding value in GR for a star having the same mass and the same radius, namely

[TABLE]

As a supplementary measure of the strong-gravity nature of our torsion-bigravity star model, let us also cite the value of the geometric-deformation quantity $1-e^{-2\Lambda(R_{s})}$ (which is also equal to the compactness $\frac{2m_{S}}{R_{s}}$ in GR)

[TABLE]

This value confirms that our torsion bigravity model induces large deformations of the geometry.

Another quantity of direct observational significance is the radius of the innermost (or last) stable circular orbit (LSO). From Eq. (4.14) of Ref. Damour:2009sm , the condition defining the LSO reads (in terms of the variables $A\equiv e^{2\Phi}$ and $u\equiv\frac{1}{r}$) ,

[TABLE]

Transcribed in terms of the function $F(r)$, this yields

[TABLE]

Solving this equation gave us

[TABLE]

Note that the ratio $R_{\rm LSO}/m_{S}\approx 15.4$ is about 2.57 larger than the well-known corresponding GR value $R^{\rm GR}_{\rm LSO}/m_{S}=6$. This difference is linked to the fact (already apparent in Fig. 2) that the gravitational field near a torsion bigravity star (of a given Keplerian mass) is significantly more attractive than in GR. [This increase in the strength of the gravitational attraction is essentially due to the extra (short-range) attraction provided by the massive spin-2 excitation.] Note that the value of the ratio $R_{\rm LSO}/m_{S}$ is in principle extractable from the observation of an accretion disk around a neutron star.

In the following section, we shall discuss more potential phenomenological aspects of torsion bigravity.

X Phenomenology of torsion bigravity

We present a preliminary analysis of the phenomenology of torsion bigravity based on the first two orders of perturbation theory, and focussing on solar-system tests of gravity.

X.1 Assuming $r\sim 1/\kappa$

Let us first consider the case where $r\sim 1/\kappa$, i.e. when the exponential decrease of the massive spin-2 excitation is important in the considered physical situation. In that case, torsion bigravity already introduces a modification of Einstein’s (purely massless) theory at the Newtonian level, i.e. when considering the linearized-gravity interaction between two slowly moving massive objects. As already mentioned, previous studies of the linearized approximation Hayashi:1980ir ; Nikiforova:2009qr have shown that the linearized interaction between two massive objects (with stress-energy tensor $T_{\mu\nu}$) involves the exchange of two fields: a massless Einsteinlike gravitational field $h^{*}_{\mu\nu}$, and a massive spin-2 field (contained within the 24 components of the contorsion tensor). The massless field $h^{*}_{\mu\nu}$ couples to $T_{\mu\nu}$ with the Newtonianlike coupling constant

[TABLE]

while the massive spin-2 excitation couples to $T_{\mu\nu}$ with the effective Yukawa-Newtonian coupling constant

[TABLE]

This means that the gravitational interaction term of the source $T_{\mu\nu}$ with itself (after integrating out the field degrees of freedom) reads

[TABLE]

with

[TABLE]

Here the extra numerical prefactors $2$ and $\frac{3}{2}$ are such that the interaction between two nonrelativistic ($T_{\mu\nu}=T_{00}\delta^{0}_{\mu}\delta^{0}_{\nu}$) stationary ($\Box=\Delta$) sources read

[TABLE]

If we consider the interaction between a test particle of mass $M_{2}$ and a spherical object (say a nonrelativistic star) of constant density $e_{0}$ and total mass $M_{1}=\int d^{3}xe_{0}$, separated by a distance $r_{12}$ (between their centers of mass), the above formulas yield an interaction potential $V_{\rm int}=-\int d^{3}xL_{\rm int}$

[TABLE]

Here the form factor $\mathcal{F}(\kappa R_{1})$ (where $R_{1}$ denotes the radius of the object $M_{1}$) is the (normalized) one introduced in Eq. (131). [If we were considering the interaction between two constant-density spherical objects, we should include two form factors: $\mathcal{F}(\kappa R_{1})\mathcal{F}(\kappa R_{2})$. In the case of a test particle considered here, we have $\mathcal{F}(\kappa R_{2})\to 1$.] It is easily checked that the radial force $F_{\rm int}=-\partial V_{\rm int}^{\rm Newtonian}/\partial r_{12}$ deduced from the interaction potential is simply equal to (setting $z_{s}=\kappa R_{1}$)

[TABLE]

where the function $F_{1}(r)$ denotes the external value of our linearized variable $F(r)=\Phi^{\prime}(r)$, as obtained in Eq. (135) above. This is a direct check of the superposition of massless and massive spin-2 excitations in the Newtonianlike potential $\Phi=\frac{1}{2}\ln(-g_{00})$.

There are many experimental data that have set upper limits on the existence, in addition to the Newtonian $1/r$ interaction, of a Yukawa-type interaction $\alpha\,e^{-\kappa r}/r$ coupled with gravitational strength to matter. See Refs Fischbach:1999bc ; Adelberger:2003zx for reviews of the experimental situation. [Note that, when considering non spin-polarized sources, the torsion bigravity interaction respects the equivalence principle, as assumed in the presently considered composition-independent limits.] The Yukawa strength parameter $\alpha$ entering these limits is simply $\alpha=G_{m}/G_{0}=\frac{4}{3}\eta$. The experimental limits on $\alpha$, as a function of $\lambda\equiv 1/\kappa$ are summarized in Fig. 2.13 of Fischbach:1999bc and Fig. 4 of Adelberger:2003zx (for the range $10^{-3}{\rm m}<\lambda<10^{+15}{\rm m}$). We note that the less stringent upper limits apply in the geophysical range (i.e. for $1\,{\rm m}\lesssim\kappa^{-1}\lesssim 10\,{\rm km}$) and roughly limits $\eta=\frac{3}{4}\alpha$, to be

[TABLE]

A range of order $\kappa^{-1}\sim 10\,{\rm km}$ is interesting to consider if one wishes to discuss possible deviations from GR in the physics of neutron stars and black holes.

X.2 PPN parametrization of the second-order torsion bigravity metric when assuming $R_{s}<r\ll 1/\kappa$

Let us now consider the other phenomenological situation where the massive-gravity range is much larger than all the length scales of our system. [We exclude from our consideration the case where $1/\kappa$ is roughly between $10\,{\rm km}$ and $10^{+11}{\rm km}$, for which there are very stringent limits on $\eta$ coming from Earth-satellite, lunar and planetary data.]

If we consider the motion of classical, non-spin-polarized, test masses in our second-order torsion bigravity spacetime (endowed with the metric $g_{\mu\nu}$ and the connection $A_{ij\mu}$), it is given (as shown in Ref. Hayashi:1980av ) by geodesics of the metric $g_{\mu\nu}$. The observational differences (say for the motion of the planets around the Sun) between torsion bigravity and GR are then encapsulated in the difference between our spherically symmetric metric

[TABLE]

and the usual Schwarzschild metric. As is well-known, solar-system experiments are primarily sensitive only to the first post-Newtonian approximation to the metric in the solar system, which is described by the Eddington PPN parameters $\beta$ and $\gamma$. When using (as we do) a Schwarzschildlike radial coordinate, the PPN parameters are defined by writing the first post-Newtonian metric as (see, e.g., Weinberg:1972kfs )

[TABLE]

where $m_{0}=G_{N}M_{0}$ is some observable Keplerian mass parameter. Such an expansion assumes the presence of only power-law deviations from Einstein’s theory. In order to be consistent with it, we shall therefore assume in the present subsection that the Compton wavelength $1/\kappa$ is much larger than the length scales that are being experimentally probed.

The first equation (X.2) implies the following second-order expansions for $\Phi$ and its radial derivative $F=\Phi^{\prime}$:

[TABLE]

and

[TABLE]

Similarly one gets

[TABLE]

Let us now compare these expansions to the corresponding $\kappa\to 0$ limits of the torsion bigravity variables $F=F_{1}+F_{2}$ and $\Lambda_{1}$. According to (164) and (186), we have the following result

[TABLE]

where

[TABLE]

In addition, from Eq. (139) one gets the following $\kappa\to 0$ solution for $\Lambda_{1}$,

[TABLE]

One should identify the observable Keplerian mass $m_{0}$ with the mass parameter $m_{F}$ (which includes self-gravity effects). Then one can conclude from the last equality and Eq. (224) that we can indeed parametrize the linearized torsion-bigravity metric by an Eddington $\gamma$ parameter equal to

[TABLE]

where we consistently neglected the $O(\eta m_{1}/R_{s})$ nonlinear, gravitational binding energy correction term.

The expression (228) for $\gamma$ encapsulates two main facts related to a theory involving both a massless graviton and a massive one. We recall that $\eta$ measures the ratio between the coupling of the massive graviton to that of the massless one, see Eq. (II). When $\eta\to 0$, $\gamma\to 1$, which is the usual Einstein value, while when $\eta\to\infty$, $\gamma\to\frac{1}{2}$, which is the value corresponding to pure massive gravity Boulware:1973my .

There are stringent limits on the deviation $\gamma-1$ between the PPN parameter $\gamma$ and its Einstein value, see notably Refs. Bertotti:2003rm ; Fienga2015 . Note that the Einstein value $\gamma=1$ is obtained for $\eta\to 0$ and that $\gamma=1-\frac{2\eta}{3}+O(\eta^{2})$ as $\eta\to 0$. Using the limits from Ref. Bertotti:2003rm we see that, in the case where $\kappa^{-1}$ is very large, the allowed upper limit on $\eta$ is of order

[TABLE]

Coming back to the second-order terms in $F$, Eq. (186), we see that there are two types of deviations from Einstein’s theory. First, there is a term parametrizable by the PPN parameter $\beta$ (see (223)) with

[TABLE]

Using the fact that $m_{0}=m_{F}=\left(1+\frac{4}{3}\eta\right)m_{1}$, we get the following value of $\beta$ in torsion bigravity

[TABLE]

Note that $\beta\to 23/32$ as $\eta\to\infty$, while, in the limit $\eta\to 0$ we have

[TABLE]

Therefore the upper limit (229) on $\eta$ suffices to guarantee that $\beta-1\lesssim 10^{-5}$, which is more than sufficient to be compatible with the planetary limits on $\beta-1$ Fienga2015 .

Concerning the remaining second-order contribution $\propto\eta(1+\eta)m_{1}^{2}R_{s}^{2}/r^{5}$ in Eq. (186), we note that it is smaller than the non-Einsteinian term $-2(\beta-\gamma)m_{0}^{2}/r^{3}$ by a factor (when $\eta\to 0$) of order $(R_{s}/r)^{2}$, which is much smaller than 1 in all planetary tests. It can therefore be neglected with respect to the usual PPN terms. One should take it into account only when discussing relativistic-gravity tests for near-Earth satellites.

Let us finally recall that the results of the present section have been deduced from the assumption that the second-order perturbation theory of torsion bigravity yields a sufficiently accurate description of the deviations from GR. In our conclusions, we will discuss what modifications might exist if higher-order terms in the perturbation expansion introduce new features in the $\kappa\to 0$ limit.

XI Conclusions

We studied the spherically symmetric (and static) sector of torsion bigravity theories, i.e. the four-parameter class of Einstein-Cartan-type theories (with dynamical torsion) that contain only two physical excitations (around flat spacetime): a massless spin-2 excitation and a massive spin-2 one (of mass $\kappa$). We found that this sector of torsion bigravity has the same number of degrees of freedom (as counted by the total differential order of the equations, after discounting algebraic identities) as their analogs in ghost-free bimetric gravity theories, defined à la DeRham-Gabadadze-Tolley-Hassan-Rosen (see Eqs. (VI)). Knowing that, by contrast, spherically symmetric solutions in generic (ghost-full) bimetric gravity theories exhibit one more degree of freedom (corresponding to the Boulware-Deser ghost), this finding suggests that torsion bigravity might preserve its good ($2+5$) number of degrees of freedom in the full nonlinear regime.

Another remarkable feature of torsion bigravity concerns its behavior in the limit where the mass of the spin-2 excitation tends to zero ($\kappa\to 0$). Contrary to what happens in all bimetric gravity theories (where ordinary perturbation theory is marred by the presence of powers of $\kappa^{-2}$ that increase at each order of perturbation theory, see, e.g., Eqs. (80), (185)), we found that the perturbation theory (around flat space) of torsion bigravity involves no powers of $\kappa^{-2}$ at the first two orders of perturbation theory.

We numerically constructed a high-compactness ($|g_{00}+1|_{\rm surface}=0.45$) (asymptotically flat) star model in torsion bigravity and showed that its physical properties are significantly different from those of a general relativistic star having the same observable Kepler-Schwarzschild mass. See, e.g., Eqs. (206), (207) and equations around. We emphasized that, contrary to the Einstein-Cartan theory (where the torsion does not propagate), the dynamical torsion present in torsion bigravity is generated by the stress-energy tensor $T_{\mu\nu}$ of matter (even in absence of a spin-density distribution) and can lead (when $\eta=1$) to significant differences between the Levi-Civita connection, and the torsionfull one. See Fig. 1.

We also briefly discussed (in section X) possible phenomenologies of torsion bigravity (depending on the considered range $\kappa^{-1}$ of the massive excitation, and on the value of the ratio $\eta$ between the coupling $G_{m}$ of the massive graviton to that, $G_{0}$, of the massless one). As we are not assuming in this work that an analog of the Vainshtein mechanism might be at work in torsion bigravity, we relied on the fact that the physical effects of torsion (for non spin-polarized bodies) disappear in the $\eta\to 0$ limit to give upper limits on $\eta$ making torsion bigravity compatible with existing solar-system tests of GR. We leave to future work an analysis of the compatibility of torsion bigravity with other tests of GR, notably in binary-pulsar data and gravitational-wave data.

As already mentioned, remarkable cancellations of $1/\kappa^{2}$ factors take place at the first two orders of the perturbation theory of torsion bigravity. If these cancellations continued at all orders, one could use torsion bigravity to define an infrared modification of gravity and consider its cosmological applications (as was already attempted in previous work). On the other hand, if $1/\kappa^{2}$ factors arise at the third order of perturbation theory, a preliminary analysis suggests that they could generate contributions to the gravitational acceleration field $F=\Phi^{\prime}$ (with $g_{00}=-e^{2\Phi}$) of the type

[TABLE]

Compared to the first-order result $F_{1}\sim\frac{m}{r^{2}}$ this would mean that perturbation theory might lose its validity below a Vainshteinlike radius which could either be

[TABLE]

If we wished to consider a range $1/\kappa$ of cosmological magnitude, both possibilities would be problematic for the phenomenological consequences we deduced above from second-order perturbation theory. This would then raise the issue of whether a Vainshteinlike mechanism might be at work in torsion bigravity. We leave to future work a discussion of this issue, which is expected to be quite different from the discussion of the $\kappa^{2}\to 0$ limit in ordinary Fierz-Pauli-type massive gravity models because $\kappa^{2}$ enters the torsion bigravity action directly as a denominator (via $c_{F^{2}}=\frac{\eta\,{\lambda}}{\kappa^{2}}$), while Fierz-Pauli-type actions contain $\kappa^{2}$ in the numerator.

We wish, however, to recall that the issue of an eventual bad behavior in the $\kappa^{2}\to 0$ limit is separate from the issue of absence of a sixth degree of freedom, and of ghost-freeness, in the nonlinear regime. In addition, it is only relevant if one wishes to consider a range $1/\kappa$ of cosmological magnitude. We are currently more interested in considering ranges of relevance for modifying the gravitational interaction of compact objects (neutron stars or black holes).

Our hope is that torsion bigravity might define a theoretically healthy alternative to GR that could lead to an interesting modified phenomenology for the gravitational-wave physics of coalescing binary systems of black holes or neutron stars. The present work is just a first step in this programme. In particular, we have shown the existence of high-compacness star models. In the present work, we have only exhibited one model based on an unrealistic constant-density equation of state, but we have also constructed neutron-star models based on more realistic nuclear equations of state (with a range $\kappa^{-1}\sim 10\;{\rm km}$).

We have also noted that the exterior Schwarzschild solution defines a black hole solution in torsion bigravity. We leave to future work the issue of whether this is the unique spherically-symmetric black hole solution of torsion bigravity, or whether there exist black holes with torsion hair. Our hope is that the different Young tableau description of the massive gravity field might allow for black-hole hair.

We leave also to future work an Hamiltonian analysis of torsion bigravity to examine whether its good linearization properties around simple backgrounds, together with the good degree-of-freedom count in fully nonlinear static-spherically-symmetric solutions, are sufficient to ensure ghost-freeness (and mathematical well-posedness) in the full nonlinear theory.

Appendix A Reminders on the Einstein-Cartan formalism

In this Appendix we recall some of the basic technicalities of the Einstein-Cartan(-Weyl-Sciama-Kibble) formalism (also called Poincaré gauge theory). We generally follow the notation of Hayashi:1979wj ; Hayashi:1980av ; Hayashi:1980ir ; Hayashi:1980qp , and of the papers Nair:2008yh ; Nikiforova:2009qr ; Nikiforova:2016ngy ; Nikiforova:2017saf ; Nikiforova:2017xww ; Nikiforova:2018pdk , except for the notation used for the parameters entering into the action. We use a mostly plus signature and distinguish Lorentz-frame indices ($i,j,k,\ldots=0,1,2,3$) from coordinate ones $\mu,\nu,\ldots=0,1,2,3$. The co-frame (inverse of the vierbein) is denoted ${e^{i}}_{\mu}$ (i.e. $g_{\mu\nu}\equiv\eta_{ij}{e^{i}}_{\mu}{e^{j}}_{\nu}$), while the independent (but metric-preserving) SO$(3,1)$ connection is denoted ${A^{i}}_{j\mu}$. These fields respectively define the one-forms $e^{i}=e^{i}_{\mu}dx^{\mu}$ and ${{{\cal A}}^{i}}_{j}={A^{i}}_{j\mu}dx^{\mu}$. In turn, the basic Cartan formulas defining the (torsionless) Levi-Civita connection ${\omega^{i}}_{j}\equiv{\omega^{i}}_{j\mu}dx^{\mu}$ (often called the Riemannian spin-connection), the Riemann curvature of $e^{i}$, the torsion two-form, and curvature two-form of ${{{\cal A}}^{i}}_{j}$ respectively read:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The frame components ${T^{i}}_{[jk]}=-{T^{i}}_{[kj]}$ of the torsion tensor can be written as

[TABLE]

where $C_{i[jk]}=-C_{i[kj]}$ are the structure constants of the vierbein, defined as

[TABLE]

Here, frame indices $i,j,k$ are moved by $\eta_{ij}$.

The explicit links between the contorsion tensor $K_{ijk}=-K_{jik}$ (defined as $K_{ijk}\equiv A_{ijk}-\omega_{ijk}$) and the torsion tensor are

[TABLE]

[TABLE]

Let us also mention the expression of the Riemannian spin-connection in terms of the vierbein and its derivatives

[TABLE]

The frame components of the two curvature tensors, namely ${R^{i}}_{jkl}\equiv{R^{i}}_{j\mu\nu}{e_{k}}^{\mu}{e_{l}}^{\nu}$ and ${F^{i}}_{jkl}\equiv{F^{i}}_{j\mu\nu}{e_{k}}^{\mu}{e_{l}}^{\nu}$, can then be explicitly written (in their “all indices down” forms: $R_{ijkl}\equiv\eta_{ii^{\prime}}{R^{i^{\prime}}}_{jkl}$ and $F_{ijkl}\equiv\eta_{ii^{\prime}}{F^{i^{\prime}}}_{jkl}$) as

[TABLE]

[TABLE]

The tensor and scalar curvatures with contracted indices are defined as follows,

[TABLE]

[TABLE]

Appendix B Link with the notation used in Nikiforova:2018pdk

In our previous paper Ref. Nikiforova:2018pdk we considered the more general action

[TABLE]

containing a cosmological constant $c_{2}$ and five coupling constants $\widetilde{\alpha},\overline{\alpha},c_{3},c_{4},c_{5},c_{6}$. In order to avoid pathologies around flat spacetime, these parameters must satisfy the equation

[TABLE]

and the inequalities

[TABLE]

The field content of such a model around flat space consists of a massless spin-2, a massive spin-2 and a massive pseudoscalar field. The corresponding masses are Nikiforova:2009qr

[TABLE]

while that of the pseudoscalar field is

[TABLE]

We define torsion bigravity by setting $c_{6}=0$ so as to “freeze out” the pseudoscalar field (which becomes infinitely massive). We also set for simplicity the bare cosmological constant $c_{2}$ to zero. This leaves us with only four independent parameters: ${\widetilde{\alpha}}$, ${\overline{\alpha}}$, $c_{3}$ and $c_{4}$.

We then find convenient to change the notation of the parameters and to introduce

[TABLE]

In terms of these parameters, and of the symmetric ($F_{(ij)}$) and antisymmetric ($F_{[ij]}$) parts of $F_{ij}=F_{(ij)}+F_{[ij]}$, the torsion bigravity Lagrangian density reads

[TABLE]

This model contains only a massless spin-2 and a massive one of squared mass

[TABLE]

Acknowledgments

The authors thank Cedric Deffayet, Jürg Fröhlich, Patrick Iglesias-Zemmour, Sergiu Klainerman and Mikhail Volkov for useful discussions.