Unusual square roots in the ghost-free theory of massive gravity

TL;DR

This paper investigates the mathematical intricacies of the square root matrix in ghost-free massive gravity, proposing a new eigenvalue-based formulation to address issues of existence, uniqueness, and perturbation theory.

Contribution

It introduces an eigenvalue-based approach to massive gravity that overcomes limitations of matrix perturbation theory related to square root ambiguities.

Findings

Discrete and continuous parts of square root freedom identified

Eigenvalue formulation enables meaningful perturbation analysis

Addresses non-analytic features in the theory

Abstract

A crucial building block of the ghost free massive gravity is the square root function of a matrix. This is a problematic entity from the viewpoint of existence and uniqueness properties. We accurately describe the freedom of choosing a square root of a (non-degenerate) matrix. It has discrete and (in special cases) continuous parts. When continuous freedom is present, the usual perturbation theory in terms of matrices can be critically ill defined for some choices of the square root. We consider the new formulation of massive and bimetric gravity which deals directly with eigenvalues (in disguise of elementary symmetric polynomials) instead of matrices. It allows for a meaningful discussion of perturbation theory in such cases, even though certain non-analytic features arise.

Equations

$$\begin{gathered} S=-\frac{M_g^2}{2}\int d^{4}x\sqrt{-g}{R}_g-\frac{M_f^2}{2}\int d^{4}x\sqrt{-f}{R}_f \\ +m^2{M}_g^2\int d^{4}x\sqrt{-g}\sum _{n=0}^4{\beta }_n{e}_n\left(\sqrt{g^{-1}f}\right)+\int d^{4}x\sqrt{-g}{\mathcal{L}}_m(g,\Phi ) \end{gathered}$$

$$e_n(\mathcal{M})=\frac{1}{n}\sum _{i=1}^n(-1{)}^{i-1}[{\mathcal{M}}^i]\cdot e_{n-i}(\mathcal{M})$$

$$e_1(\sqrt{\mathbb{X}})$$

$$e_2(\sqrt{\mathbb{X}})$$

$$e_3(\sqrt{\mathbb{X}})$$

$$e_4(\sqrt{\mathbb{X}})$$

$$e_4(\mathcal{M})=\det \left(\mathcal{M}\right)$$

$$g\leftrightarrow f;{\beta }_n\leftrightarrow {\beta }_{4-n};M_g\leftrightarrow M_f;m^2\leftrightarrow m^2\frac{M_g^2}{M_f^2}.$$

$$\delta\left[\Big{(}\sqrt{g^{-1}f}\Big{)}^{n}\right]=\frac{n}{2}\cdot\left[g\Big{(}\sqrt{g^{-1}f}\Big{)}^{n}\delta g^{-1}\right]$$

$$\frac{2}{\sqrt{-g}}\delta\left(\sqrt{-g}e_{n}\left(\sqrt{g^{-1}f}\right)\right)=\sum^{n}_{m=0}(-1)^{m+1}\left[g\Big{(}\sqrt{g^{-1}f}\Big{)}^{m}\delta g^{-1}\right]\cdot e_{n-m}\left(\sqrt{g^{-1}f}\right).$$

$$\stackrel{g}{R_{\mu \nu }}-\frac{1}{2}\stackrel{g}{R}{g}_{\mu \nu }+\frac{m^2}{2}\sum _{n=0}^3(-1{)}^n{\beta }_n\left(g_{\mu \lambda }{{\mathbb{Y}}_{(n)}}_{\nu }^{\lambda }\left(\sqrt{g^{-1}f}\right)+g_{\nu \lambda }{{\mathbb{Y}}_{(n)}}_{\mu }^{\lambda }\left(\sqrt{g^{-1}f}\right)\right)=\frac{1}{M_g^2}{T}_{\nu }^{\mu }$$

$${\mathbb{Y}}_{(0)}\left(\sqrt{\mathbb{X}}\right)=\mathbb{I},$$

$${\mathbb{Y}}_{(1)}\left(\sqrt{\mathbb{X}}\right)=\sqrt{\mathbb{X}}-\mathbb{I}\cdot \left[\sqrt{\mathbb{X}}\right],$$

$${\mathbb{Y}}_{(2)}\left(\sqrt{\mathbb{X}}\right)={\sqrt{\mathbb{X}}}^2-\sqrt{\mathbb{X}}\cdot \left[\sqrt{\mathbb{X}}\right]+\frac{1}{2}\mathbb{I}\cdot \left({\left[\sqrt{\mathbb{X}}\right]}^2-\left[{\sqrt{\mathbb{X}}}^2\right]\right),$$

$$\begin{gathered} {\mathbb{Y}}_{(3)}\left(\sqrt{\mathbb{X}}\right)={\sqrt{\mathbb{X}}}^3-{\sqrt{\mathbb{X}}}^2\cdot \left[\sqrt{\mathbb{X}}\right]+\frac{1}{2}\sqrt{\mathbb{X}}\cdot \left({\left[\sqrt{\mathbb{X}}\right]}^2-\left[{\sqrt{\mathbb{X}}}^2\right]\right) \\ -\frac{1}{6}\mathbb{I}\cdot \left({\left[\sqrt{\mathbb{X}}\right]}^3-3\left[\sqrt{\mathbb{X}}\right]\cdot \left[{\sqrt{\mathbb{X}}}^2\right]+2\left[{\sqrt{\mathbb{X}}}^3\right]\right), \end{gathered}$$

$${\mathbb{Y}}_{(n)}\left(\sqrt{\mathbb{X}}\right)=\sum _{i=0}^n(-1{)}^i{e}_i\left(\sqrt{\mathbb{X}}\right)\cdot {\left(\sqrt{\mathbb{X}}\right)}^{n-i}.$$

$$\stackrel{f}{R_{\mu \nu }}-\frac{1}{2}\stackrel{f}{R}{f}_{\mu \nu }+\frac{m^2}{2M_{\ast }^2}\sum _{n=0}^3(-1{)}^n{\beta }_{4-n}\left(f_{\mu \lambda }{{\mathbb{Y}}_{(n)}}_{\nu }^{\lambda }\left(\sqrt{f^{-1}g}\right)+f_{\mu \lambda }{{\mathbb{Y}}_{(n)}}_{\mu }^{\lambda }\left(\sqrt{f^{-1}g}\right)\right)=0$$

$${\stackrel{g}{\nabla }}_{\mu }\sum _{n=0}^3(-1{)}^n{\beta }_n\left({{\mathbb{Y}}_{(n)}}_{\nu }^{\mu }\left(\sqrt{g^{-1}f}\right)+g^{\mu \alpha }{g}_{\nu \lambda }{{\mathbb{Y}}_{(n)}}_{\alpha }^{\lambda }\left(\sqrt{g^{-1}f}\right)\right)=0.$$

$${\stackrel{f}{\nabla }}_{\mu }\sum _{n=0}^3(-1{)}^n{\beta }_{4-n}\left({{\mathbb{Y}}_{(n)}}_{\nu }^{\mu }\left(\sqrt{f^{-1}g}\right)+f^{\mu \alpha }{f}_{\nu \lambda }{{\mathbb{Y}}_{(n)}}_{\alpha }^{\lambda }\left(\sqrt{f^{-1}g}\right)\right)=0.$$

$$\mathbb{X}\equiv g^{\mu \alpha }{\eta }_{\alpha \nu }=\frac{1}{N^2}\left(\begin{array}{cc}1& N^l{\delta }_{lj}\\ -N^j& (N^2{\gamma }^{il}-N^i{N}^l){\delta }_{lj}\end{array}\right)$$

$${\mathcal{L}}_{{\beta }_1}^{\text{dRGT}}={\pi }^{ij}{\partial }_t{\gamma }_{ij}+N{R}^0+N^i{R}_i+2m^2\sqrt{\gamma }N\left(3-\left[\sqrt{\mathbb{X}}\right]\right)$$

$$N^i=({\delta }_j^i+N{D}_j^i(n,\gamma ))\cdot n^j$$

$$N\sqrt{\mathbb{X}}=\mathbb{A}+N\mathbb{B}.$$

$${\mathbb{A}}_k^i=\frac{1}{\sqrt{1-n^r{\delta }_{rs}{n}^s}}\left(\begin{array}{cc}1& n^r{\delta }_{rk}\\ -n^i& -n^i{n}^r{\delta }_{rk}\end{array}\right);\mathbb{B}=\left(\begin{array}{cc}0& 0\\ 0& \sqrt{({\gamma }^{is}-D_r^i{n}^r{n}^l{D}_l^s){\delta }_{sk}}\end{array}\right);$$

$$\left(\sqrt{1-n^r{\delta }_{rs}{n}^s}\right)D_k^i=\sqrt{\left({\gamma }^{is}-D_r^i{n}^r{n}^l{D}_l^s\right){\delta }_{sk}}$$

$$\left(R_i+\frac{2m^2\sqrt{\gamma }{n}^l{\delta }_{li}}{\sqrt{1-n^r{\delta }_{rs}{n}^s}}\right)\cdot \left({\delta }_k^i+N\frac{\partial }{\partial n^k}(D_j^i{n}^j)\right)=0$$

$$n^i=-R_j{\delta }^{ji}{\left[4m^4\text{det}\gamma +R_k{\delta }^{kl}{R}_l\right]}^{-1/2}.$$

$$R^0+R_i{D}_j^i{n}^j+2m^2\sqrt{\gamma }\left[3-\sqrt{1-n^r{\delta }_{rs}{n}^s}{D}_k^k\right]=0$$

$$\mathbb{A}\cdot \mathbb{A}=\mathbb{X}$$

$${\left(\mathcal{S}\mathbb{A}{\mathcal{S}}^{-1}\right)}^2=\mathcal{S}{\mathbb{A}}^2{\mathcal{S}}^{-1}$$

$$\mathbb{X}={\mathcal{S}}^{-1}\cdot \mathrm{diag}({\lambda }_1,{\lambda }_2,\dots {\lambda }_n)\cdot \mathcal{S}$$

$$\sqrt{\mathbb{X}}={\mathcal{S}}^{-1}\cdot \mathrm{diag}\left(\pm \sqrt{{\lambda }_1},\pm \sqrt{{\lambda }_2},\dots \pm \sqrt{{\lambda }_n}\right)\cdot \mathcal{S}$$

Full text

1 Introduction

Over the last century the theory of general relativity has been confronted with experiment at ever increasing precision in a big variety of scales. It works very well in many different circumstances, from Solar System tests and even GPS technology to cosmological perturbation theory. However, in the latter case there is a number of issues that might hint towards modifications of GR. Those include the missing mass and the accelerated expansion problems. It perfectly motivates a whole variety of approaches to modify or extend the theory of gravitational interactions. Massive (and bimetric) gravity represent a rather popular option for that; one which is theoretically challenging but promising and relatively well studied.

The idea of a massive graviton dates back to classical work of Fierz and Pauli in 1939 [1] which at the time was more of a theoretical excercise. At later times, the naive expectation was that the Yukawa exponent should give rise to an effective cut-off for the gravitational force on large scales leading to possible explanation of the Dark Energy phenomenon. In reality it turned out not that simple, of course. And, even at the most formal level, one of the earliest issues was the infamous vDVZ discontinuity (van Dam-Veltman-Zakharov [2, 3]) discovered in 1970: when the graviton mass tends to zero, the (by definition, linearised) Fierz-Pauli theory does not restore massless GR since the scalar mode does not completely decouple. This problem can have a resolution via non-linear effects known under the name of Vainshtein mechanism [4] which has been proposed in 1972, see also Ref. [5].

However, the second significant problem was precisely about building a non-linear realisation of the Fierz-Pauli model. In GR ten independent variables are reduced to two degrees of freedom because the lapse and shifts play the role of four Lagrange multipliers enforcing four first class constraints in the spatial sector of four otherwise independent variables. If one modifies the theory, the non-dynamical variables generically would enter the action non-linearly leading to disappearance of constraints in the spatial sector and to six degrees of freedom being left. Five of them are the usual polarisations of massive spin two, while the remaining one is a ghost, the Boulware-Deser ghost [6]. Therefore, a modified theory needs to have a constraint for elimination of the ghost. A beautiful solution was suggested in 2010 by de Rham, Gabadadze and Tolley [7, 8]. It has a pair of second class constraints for the spatial sector [9, 10].

Obviously, such models require a second metric on top of the physical one (in order to have an invariant non-derivative interaction potential, in the first place). Initially, the model made use of Minkowski metric $\eta_{\mu\nu}$ which serves as a proper background for perturbation theory in the weak gravity limit. However, it was quickly generalised to an arbitrary non-dynamical symmetric tensor $f_{\mu\nu}$ [11, 12], and even to the full-fledged bimetric gravity [12, 13] with two dynamical metrics $g_{\mu\nu}$ and $f_{\mu\nu}$. Massive gravity and its various extensions have been widely applied in cosmological settings. Having a viable cosmology with massive gravity appears non-trivial [14] but not impossible [15, 16]. Actually, it might even be good enough for reproducing effects of Dark Matter [17]. Nice reviews of massive and bimetric gravity can be found in Refs. [18, 19].

As we have already stated above, massive gravity has two metrics, one of which ($g_{\mu\nu}$) is the physical one which couples to matter, while the second one is an auxiliary (fiducial) metric which can either be fixed of dynamical. They go with an interaction potential made of $\sqrt{g^{-1}f}$, the square root of the matrix ${\mathbb{X}}^{\mu}_{\nu}\equiv g^{\mu\alpha}f_{\alpha\nu}$. Among very many papers devoted to different aspects of bimetric gravity, only very few [20, 21, 22, 23, 24] deal with the formal aspects of solving the matrix equation $\mathbb{A}^{2}=\mathbb{X}$ that arises in theory and defines the interaction between metrics. This paper aims at filling this yawning gap based on our new approach [25] which deals with eigenvalues instead of matrix components.

The paper is organised as follows. In Section 2 we review the bimetric gravity. In Section 3 (and Appendix A) we describe square roots of matrices and related problems of massive gravity. In Section 4 we describe our new formulation of the model. In Section 5 we apply it to perturbation theory around unusual square roots with continuous ambiguity. Finally, in Section 6 we discuss implications and conclude.

2 Bimetric theory

Let us first review the basic foundations of bimetric theory. The action can be written (in four dimensions) as

[TABLE]

with the elementary symmetric polynomials $e_{n}\equiv\sum\limits_{i_{1}<i_{2}<\ldots<i_{n}}\lambda_{i_{1}}\lambda_{i_{2}}\cdots\lambda_{i_{n}}$ of the eigenvalues $\lambda_{i}$ of the matrix $\sqrt{g^{-1}f}\equiv\sqrt{\mathbb{X}}$. By definition we take $e_{0}(\mathcal{M})=1$, and then the following recurrent relation can be proven:

[TABLE]

for any matrix $\mathcal{M}$ where $\left[\mathcal{M}\right]$ denotes the trace of $\mathcal{M}$. In particular,

[TABLE]

and in four dimensions we automatically have

[TABLE]

and $e_{n}=0$ for $n>4$. The field $\Phi$ stands in the action (1) for unspecified matter fields which couple to the physical metric $g_{\mu\nu}$. If needed, the number of spacetime dimensions can be changed straightforwardly.

The mass parameter $m$ is redundant with parameters $\beta_{i}$, where $\beta_{0}$ and $\beta_{4}$ (or $\beta_{D}$ in $D$ dimensions) give rise only to cosmological constants for metrics $g_{\mu\nu}$ and $f_{\mu\nu}$ respectively while three ($D-1$) remaining parameters correspond to three ($D-1$) different interaction terms. The action is invariant under the following simultaneous transformations:

[TABLE]

In order to derive the equations of motion, let us notice that from $\delta\mathbb{X}=\delta\left(\sqrt{\mathbb{X}}\cdot\sqrt{\mathbb{X}}\right)$ follows $\delta\mathbb{X}=\left(\delta\sqrt{\mathbb{X}}\right)\cdot\sqrt{\mathbb{X}}+\sqrt{\mathbb{X}}\cdot\left(\delta\sqrt{\mathbb{X}}\right)$, or equivalently $\sqrt{\mathbb{X}}^{-1}\cdot\delta\mathbb{X}=\sqrt{\mathbb{X}}^{-1}\cdot\left(\delta\sqrt{\mathbb{X}}\right)\cdot\sqrt{\mathbb{X}}+\delta\sqrt{\mathbb{X}}$. Taking the trace of this equation we get $\left[\delta\sqrt{\mathbb{X}}\right]=\dfrac{1}{2}\left[\sqrt{\mathbb{X}}^{-1}\cdot\delta\mathbb{X}\right]$. More generally, we have the following identity [9]:

[TABLE]

which gives

[TABLE]

Varying the action (1) with respect to $g_{\mu\nu}$ we obtain by virtue of (4):

[TABLE]

where $\mathbb{Y}_{(n)}$ are defined as:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

or in general

[TABLE]

Analogously, varying the action (1) with respect to $f_{\mu\nu}$, one can get

[TABLE]

where the dimensionless ratio of two Planck masses has been introduced: $M_{*}=\frac{M_{f}}{M_{g}}$. The quantities with letter $g$ (respectively $f$) on top are defined with respect to the metric $g_{\mu\nu}$ (respectively $f_{\mu\nu}$).

For completeness, let us mention that, as a consequence of the Bianchi identity and the covariant conservation of $T_{\mu\nu}$, the $g_{\mu\nu}$ equation of motion (5) leads to the Bianchi constraint

[TABLE]

Similarly for $f_{\mu\nu}$ we find from eq. (7) that

[TABLE]

In this paper we consider the massive gravity regime with the auxiliary metric $f_{\mu\nu}=\eta_{\mu\nu}$, and therefore the $R_{f}$ term in the action should be ignored. However, the issues we discuss are generic for any situations in massive and bimetric gravity in metric formulation.

2.1 The constraint analysis

We are interested in how far the good properties of the model, including the absence of Boulware-Deser ghost, can go with unusual square roots. Initially, the degrees of freedom count has been done in perturbation theory around the trivial square root of unit matrix. However, it left open the question of self-consistency beyond perturbative level which required methods of non-perturbative Hamiltonian analysis.

The latter was first done by Hassan and Rosen in a series of papers starting from [10] with the simplest case of $\beta_{1}$ model and Minkowski fiducial metric which is enough for the purposes of our discussion. They used the ADM decomposition for $g_{\mu\nu}$:

[TABLE]

where the lapse $N$, the shift $N_{i}$, and the spatial metric $\gamma_{ij}$ are the usual ADM variables for $g_{\mu\nu}$, and $\delta_{ij}$ is the Euclidean spatial metric (as a part of fiducial Minkowski $\eta_{\mu\nu}$).

The first order Lagrangian of the $\beta_{1}$ model can be written in these variables as:

[TABLE]

where the number $3$ in the potential stands solely for compensating the effective cosmological constant from the $\beta_{1}$ term. $R^{0}$, $R_{i}$ are the usual constraints of GR, and $\pi^{ij}$ is the canonical conjugate momentum of $\gamma_{ij}$. The crucial feature of generic massive gravity is that the lapse and shift enter the potential term non-linearly, and therefore they fail to serve as Lagrange multipliers. Their variations give algebraic equations for themselves instead of constraining the spatial sector. That’s how we get $6$ degrees of freedom: $5$ of massive spin two, and $1$ of the scalar ghost.

A way out would be to have a very special potential which allows one (primary) constraint to survive in the spatial sector. Hassan and Rosen have shown that one can make a (linear in lapse) redefinition of the shift vector

[TABLE]

such that the Lagrangian density of the $\beta_{1}$ model is rendered linear in the lapse,

[TABLE]

It allows the lapse variable to serve as Lagrange multiplier and constrain the spatial sector while the new shift $n^{i}$ would be determined by its own (algebraic) equation of motion.

Namely, one can check by direct calculation [10] that the matrices $\mathbb{A}$ and $\mathbb{B}$ are given by

[TABLE]

provided that the transformation of shifts is obtained by solving equation

[TABLE]

The equation of motion for $n^{i}$ reads

[TABLE]

where the second factor is the Jacobian of the transformation (8) and cannot be equal to zero, while vanishing of the first factor determines the shift:

[TABLE]

The variation with respect to $N$ gives then:

[TABLE]

which, after substitution of the shifts, is the constraint equation for the spatial sector. Its conservation in time would amount to another constraint of a second class pair. At the final step of the Hamiltonian analysis the value of the lapse would be determined.

Let us note that this proof allows for at least some freedom in choosing the square root $\mathbb{X}$. Indeed, the $3\times 3$-dimensional square root remains unspecified. And the freedom of changing the sign of the temporal eigenvalue can be realised by simply changing the overall sign of the $\sqrt{\mathbb{X}}$ together with explicit three-dimensional freedom. However, we see that changing the signs of $\mathbb{A}$ and $\mathbb{B}$ does not alter the conclusions more than simply changing the sign of the shift: $n^{i}=+R_{j}\delta^{ji}\left[4m^{4}\text{det}\gamma+R_{k}\delta^{kl}R_{l}\right]^{-1/2}.$ Therefore, it seems that the proof of Hassan and Rosen has pretty much generic consequences, beyond the standard choice of the square roots.

One can also give independent arguments which look even more generic [28]. It is possible to write down the $\beta_{1}$ term as $\beta_{1}\sqrt{\gamma}\left[\mathbb{F}\right]$ with the constraint, enforced by Lagrange multipliers, that ${\mathbb{F}}^{2}=N\mathbb{X}$. Then one can show that there appears a constraint in the spatial sector [28]. Therefore it might be sensible to consider massive gravity with non-standard square roots. And we must understand the mathematical structures behind them.

3 The structure of square roots

The mathematics of extracting square roots of matrices amounts to solving a simple matrix equation

[TABLE]

where $\mathbb{X}$ is a given $n\times n$ matrix ($g^{-1}f$ in our case), and $\mathbb{A}$ is an unknown which we denote by $\mathbb{A}\equiv\sqrt{\mathbb{X}}$.

Note that the procedure of squaring a matrix is respected by similarity transformations in the sense that

[TABLE]

for any non-singular matrix $\mathcal{S}$. As usual, it is more convenient to start the discussion over the field of complex numbers. As is well known [27], any matrix can be brought to the Jordan normal form by a similarity transformation. Therefore, it is enough to understand how to take square roots of matrices in the Jordan normal form, for otherwise a suitable similarity transformation can be applied.

It can be shown, see Appendix A and Ref. [27], that a square root of a Jordan normal form is obtained by taking square roots of separate Jordan blocks111This is the only possibility if the eigenvalues of different Jordan blocks are not equal to each other. Otherwise there is also an additional continuous freedom in choosing a square root, see below. which is always possible (recall that we are working with complex numbers) apart from the case of degenerate ($\lambda=0$) Jordan blocks of dimension greater than one. The latter case does not concern us much since our matrices are by definition non-degenerate (see however [21, 22]). The details are given in the Appendix A, while here we assume for simplicity that the matrix $\mathbb{X}$ is diagonalisable, i.e. all its Jordan blocks are $1$-dimensional.

For any diagonalisable matrix

[TABLE]

we can take

[TABLE]

which features discrete freedom which corresponds to $2^{n}$ elements in $\left(\mathbb{Z}_{2}\right)^{n}$.

If some eigenvalues are equal then also a larger freedom appears. Indeed, if we have chosen $+\sqrt{\lambda}$ in one place and $-\sqrt{\lambda}$ in another, then a non-trivial similarity transformation in the subspace spanned by corresponding eigenvectors would change the square root $\sqrt{\mathbb{X}}$ but not the matrix $\mathbb{X}$ itself, for the latter is proportional to identity in this subspace.

Suppose that the continuous freedom is present. With a given choice of eigenvalues for $\sqrt{\mathbb{X}}$, the set of different solutions for $\sqrt{\mathbb{X}}$ naturally has a structure of smooth manifold. Indeed, the group of similarity transformations can be viewed as the adjoint representation of the general linear group $GL(n,\mathbb{C})$ given by ${\mathfrak{L}}_{\mathcal{S}}:\ {\mathbb{X}}\longrightarrow{\mathcal{S}}{\mathbb{X}}{\mathcal{S}}^{-1}$. The stabiliser subgroup (or centraliser) of a given matrix $\mathbb{X}$ is, by definition, a subgroup which maps $\mathbb{X}$ to itself, ${\mathfrak{S}}_{\mathbb{X}}\equiv\{{\mathcal{S}}\in GL(n,\mathbb{C}):\ {\mathfrak{L}}_{\mathcal{S}}\mathbb{X}=\mathbb{X}\}$. Obviously, ${\mathfrak{S}}_{\sqrt{\mathbb{X}}}\subset{\mathfrak{S}}_{\mathbb{X}}$, and the continuous freedom arises when ${\mathfrak{S}}_{\sqrt{\mathbb{X}}}\neq{\mathfrak{S}}_{\mathbb{X}}$ which is true if and only if equal eigenvalues of $\mathbb{X}$ were chosen to have unequal square roots. In this case ${\mathfrak{S}}_{\sqrt{\mathbb{X}}}$ naturally acts on ${\mathfrak{S}}_{\mathbb{X}}$, and the freedom of choosing the square root is parametrised by the homogeneous space [26] ${\mathfrak{S}}_{\mathbb{X}}/{\mathfrak{S}}_{\sqrt{\mathbb{X}}}$.

3.1 On conditions of reality

Of course, we are interested in classifying the real square roots. In this case, one has to bear in mind that a real matrix can have a pair of complex eigenvalues, complex conjugate to each other. As an example, note that a matrix $\begin{pmatrix}a&b\\ -b&a\\ \end{pmatrix}$ has the normal form which reads $\begin{pmatrix}a+ib&0\\ 0&a-ib\\ \end{pmatrix}$. In this case we need to choose the square roots $\sqrt{a\pm ib}$ to be complex conjugate to each other if we want the inverse similarity transformation to bring it back to real values. It narrows the freedom.

Note also that if the matrix $\mathbb{X}$ has negative eigenvalues, then their square roots are imaginary. However, if they go in pairs, or more generally if there are even numbers of identical Jordan blocks with negative eigenvalues, then the $\sqrt{\lambda}$-s can be chosen in complex conjugate pairs (which of course invokes continuous freedom), and the matrix $\sqrt{\mathbb{X}}$ would become real after a suitable similarity transformation from what is allowed by the continuous freedom. Otherwise, a real square root does not exist. See also the Ref. [29] for geometrical and physical understanding of reality conditions.

3.2 The problem with variations

In the standard approach to massive gravity we make an assumption that the object $\delta\sqrt{\mathbb{X}}\equiv\sqrt{\mathbb{X}+\delta\mathbb{X}}-\sqrt{\mathbb{X}}$ does exist as an infinitesimal variation of $\sqrt{\mathbb{X}}$. This assumption can break down with non-standard square roots [23]. This is somewhat unusual. Indeed, usually we can even have a series expansion in powers of $\sqrt{\mathbb{X}}$, at least around a solution with proportional metrics for which $\mathbb{X}=c\cdot\mathbb{I}$ with $c\in\mathbb{R}$. However, even for proportional backgrounds, if we choose $\sqrt{\mathbb{I}}\neq\mathbb{I}$ then a generic $\delta\mathbb{X}$ does not commute with $\sqrt{\mathbb{I}}$, and we cannot write the usual expansion for $\sqrt{\mathbb{I}+\delta\mathbb{X}}$ in terms of powers of $\delta\mathbb{X}$.

It makes the standard variations problematic. And actually, in certain cases, the problem is more profound than just a mere technicality. In cases of continuous freedom a smooth variation of $\sqrt{\mathbb{X}}$ does not exist at all. Let us consider the $2\times 2$ unit matrix $\mathbb{I}=\left(\begin{matrix}1&0\\ 0&1\end{matrix}\right)$. It has two isolated roots, $\sqrt{\mathbb{I}}_{++}=\mathbb{I}=\left(\begin{matrix}1&0\\ 0&1\end{matrix}\right)$ and $\sqrt{\mathbb{I}}_{--}=-\mathbb{I}=\left(\begin{matrix}-1&0\\ 0&-1\end{matrix}\right)$, and also a two-dimensional manifold of roots $\sqrt{\mathbb{I}}_{ab}=\mathbb{I}=\left(\begin{matrix}a&b_{1}\\ b_{2}&-a\end{matrix}\right)$ where $a^{2}=1-b_{1}b_{2}$. However, even an infinitesimal change of the initial matrix leaves only two points from this manifold222Introducing new variables $b_{\pm}\equiv\frac{1}{\sqrt{2}}(b_{1}\pm b_{2})$, we see that the manifold is the surface $a^{2}+b_{+}^{2}=1+b_{-}^{2}$ in three-dimensional space.. Indeed, for $\mathbb{M}_{\epsilon}=\left(\begin{matrix}1&0\\ 0&1+\epsilon\end{matrix}\right)$ we have only four isolated roots: $\sqrt{\mathbb{M}_{\epsilon}}_{++}=\left(\begin{matrix}1&0\\ 0&\sqrt{1+\epsilon}\end{matrix}\right)$,  $\sqrt{\mathbb{M}_{\epsilon}}_{--}=\left(\begin{matrix}-1&0\\ 0&-\sqrt{1+\epsilon}\end{matrix}\right)$,  $\sqrt{\mathbb{M}_{\epsilon}}_{+-}=\left(\begin{matrix}1&0\\ 0&-\sqrt{1+\epsilon}\end{matrix}\right)$, and $\sqrt{\mathbb{M}_{\epsilon}}_{-+}=\left(\begin{matrix}-1&0\\ 0&\sqrt{1+\epsilon}\end{matrix}\right)$.

Let us recapitulate the situation. If the freedom is only of discrete type, then it is basically a technical issue of finding a proper series expansion when $\mathbb{X}$ and $\delta\mathbb{X}$ do not commute [30]. For example, the linear correction solves the equation of Sylvester type: $\sqrt{\mathbb{X}}\cdot\delta\sqrt{\mathbb{X}}+\delta\sqrt{\mathbb{X}}\cdot\sqrt{\mathbb{X}}=\delta\mathbb{X}$. And the problem is to find $\delta\sqrt{\mathbb{X}}$ explicitly. It can be done at the expense of using clever field redefinitions [31].

Continuous freedom presents a much worse problem. Typically in symmetric set-ups, we can have coincident eigenvalues in the background matrix $\mathbb{X}$, and therefore the freedom of continuous type in $\sqrt{\mathbb{X}}$. Generic perturbation $\delta\mathbb{X}$ would split the eigenvalues, and only some discrete subset of possible solutions for $\sqrt{\mathbb{X}}$ would be smoothly related to those for $\sqrt{\mathbb{X}+\delta\mathbb{X}}$. Perturbation theory in terms of matrices is not well defined. Note that the approach with the Sylvester equation works well if and only if the spectra of $\sqrt{\mathbb{X}}$ and $-\sqrt{\mathbb{X}}$ have empty intersection [31, 32]. Obviously, it fails precisely when some equal eigenvalues have different square roots and continuous freedom exists.

Note though that the continuous freedom does not affect the eigenvalues which are the only thing which enters the action. Therefore, it should be better to perform the variations at the level of eigenvalues. We have recently proposed one such method [25]. In this approach we are not interested in square root matrices, instead we define the polynomials $e_{i}$ in the action (1) as solutions of algebraic equations in terms of elementary symmetric polynomials of the $g^{-1}f$ or $f^{-1}g$ matrix. Effectively, it eliminates the continuous freedom since the eigenvalues are not changing under the similarity transformations.

4 New approach to the action principle

As we have seen, working in terms of eigenvalues could allow to overcome the problems with perturbation theory which arise when continuous freedom is present. However, explicitly computing eigenvalues for a generic matrix in massive gravity might be a very unpleasant task. In our previous paper [25] we offered a better method. The action (1) depends only on elementary symmetric polynomials of eigenvalues. And we used the following observation:

[TABLE]

to find how the symmetric polynomials of $\mathbb{X}$ and $\sqrt{\mathbb{X}}$ are related to each other.

From now on we denote for brevity

[TABLE]

and get the desired relation from (11) with $\mathcal{M}=\sqrt{\mathbb{X}}$:

[TABLE]

In four dimension we have more explicitly:

[TABLE]

The interaction Lagrangian takes the form of $\sum\beta_{i}e_{i}$. And we can perform the variation as:

[TABLE]

where one easily gets

[TABLE]

Indeed, for $n=1$ we trivially see that $g^{\mu\alpha}\frac{\partial}{\partial g^{\alpha\nu}}\left[g^{-1}\eta\right]={\mathbb{X}}^{\mu}_{\nu}$ and more generally $g^{-1}\frac{\partial}{\partial g^{-1}}\left[\mathbb{X}^{n}\right]=n{\mathbb{X}}^{n}$, and then go by induction using $p_{n}=\frac{1}{n}\sum\limits_{i=1}^{n}(-1)^{i-1}[\mathbb{X}^{i}]\cdot p_{n-i}$.

In any case when the standard approach in terms of matrices works well, it is guaranteed that these equations of motion coincide with the usual ones since we are doing literally the same variation though without having the matrix elements explicitly. Let us illustrate how it works in the simplest case of two dimensions where equations (12) take the form of

[TABLE]

and therefore

[TABLE]

After a simple calculation we obtain from (17) and (18)

[TABLE]

Multiplying this relation by $g^{-1}$, we need to prove that ${\mathcal{M}}^{2}=\mathbb{X}$ for

[TABLE]

where we have used the Cayley-Hamilton theorem (see Appendix A.5), $\mathbb{X}^{2}-p_{1}\mathbb{X}+p_{2}\mathbb{I}=\mathbb{O}$, for the matrix $\mathbb{X}$. It can be done by direct calculation employing again the Cayley-Hamilton theorem.

In four dimensions the things become rather bulky, and we have

[TABLE]

where

[TABLE]

and for the $\beta_{1}$ term we have to check again that ${\mathcal{M}}^{2}=\mathbb{X}$ with

[TABLE]

where we have used, and need to further use in ${\mathcal{M}}^{2}$, the Cayley–Hamilton theorem which in four dimensions reads $\mathbb{X}^{4}-p_{1}\mathbb{X}^{3}+p_{2}\mathbb{X}^{2}-p_{3}\mathbb{X}+p_{4}\mathbb{I}=\mathbb{O}$. After we have got an expression for $\sqrt{\mathbb{X}}$, all other terms can be considered straightforwardly.

Obviously, it can be done in any dimension, and works well as long as ${\rm det}\left(\frac{\partial p}{\partial e}\right)\neq 0$. However, we will not go into this combinatorics, and we only write down the determinant in the three dimensional case

[TABLE]

for future use.

4.1 Perturbation theory

Obviously, it is not always possible to have an exact solution, and one often uses perturbation theory. In particular, for the weak gravity limit the background solution corresponds to $\mathbb{X}=\mathbb{I}$ and values of polynomials $p_{1}=4,\quad p_{2}=6,\quad p_{3}=4,\quad p_{4}=1$. With the trivial choice of the square root $\sqrt{\mathbb{X}}=\mathbb{I}$ we have $e_{i}=p_{i}$, and one can check that ${\rm det}\left(\frac{\partial p}{\partial e}\right)\neq 0$.

For perturbative calculations around Minkowski space in massive gravity it is convenient to use the symmetry of massive gravity (3) and work with $\eta^{-1}g\equiv\mathbb{I}+\mathcal{H}$ matrix instead of $g^{-1}\eta$ since $\sqrt{-g}e_{i}\left(\sqrt{g^{-1}\eta}\right)=e_{4-i}\left(\sqrt{\eta^{-1}g}\right)$, and having worked out a symmetric polynomial one does not need to multiply its perturbative expansion by expansion for $\sqrt{-g}$.

In the previous paper [25], we have shown how it goes, and reproduced the standard result that Minkowski space is a vacuum solution if and only if

[TABLE]

and the quadratic action is then

[TABLE]

which is the Fierz-Pauli term.

5 Perturbation theory around unusual roots

Now let us consider perturbation theory around unusual roots. Note that in those cases when odd numbers of eigenvalues are taken with reversed signs, we need to be more accurate about the symmetry transformation (3). Since any such choice entails ${\rm det}\left(\sqrt{\eta^{-1}g}\right)<0$, then we need to take $\beta_{i}\leftrightarrow-\beta_{D-i}$ instead of $\beta_{i}\leftrightarrow\beta_{D-i}$ if $\sqrt{-g}$ is kept positive. This is not a big deal, but unfortunately it is not the end of the story.

Unusual roots are precisely the cases when ${\rm det}\left(\frac{\partial p}{\partial e}\right)=0$. And perturbation theory ceases to be analytic and nice. Let us first look at a toy example of $2D$. Unusual square roots are given by

[TABLE]

and all its similarity transformations. In other words, we have $e_{1}=0,\quad e_{2}=-1$ in the background. Then we have $p_{1}=2+e_{1}(\mathcal{H})$ and $p_{2}=1+e_{1}(\mathcal{H})+e_{2}(\mathcal{H})$ with perturbation $\mathcal{H}\equiv\eta^{-1}\cdot\delta g$. Equations for polynomials with our background choice can be solved exactly as

[TABLE]

We see that there are two solutions for $e_{1}$. Moreover, at the lowest order we have

[TABLE]

which scales linearly with $\mathcal{H}$ but not analytic.

The reason is simple. For a generic perturbation $\mathcal{H}$ two eigenvalues are different, and the sign of $e_{1}$ depends on which of them (the smaller or the larger) is taken with the minus sign.

In the interaction potential we have then ${\mathcal{L}}_{int}=-\beta_{0}e_{2}-\beta_{1}e_{1}$ which has at the linear level $e_{2}=-1-\frac{1}{2}{\mathcal{H}}^{\mu}_{\mu}+{\mathcal{O}}\left({\mathcal{H}}^{2}\right)$ and $e_{1}=\pm\sqrt{\frac{1}{2}{\mathcal{H}}^{\mu\nu}{\mathcal{H}}_{\mu\nu}-\frac{1}{4}\left({\mathcal{H}}^{\mu}_{\mu}\right)^{2}}\cdot\left(1+{\mathcal{O}}({\mathcal{H}})\right)$. We see that vanishing of the linear variation with arbitrary ${\mathcal{H}}$ around this solution requires two conditions, for analytic and non-analytic parts separately. In $2D$ it is too restrictive, and sets both $\beta$-s to zero.

Obviously, non-analytic features should be seen in more complicated cases, too. The number of solutions is given by the number of possible ways to choose the distribution of signs among the square roots of eigenvalues. We will now check it explicitly in three and four dimensions.

5.1 The 3D case

Let us now turn to three dimensions. Since we are interested in Minkowski background, we need to discuss square roots of three-dimensional unit matrix. In terms of eigenvalues, we choose the signs of three square roots of unity. Barring the overall sign of the potential term, there are two distinct options: either all three square roots have equal signs which is the standard solution without continuous freedom, or we have one eigenvalue with a different sign. We see that in three dimensions it suffices to study the following unusual square root of $\mathbb{I}$:

[TABLE]

and all its similarity transformations.

We have the interaction term

[TABLE]

where the prefactor of $m^{2}M_{Pl}^{2}$ is omitted (or absorbed into the $\beta$-s). Also the $\beta_{3}$ term has been omitted since it does not depend on $g$. The difference from the standard approach is that we now treat the quantities $e_{i}$ as solutions of algebraic equations (12) with no regard to square root matrices which effectively eliminates the continuous freedom of similarity transformations. Equations (12) take the form

[TABLE]

with $p_{i}\equiv e_{i}\left(\eta^{-1}g\right)=e_{i}\left(\mathbb{I}+\mathcal{H}\right)$:

[TABLE]

where $h_{i}\equiv e_{i}(\mathcal{H})$ are elementary symmetric polynomials of the matrix ${\mathcal{H}}^{\mu}_{\nu}=\eta^{\mu\alpha}\cdot\delta g_{\alpha\nu}$.

We now build the perturbation theory in the form of

[TABLE]

with $e^{(i)}_{k}$ being the $i$-th order (which scales as the $i$-th power of $\mathcal{H}$) correction to the $k$-th polynomial. The [math]-th orders are by definition the background values: $e^{(0)}_{1}=1$, $e^{(0)}_{2}=-1$, $e^{(0)}_{3}=-1$ for which the determinant in the formula (20) is obviously equal to zero. Nevertheless, the third polynomial can be directly solved for according to (24):

[TABLE]

Now we substitute the expansions (25) and (26) into the system of equations (22), (23) and keep only terms of the first order:

[TABLE]

The system is obviously degenerate. We have

[TABLE]

with one first order variable undetermined at the linear level.

Taking this (27) into account, we get again two coincident equations (22), (23) at quadratic order:

[TABLE]

which can be solved as

[TABLE]

now leaving one second order variable undetermined.

Let us now turn to the third order. The first equation (22) acquires the form

[TABLE]

and, after taking it into account, the second one (23) finally allows to know the first order completely:

[TABLE]

This is a generic cubic equations (the three coefficients depend on three parameters $h_{i}$). Its three roots correspond to freedom of choosing which one of three eigenvalues of perturbed matrix $\mathbb{X}$ would go with the reversed sign.

It is worth to perform a consistency check here. Let $\lambda_{i}$, $i=1,2,3$ be eigenvalues of $\mathcal{H}$. Then we have $h_{1}=\sum\limits_{i}\lambda_{i}$, $h_{2}=\sum\limits_{i<j}\lambda_{i}\lambda_{j}$, $h_{3}=\sum\limits_{i<j<k}\lambda_{i}\lambda_{j}\lambda_{k}=\lambda_{1}\lambda_{2}\lambda_{3}$. One can easily check by Viet theorem that the three solutions of the cubic equation (29) are $\frac{1}{2}\left(-\lambda_{1}+\lambda_{2}+\lambda_{3}\right)$, $\frac{1}{2}\left(-\lambda_{2}+\lambda_{3}+\lambda_{1}\right)$, and $\frac{1}{2}\left(-\lambda_{3}+\lambda_{1}+\lambda_{2}\right)$. These are indeed the three possible first order contributions to $e_{1}$ with eigenvalues $\pm\sqrt{1+\lambda_{i}}$ and our choice of the square root.

After we have reproduced the full freedom of choosing the eigenvalues, the system (22), (23) becomes non-degenerate giving two new pieces of information at each order. In particular, at the fourth order the first equation (22) gives

[TABLE]

At the same time, the second equation (23) contains $e^{(2)}_{1}$ both quadratically and linearly. However, after taking the first equation (30) into account, one gets a linear equation for $e^{(2)}_{1}$ which can be unequivocally solved. Explicit expressions are rather cumbersome, and we don’t go into that in this paper333To avoid confusion, let us mention that there would be a denominator which can take the zero value. However, playing with $\lambda$-s one can check that it happens only when the enumerator is zero, too..

The first order contribution to the interaction Lagrangian is:

[TABLE]

If Minkowski space is to be a solution of equations of motion, then we need to impose $\delta^{(1)}\mathcal{L}_{\rm int}=0$. Since $e^{(1)}_{1}$ is generically irrational for rational values of $h_{i}$, we need the two parts, with $h_{1}$ and with $e^{(1)}_{1}$, to vanish independently:

[TABLE]

It is interesting to note that one can formally plug the Minkowski space with our choice (21) of $\sqrt{g^{-1}\eta}$ into the naive equation of motion (5),

[TABLE]

and deduce the same condition (31) for $\beta$-s.

Now we look at the second order Lagrangian for perturbations:

[TABLE]

We see that it contains the root $e^{(1)}_{1}$ of the cubic equation (29), and therefore is not analytic in $\mathcal{H}$. If we enumerate the eigenvalues of $\mathcal{H}$ such that $e^{(1)}_{1}=\dfrac{1}{2}(-\lambda_{1}+\lambda_{2}+\lambda_{3})$, we get:

[TABLE]

The three branches of solution correspond to three possible choices of two out of three eigenvalues for putting them into $\delta^{(2)}\mathcal{L}_{\rm int}$.

5.2 Square roots in the physical dimension

In four dimensions the system of equations (12) would be slightly more complicated:

[TABLE]

Now we basically have two possible choices of unusual roots $\sqrt{\mathbb{I}}$: with one reversed-sign eigenvalue or with two of them (and again with continuous freedom of similarity transformations which can be disregarded when working directly with eigenvalues).

Let us start with the first option

[TABLE]

which goes pretty much analogous to the three dimensional case. The background values for this choice are $e^{(0)}_{1}=2$, $e^{(0)}_{2}=0$, $e^{(0)}_{3}=-2$, $e^{(0)}_{4}=-1$, and one easily checks that the Jacobian (19) of transformation from $e_{i}$ to $p_{i}$ vanishes.

As before, the last polynomial (35) is easy to find:

[TABLE]

The remaining equations (32) - (34) give at the first order:

[TABLE]

The system is degenerate, and we can only find that

[TABLE]

leaving one variable undetermined.

At the second order we get a similar degeneracy, and find

[TABLE]

The degeneracy gets broken at the fourth order level producing a generic (in terms of independent parameters) quartic equation for $e^{(1)}_{1}$.

Let us now have a look at the interaction Lagrangian:

[TABLE]

the first order contribution to which reads

[TABLE]

Again, the Minkowski space is a solution when $\delta^{(1)}\mathcal{L}_{\rm int}=0$. And we require the two parts, with $h_{1}$ and $e^{(1)}_{1}$, to vanish independently:

[TABLE]

As in $3D$, one can formally plug the Minkowski space with (36) into the naive equations of motion $\sum\limits_{i=0}^{3}(-1)^{i}\beta_{i}\mathbb{Y}_{(i)}=0$:

[TABLE]

and deduce the same condition (37) for $\beta$-s.

Now we look at the second order Lagrangian for perturbations:

[TABLE]

If we enumerate the eigenvalues of $\mathcal{H}$ such that $e^{(1)}_{1}=\dfrac{1}{2}(-\lambda_{1}+\lambda_{2}+\lambda_{3}+\lambda_{4})$, then we get:

[TABLE]

The structure is basically the same as for the Fierz-Pauli term, though with only a part of the full set of eigenvalues. We remind to the reader that for such square roots the usual perturbation theory in terms of matrices is not well defined unless the perturbation respects separation into two subspaces (of positive and negative eigenvalues). In our case, of this and the previous subsections, the latter occurs if $h_{1i}=0$ for $i\geq 2$. Then the structure of quadratic Lagrangian would be $\sum\limits_{i,j\geq 2}h_{ij}h^{ij}-\left(\sum\limits_{i\geq 2}h^{i}_{i}\right)^{2}$. Otherwise the perturbation theory is not analytic and has three (respectively four, or $D$) branches according to solutions of the cubic (respectively quartic, or power $D$) equation.

An interesting fact is that the background equation of motion is correctly reproduced by naively applying formula (5) which is not even well defined in such cases. Apparently, it has to do with the fact that $\mathbb{X}\ \propto\ \mathbb{I}$ can be treated as a limiting case of a cosmological matrix ${\rm diag}\left(N(t),a(t),\ldots,a(t)\right)$ for which the perturbations are well defined. Note though that this limit is singular in the sense that it features appearance of non-analyticity in perturbations, and is not defined at all in the language of matrices (compare to [23] where they also mention non-analyticity in the case of non-zero ”new shift” variable which is probably the first indication in the literature towards the phenomena we have thoroughly studied in the present paper).

5.3 4D: the second case

Another option available in four dimensions can be represented by the root

[TABLE]

for which the background values are $e^{(0)}_{1}=0$, $e^{(0)}_{2}=-2$, $e^{(0)}_{3}=0$, $e^{(0)}_{4}=1$; they correspond to vanishing determinant in equation (19). And the last polynomial (16) is

[TABLE]

The first order equations (13) - (15) are particularly simple and degenerate:

[TABLE]

which only give

[TABLE]

leaving two other variables undetermined.

At the second order the equations (13) - (15) are

[TABLE]

and we deduce

[TABLE]

At the third order, the first equation (13) gives

[TABLE]

and after that the second (14) and the third (15) equations are degenerate with each other and give an interesting result:

[TABLE]

At this point one might start worrying about the self-consistency. Wouldn’t the perturbation theory break down at the second order when $e^{(1)}_{1}=0$? Indeed, given some eigenvalues $\lambda_{i}$, there should be six values of $e^{(1)}_{1}$ of the form $\frac{1}{2}\left(\lambda_{1}+\lambda_{2}-\lambda_{3}-\lambda_{4}\right)$. Can a third order symmetric polynomial $h_{3}-\frac{1}{2}h_{1}h_{2}+\frac{1}{8}h_{1}^{3}$ be divisible by all six of them? Actually, it can since those six values come in three pairs of opposite sign numbers. And one can easily check that

[TABLE]

which gives a finite value for

[TABLE]

Then, at the fourth order, the system of equations (13) - (15) becomes non-degenerate, and the first equation (13) reads

[TABLE]

Two other equations are very cumbersome. However, subtracting the second (14) from the third (15), and using the first one (41), we get a very nice equation which fully determines the first order of perturbations:

[TABLE]

Together with the result (39) for $e^{(2)}_{1}+e^{(2)}_{3}$, it amounts to cubic equation for $\left(e^{(1)}_{1}\right)^{2}$ precisely as it should be.

Moreover, taking again $e^{(1)}_{1}=\frac{1}{2}\left(\lambda_{1}+\lambda_{2}-\lambda_{3}-\lambda_{4}\right)$ we easily see that

[TABLE]

and, given also (40), those values of $e^{(1)}_{1}$ are indeed the solutions of (42).

Let us now have a look at the interaction Lagrangian:

[TABLE]

with the first order expression being

[TABLE]

We look for Minkowski space to be a solution: $\delta^{(1)}\mathcal{L}_{\rm int}=0$, and demand

[TABLE]

Probably, at this point the reader would not be much surprised if we say that the naive equation of motion $\sum\limits_{i=0}^{3}(-1)^{i}\beta_{i}\mathbb{Y}_{(i)}=0$ formally gives

[TABLE]

which entails the same condition (44) for $\beta$-s.

The second order Lagrangian for perturbations is

[TABLE]

which, given the identifications (40) and (43) above, is equivalent to

[TABLE]

The structure is basically the same as for the Fierz-Pauli term, though with two separate parts of the full set of eigenvalues. As we have already mentioned before, the usual perturbation theory in terms of matrices is not well defined unless the perturbation respects separation into two subspaces (of positive and negative eigenvalues). In our case the latter occurs if $h_{ij}=0$ for $i\leq 2$ and $j\geq 3$ or vice versa. Then the structure of quadratic Lagrangian has two Fierz-Pauli-like terms $\sum\limits_{i,j\leq 2}h_{ij}h^{ij}-\left(\sum\limits_{i\leq 2}h^{i}_{i}\right)^{2}$ and$\sum\limits_{i,j\geq 3}h_{ij}h^{ij}-\left(\sum\limits_{i\geq 3}h^{i}_{i}\right)^{2}$. Otherwise the perturbation theory is not analytic. The six branches of it correspond to different separations of four eigenvalues into two groups of two.

Note that this kind of unusual root can also be related to bimetric cosmology. However, making it regular would require at least an axially symmetric Bianchi I spacetime. In bimetric Friedmann cosmology given by

[TABLE]

with $\mathbb{X}=g^{-1}f=\{W^{2},Y^{2}/a^{2},Y^{2}/a^{2},Y^{2}/a^{2}\}$, the classical equations of motion are well defined only for square roots of $\sqrt{g^{-1}f}=\pm\{\pm W,Y/a,Y/a,Y/a\}$ type since they do not make equal eigenvalues different. But we must be careful if choosing another square root: $\sqrt{g^{-1}f}=\{W,Y/a,-Y/a,-Y/a\}$ which would give rise to non-analyticity of the type which has been considered in our Subsection 5.1. The standard perturbation theory would not be well defined (compare to [23]).

In higher dimensions the concrete calculations become more and more involved. We don’t aim at tackling this issues in the present paper. However, it seems reasonable to conjecture that the quadratic Lagrangian around Minkowski would always look like

[TABLE]

where $C_{1}$ and $C_{2}$ are constants, and ${\mathcal{I}}_{-}$ and ${\mathcal{I}}_{+}$ are the sets of indices corresponding to negative and positive eigenvalues of $\sqrt{\mathbb{I}}$ respectively.

6 Discussion and Conclusions

It is amazingly difficult to promote the simple action of Fierz and Pauli to a full-fledged non-linear theory of massive graviton. In the metric formulation one needs to compute the square root function of the matrix $g^{-1}f$. This procedure is definitely not among the simplest recipes of modern theoretical physics. And even in perturbation theory it is not so easy to deal with perturbative expansion of the square root, unless around a matrix proportional to identity. In some cases, it is not a mere technicality but rather a severe obstacle for correctly introducing small perturbations.

There is some freedom in taking the square root of a matrix. It can be of two types. The first type is absolutely generic. It gives discrete freedom which basically refers to the usual two branches of the square root function in complex analysis. The peculiarity is that those branches can be independently chosen for different Jordan blocks. The second type is continuous and happens only in special cases. If we have taken different branches for Jordan blocks with equal eigenvalues, then a manifold of solutions arises given by those similarity transformations which change the square root without affecting the initial matrix.

In presence of continuous freedom, the naive perturbation theory in terms of matrices breaks down. Generic perturbation collapses the manifold to a finite number of points which depend on the chosen perturbation. Perturbation theory around such square roots is not well defined. A positive observation on top of that is that the similarity transformations which generate the manifold do not change the eigenvalues, and the eigenvalues are the only quantities which affect the Lagrangian. Of course, working directly with eigenvalues is very difficult. But fortunately everything can be done in terms of their elementary symmetric polynomials, and the mission is rendered accomplishable.

In this paper we have described the mathematical issues of working with square roots in massive gravity and developed the new perturbation theory in terms of eigenvalues and their elementary symmetric polynomials. The latter appeared to be a suitable tool for perturbation theory, even when around unusual roots. Our method allows to meaningfully define perturbations in cases of continuous freedom, even though with non-analyticity. Note that the origins of non-analyticity are very clear, however its implications remain to be understood, as well as its possible relations to the vielbein formulations of massive gravity.

Acknowledgements

FS thanks Nicolay Gordeev and Yelena Yashina for useful discussions about some mathematical questions related to this paper. AG enjoyed many inspiring discussions about the square roots and other topics of massive gravity with Fawad Hassan and Mikica Kocic. AG is grateful to the Dynasty Foundation and to the Saint Petersburg State University travel grant 11.42.687.2017 for support; and also support from the Russian Science Foundation (project 16-11-10218) is gratefully acknowledged.

Appendix A Some known facts from the theory of matrices

In this Appendix we briefly review some well known facts from the theory of matrices over the field of complex numbers.

A.1 Spectral theorem

In finite-dimensional spaces one can give the full spectral classification of linear operators, without any requirements of self-adjointness or unitarity. In the language of matrices, an $n\times n$ matrix $\mathbb{X}$ can be brought by a similarity transformation to the Jordan normal form which is a block-diagonal form ${\mathcal{S}}{\mathbb{X}}{\mathcal{S}}^{-1}=J_{k_{1}}(\lambda_{1})\bigoplus J_{k_{2}}(\lambda_{2})\bigoplus\ldots$ with $k_{i}$ standing for the dimensionality of the $i$-th block, and the latter is given by

[TABLE]

with the complex number $\lambda$ being the eigenvalue (of geometric multiplicity $1$ and algebraic multiplicity $k$ in this particular block). In components, we have non-zero elements $J_{i,i}=\lambda$ and $J_{i,i+1}=1$. If all Jordan blocks are $1$-dimensional then the matrix is diagonalisable. Note also that a similarity transformation can be viewed as simply changing the basis in the linear space in which the operator acts. The matrix $\mathcal{S}$ maps one basis onto another.

Various proofs of this theorem can be found in numerous textbooks on matrices. Here we only sketch a simple one. Given an $n\times n$ matrix $\mathbb{X}$ over complex numbers, one can always find an eigenvector444Recall that the solvability condition for the eigenvector equation $\mathbb{X}e=\lambda e$ is $P(\lambda)=0$ where $P$ is the characteristic polynomial of $\mathbb{X}$. And the latter always has a solution in complex numbers., and choosing it as the $n$-th element of the basis we make ${\mathbb{X}}_{n,k}=0$ for $k<n$. After that we do the same for the upper left $(n-1)\times(n-1)$ corner of the new matrix, and so on. By induction, we bring $\mathbb{X}$ to an upper triangular form. Obviously, the numbers on the diagonal are eigenvalues, i.e. the roots of the characteristic polynomial of $\mathbb{X}$. By renumbering the elements of the basis, we can group them into sets of equal eigenvalues (if there are such). Let us choose one of the eigenvalues, $\lambda_{1}$. Then, we have $\mathbb{X}-\lambda_{1}{\mathbb{I}}=\left(\begin{matrix}{\tilde{X}}&M\\ {\mathbb{O}}&X_{0}\end{matrix}\right)$ where $\tilde{X}$ is non-degenerate and $X_{0}$ is nilpotent. If we make a similarity transformation by matrix $\left(\begin{matrix}{\mathbb{I}}&C\\ {\mathbb{O}}&\mathbb{I}\end{matrix}\right)$ then the matrix $\mathbb{X}-\lambda_{1}{\mathbb{I}}$, and also the matrix $\mathbb{X}$, becomes block-diagonal if ${\tilde{X}}C-CX_{0}=M$ which can be solved as $C=\sum\limits_{i\geq 0}{\tilde{X}}^{-i-1}MX_{0}^{i}$ which is a well-defined finite sum since $\tilde{X}$ is non-degenerate and $X_{0}$ is nilpotent. By induction, we separate all blocks with different eigenvalues. Now we only have to understand the structure of a matrix $X_{\lambda_{1}}$ with all eigenvalues equal to a given number $\lambda_{1}$. Obviously, it suffices to study the nilpotent matrix $X_{0}=X_{\lambda_{1}}-\lambda_{1}\mathbb{I}$ of a given size $k\times k$. Actually, what we need now is a very generic statement which does not rely on the nice properties of the complex numbers: for any nilpotent linear operator in any (finite-dimensional) linear space one can decompose the vector space into direct sum of cyclic subspaces. Let’s take the minimal number $k_{1}\leq k$ such that $X_{0}^{k_{1}}=0$. Then there is a vector $e_{0}$ for which $X_{0}^{k_{1}-1}e_{0}\neq 0$. This vector together with all vectors $X_{0}^{i}e_{0}$ generates a vector subspace $\mathcal{V}_{0}$ in which the matrix $X_{0}$ acts as a Jordan block with $\lambda=0$. If $k_{1}=k$, we are done. Otherwise, let’s take a suitable555If the word ”suitable” sounds as too sloppy a definition, one might think of the vector $e_{1}$ as one of those which require the maximal power of $\mathbb{X}$ to be sent to $\mathcal{V}_{0}$ (including zero). Or even more formally, one can consider the factor space of the full $\lambda=0$ block over $\mathcal{V}_{0}$, and perform the construction there. Since $\mathcal{V}_{0}$ is an invariant subspace, the action of $\mathbb{X}$ is correctly defined on the factor space, and any representative of the highest vector can serve as $e_{1}$. vector $e_{1}$ from a linear complement to $\mathcal{V}_{0}$. We want to build a new cyclic subspace by vectors $X_{0}^{i}e_{1}$. This strategy might fail if for some $i$ we get $0\neq X_{0}^{i}e_{1}\in V_{0}$. Then it means that the vector $X_{0}^{i}e_{1}$ is a linear combination of vectors $X_{0}^{j}e_{0}$ with $j\geq i$ (for otherwise we have $X_{0}^{k_{1}}e_{1}\neq 0$). Therefore, there is some $\tilde{e}\in\mathcal{V}_{0}$ such that $X_{0}^{i}\tilde{e}=X_{0}^{i}e_{1}$, and one can substitute $e_{1}$ by $e_{1}^{\prime}=e_{1}-\tilde{e}$ to overcome this difficulty. This procedure can be repeated if needed. By induction, we prove that $X_{0}$ can be brought to the form of direct sum of nilpotent Jordan blocks of sizes $k_{i}$ such that $\sum k_{i}=k$. The spectral theorem is proven.

A.2 Analytic functions of matrices

Similarity transformations preserve polynomial relations. If $P$ is a polynomial, then obviously $P({\mathcal{S}}{\mathbb{X}}{\mathcal{S}}^{-1})={\mathcal{S}}P(\mathbb{X}){\mathcal{S}}^{-1}$. Of course, it can be extended to any functions of matrices defined by norm-convergent series. The Jordan normal form gives also a very convenient tool for representing functions of matrices by

[TABLE]

One can easily compute the powers of the Jordan block and see that for an analytic function $f$ which is given by its Taylor series it gives

[TABLE]

It allows to extend the definition to the class of $\mathcal{C}^{k-1}(\mathfrak{U})$ functions where $k$ is the maximal size of the Jordan blocks, and $\mathfrak{U}\in\mathbb{R}$ is a domain which contains the full spectrum of the matrix.

In particular, since $f(x)=\sqrt{x}$ is smooth for $x\neq 0$, we have a definition of the square root of any non-degenerate matrix. Since the square root is a multivalued function, we see that the definition of $\sqrt{\mathbb{X}}$ admits the discrete freedom of choosing the sign of $\sqrt{\lambda}$ in each Jordan block. One has a matrix from the similarity class of $J_{k}(\pm\sqrt{\lambda})$ for $\sqrt{J_{k}(\lambda)}$ where $\lambda\neq 0$. We will see below that it is not the end of the story if we look at square roots as solutions of the matrix equation $\mathbb{A}\cdot\mathbb{A}=\mathbb{X}$.

Before doing that, let us comment on degenerate matrices. One can easily see that

[TABLE]

Therefore, $\left(J_{3}(0)\right)^{2}$ is equivalent to $J_{2}(0)\bigoplus J_{1}(0)$ under renumbering of the basis’ elements. The Jordan block splits into two. Obviously, the splitting is related to the fact that $f^{\prime}(\lambda)=0$. Note that it never occurs if $f^{\prime}(\lambda)\neq 0$.

Given also that $\left(J_{2}(0)\right)^{2}=\mathbb{O}$, we see that a $3\times 3$ matrix $\mathbb{A}$ such that ${\mathbb{A}}^{2}=J_{3}(0)$ does not exist. For some degenerate matrices, square roots do not exist at all, even over the field of complex numbers.

A.3 Commuting matrices

In order to determine all solutions of quadratic equation $\mathbb{A}\cdot\mathbb{A}=\mathbb{X}$, we will need to know the full set of solutions for commutation equation of the form $\mathbb{A}\mathbb{X}=\mathbb{X}\mathbb{A}$.

Let us work in the basis in which the given matrix $\mathbb{X}$ acquires the canonical Jordan form $\mathbb{X}=\bigoplus\limits_{i}J_{k_{i}}(\lambda_{i})$. This is a block diagonal form, and so let us consider the matrix $\mathbb{A}$ in a block form composed of $k_{i}\times k_{j}$ matrices $A_{i,j}$. Then the equation $\mathbb{A}\mathbb{X}=\mathbb{X}\mathbb{A}$ takes the form $A_{i,j}J_{k_{j}}(\lambda_{j})=J_{k_{i}}(\lambda_{i})A_{i,j}$ (no summation).

If $\lambda_{i}\neq\lambda_{j}$ we get $(\lambda_{i}-\lambda_{j})A_{i,j}=J_{k_{i}}(0)A_{i,j}-A_{i,j}J_{k_{j}}(0)$. The matrices $J(0)$ in the right hand side are nilpotent. We can multiply this equation by $\lambda_{i}-\lambda_{j}$ as many times as we want, and have higher powers of nilpotent matrices in the right hand side. Finally, we get $A_{i,j}=0$.

If $\lambda_{i}=\lambda_{j}$, then we have a simple equation $A_{i,j}J_{k_{j}}(\lambda)=J_{k_{i}}(\lambda)A_{i,j}$. The Jordan blocks have very simple structure, and one can easily check that if $k_{i}=k_{j}$ (for example, $i=j$ - diagonal blocks) then $A_{i,j}$ is an upper triangular matrix with components $\left(A_{i,j}\right)_{a,a+b}=C_{b}$ given by arbitrary numbers $C_{b}$ for $b=0,1,\ldots k-1$. If $k_{i}\neq k_{j}$ then the solution is either $A_{i,j}=\left(\begin{matrix}\mathbb{O}&{\tilde{A}}_{i,j}\end{matrix}\right)$ or $A_{i,j}=\left(\begin{matrix}{\tilde{A}}_{i,j}\\ \mathbb{O}\end{matrix}\right)$ where ${\tilde{A}}_{i,j}$ is an upper triangular matrix of the same form as presented above and the size of ${\rm min}(k_{i},k_{j})\times{\rm min}(k_{i},k_{j})$.

If the matrix $\mathbb{X}$ was diagonalisable, then all the blocks are $1\times 1$-dimensional, and we simply have that the matrix element $A_{ij}$ is arbitrary if $\lambda_{i}=\lambda_{j}$, and $A_{ij}=0$ if $\lambda_{i}\neq\lambda_{j}$. It is very easy to understand: in those subspaces where the matrix $\mathbb{X}$ is proportional to identity it commutes with anything. And in general we need the two matrices to be simultaneously diagonalisable which is well known in quantum mechanics as the condition for two quantities to be simultaneously measurable.

A.4 Square roots of matrices

Let us now assume we are given a non-degenerate matrix

[TABLE]

and need to solve the quadratic equation ${\mathbb{A}}^{2}=\mathbb{X}$ for $\mathbb{A}$. The matrix $\mathbb{A}$ can also be decomposed into the direct some of non-degenerate Jordan blocks which do not split under taking the square of the matrix. And therefore there is a one-to-one correspondence between the Jordan blocks of $\mathbb{X}$ and $\mathbb{A}$. And, obviously, instead of $J_{k}(\pm\sqrt{\lambda})$ we can use $\sqrt{J_{k}(\lambda)}$ in the spectral decomposition for $\mathbb{A}$:

[TABLE]

where $\mathcal{T}$ gives a similarity transformation to be found, with the simplest solution being $\mathcal{T}=\mathcal{S}$.

Given that ${\mathbb{A}}^{2}=\mathbb{X}$ we have: ${\mathcal{T}}^{-1}\left(\bigoplus\limits_{i}J_{k_{i}}(\lambda_{i})\right){\mathcal{T}}={\mathcal{S}}^{-1}\left(\bigoplus\limits_{i}J_{k_{i}}(\lambda_{i})\right){\mathcal{S}}$. It easily translates into the commutation property: ${\mathcal{S}}{\mathcal{T}}^{-1}\left(\bigoplus\limits_{i}J_{k_{i}}(\lambda_{i})\right)=\left(\bigoplus\limits_{i}J_{k_{i}}(\lambda_{i})\right){\mathcal{S}}{\mathcal{T}}^{-1}$. One solution is always $\mathcal{T}=\mathcal{S}$, while the general solution is that $\mathcal{T}{\mathcal{S}}^{-1}$ commutes with $\bigoplus\limits_{i}J_{k_{i}}(\lambda_{i})$. In other words, $\mathcal{T}$ is a product of $\mathcal{S}$ and an arbitrary matrix which commutes with $\bigoplus\limits_{i}J_{k_{i}}(\lambda_{i})$. Of course, it gives a new square root only if $\mathcal{T}{\mathcal{S}}^{-1}$ does not commute with $\bigoplus\limits_{i}\sqrt{J_{k_{i}}(\lambda_{i})}$.

Now, we can see where appears a continuous freedom in square roots. If there are two Jordan blocks in $\mathbb{X}$ with equal eigenvalues, $\lambda_{i}=\lambda_{j}$, then a similarity transformation matrix $\mathcal{S}$ which leaves $\mathbb{X}$ invariant can have non-trivial components in the $S_{i,j}$ block which were described above in A.3. At the same time, if we chose to have $\sqrt{\lambda_{i}}\neq\sqrt{\lambda_{j}}$, then such transformation does change the square root $\sqrt{\mathbb{X}}$. This is how we get the continuous freedom corresponding to the manifold of similarity transformations which change $\sqrt{\mathbb{X}}$ leaving $\mathbb{X}$ intact.

A.5 Cayley-Hamilton theorem

Let $P(\lambda)$ be the characteristic polynomial of $\mathbb{X}$. The theorem states that $P(\mathbb{X})=\mathbb{O}$. It trivially follows from our considerations above. Indeed, if there is a Jordan block $J_{k_{i}}(\lambda_{i})$ in the spectral decomposition of $\mathbb{X}$, then it means that $\lambda_{i}$ is a root of $P$ of order not less than $k_{i}$. Therefore, $P(J_{k_{i}}(\lambda_{i}))=\mathbb{O}$ for every block in the spectral representation of $\mathbb{X}$.