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

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.


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 r1s 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 η µν 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 µν [11,12], and even to the full-fledged bimetric gravity [12,13] with two dynamical metrics g µν and f µν . 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 µν ) 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 a g´1f , the square root of the matrix X µ ν " g µα f αν . 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 A 2 " 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.

Bimetric theory
Let us first review the basic foundations of bimetric theory. The action can be written (in four dimensions) as (1) with the elementary symmetric polynomials e n " ř X. By definition we take e 0 pMq " 1, and then the following recurrent relation can be proven: for any matrix M where rMs denotes the trace of M. In particular, and in four dimensions we automatically have e 4 pMq " det pMq and e n " 0 for n ą 4. The field Φ stands in the action (1) for unspecified matter fields which couple to the physical metric g µν . If needed, the number of spacetime dimensions can be changed straightforwardly.
The mass parameter m is redundant with parameters β i , where β 0 and β 4 (or β D in D dimensions) give rise only to cosmological constants for metrics g µν and f µν respectively while three (D´1) remaining parameters correspond to three (D´1) different interaction terms. The action is invariant under the following simultaneous transformations: In order to derive the equations of motion, let us notice that from δX " δ`?X¨?X˘follows Varying the action (1) with respect to g µν we obtain by virtue of (4): or in general Analogously, varying the action (1) with respect to f µν , one can get where the dimensionless ratio of two Planck masses has been introduced: M˚" Mg . The quantities with letter g (respectively f ) on top are defined with respect to the metric g µν (respectively f µν ).
For completeness, let us mention that, as a consequence of the Bianchi identity and the covariant conservation of T µν , the g µν equation of motion (5) leads to the Bianchi constraint Similarly for f µν we find from eq. (7) that In this paper we consider the massive gravity regime with the auxiliary metric f µν " η µν , 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.

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 nonperturbative Hamiltonian analysis.
The latter was first done by Hassan and Rosen in a series of papers starting from [10] with the simplest case of β 1 model and Minkowski fiducial metric which is enough for the purposes of our discussion. They used the ADM decomposition for g µν : here the lapse N , the shift N i , and the spatial metric γ ij are the usual ADM variables for g µν , and δ ij is the Euclidean spatial metric (as a part of fiducial Minkowski η µν ).
The first order Lagrangian of the β 1 model can be written in these variables as: here the number 3 in the potential stands solely for compensating the effective cosmological constant from the β 1 term. R 0 , R i are the usual constraints of GR, and π ij is the canonical conjugate momentum of γ 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 such that the Lagrangian density of the β 1 model is rendered linear in the lapse, 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 A and B are given by 1´n r δ rs n sˆ1 n r δ rḱ n i´ni n r δ rk˙; B "˜0 0 0 b pγ is´Di r n r n l D s l qδ sk¸; provided that the transformation of shifts is obtained by solving equatioń a 1´n r δ rs n s¯D i k " b`γ is´Di r n r n l D s l˘δ sk The equation of motion for n i reads 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: The variation with respect to N gives then: 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 X. Indeed, the 3ˆ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 ?
X together with explicit three-dimensional freedom. However, we see that changing the signs of A and B does not alter the conclusions more than simply changing the sign of the shift: . 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 β 1 term as β 1 ? γ rFs with the constraint, enforced by Lagrange multipliers, that F 2 " N 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 where X is a given nˆn matrix (g´1f in our case), and A is an unknown which we denote by A " ? X. Note that the procedure of squaring a matrix is respected by similarity transformations in the sense that`S AS´1˘2 " SA 2 S´1 for any non-singular matrix 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 blocks 1 which is always possible (recall that we are working with complex numbers) apart from the case of degenerate (λ " 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 X is diagonalisable, i.e. all its Jordan blocks are 1-dimensional.
For any diagonalisable matrix we can take ? X " S´1¨diag´˘aλ 1 ,˘aλ 2 , . . .˘aλ n¯¨S which features discrete freedom which corresponds to 2 n elements in pZ 2 q n . If some eigenvalues are equal then also a larger freedom appears. Indeed, if we have choseǹ ? λ in one place and´?λ in another, then a non-trivial similarity transformation in the subspace spanned by corresponding eigenvectors would change the square root ? X but not the matrix 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 ? X, the set of different solutions for ? 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 GLpn, Cq given by L S : X ÝÑ SXS´1. The stabiliser subgroup (or centraliser) of a given matrix X is, by definition, a subgroup which maps X to itself, S X " tS P GLpn, Cq : L S X " Xu. Obviously, S ? X Ă S X , and the continuous freedom arises when S ? X ‰ S X which is true if and only if equal eigenvalues of X were chosen to have unequal square roots. In this case S ? X naturally acts on S X , and the freedom of choosing the square root is parametrised by the homogeneous space [26] S X {S ? X .

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ˆa b b a˙h as the normal form which readsˆa`i b 0 0 a´ib˙.
In this case we need to choose the square roots ? a˘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 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 ? λ-s can be chosen in complex conjugate pairs (which of course invokes continuous freedom), and the matrix ? 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.

The problem with variations
In the standard approach to massive gravity we make an assumption that the object δ ? X " ?
X`δX´?X does exist as an infinitesimal variation of ? 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 ? X, at least around a solution with proportional metrics for which X " c¨I with c P R. However, even for proportional backgrounds, if we choose ? I ‰ I then a generic δX does not commute with ? I, and we cannot write the usual expansion for ?
I`δX in terms of powers of δ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 ? X does not exist at all. Let us consider the 2ˆ2 unit matrix I "ˆ1 0 0 1˙. It has two isolated roots, ?
and also a two-dimensional manifold of roots ?
here a 2 " 1´b 1 b 2 . However, even an infinitesimal change of the initial matrix leaves only two points from this manifold 2 . Indeed, for M ǫ "ˆ1 0 0 1`ǫ˙w e have only four isolated roots: ? 1`ǫ˙. 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 X and δX do not commute [30]. For example, the linear correction solves the equation of Sylvester type: ?
X¨?X " δX. And the problem is to find δ ? 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 X, and therefore the freedom of continuous type in ? X. Generic perturbation δX would split the eigenvalues, and only some discrete subset of possible solutions for ? X would be smoothly related to those for ? X`δ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 ? X and´?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 2 Introducing new variables b˘" 1 ? 2 pb1˘b2q, we see that the manifold is the surface a 2`b2 " 1`b 2 in three-dimensional space. 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´1f or f´1g 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: (11) to find how the symmetric polynomials of X and ?
X are related to each other. From now on we denote for brevity Xq " e i and get the desired relation from (11) with M " ? X: p´1q k e k e m " p´1q n p n .
In four dimension we have more explicitly: The interaction Lagrangian takes the form of ř β i e i . And we can perform the variation as: ( 17) where one easily gets Indeed, for n " 1 we trivially see that g µα B Bg αν " g´1η ‰ " X µ ν and more generally g´1 B Bg´1 rX n s " nX n , and then go by induction using p n " 1 n n ř i"1 p´1q i´1 rX i s¨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 e 2 1´2 e 2 " p 1 , e 2 2 " p 2 .
and therefore After a simple calculation we obtain from (17) and (18) B Bg µν`?´g e 1˘" Multiplying this relation by g´1, we need to prove that M 2 " X for here we have used the Cayley-Hamilton theorem (see Appendix A.5), X 2´p 1 X`p 2 I " O, for the matrix 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 and for the β 1 term we have to check again that M 2 " X with M " e 2 1 e 3 X`e 3 1 X 2´e 1 pe 2 e 4 I`pe 2 2`e 4 qX`2e 2 X 2`X3 q`e 3 pe 4 I`e 2 X`X 2 q e 1 e 2 e 3`e 2 3`e 2 where we have used, and need to further use in M 2 , the Cayley-Hamilton theorem which in four dimensions reads X 4´p 1 X 3`p 2 X 2´p 3 X`p 4 I " O. After we have got an expression for ?
X, all other terms can be considered straightforwardly.
Obviously, it can be done in any dimension, and works well as long as det´B p Be¯‰ 0. However, we will not go into this combinatorics, and we only write down the determinant in the three dimensional case for future use.

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 X " I and values of polynomials p 1 " 4, p 2 " 6, p 3 " 4, p 4 " 1. With the trivial choice of the square root ?
X " I we have e i " p i , and one can check that det´B p Be¯‰ 0.
For perturbative calculations around Minkowski space in massive gravity it is convenient to use the symmetry of massive gravity (3) and work with η´1g " I`H matrix instead of g´1η since ?´g e i´a g´1η¯" e 4´i´a η´1g¯, and having worked out a symmetric polynomial one does not need to multiply its perturbative expansion by expansion for ?´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 β 0 "´3β 1´3 β 2´β3 , and the quadratic action is then which is the Fierz-Pauli term.

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 det´aη´1g¯ă 0, then we need to take β i Ø´β D´i instead of β i Ø β D´i if ?´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 det´B p Be¯" 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 ? I "ˆ´1 0 0 1ȧ nd all its similarity transformations. In other words, we have e 1 " 0, e 2 "´1 in the background. Then we have p 1 " 2`e 1 pHq and p 2 " 1`e 1 pHq`e 2 pHq with perturbation H " η´1¨δg. Equations for polynomials with our background choice can be solved exactly as e 2 "´?p 2 "´a1`e 1 pHq`e 2 pHq, e 1 "˘ap 1`2 e 2 "˘b2`e 1 pHq´2 a 1`e 1 pHq`e 2 pHq.
We see that there are two solutions for e 1 . Moreover, at the lowest order we have hich scales linearly with H but not analytic. The reason is simple. For a generic perturbation 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 L int "´β 0 e 2´β1 e 1 which has at the linear level 2˘a nd e 1 "˘b 1 2 H µν H µν´1 4 pH µ µ q 2¨p 1`OpHqq. We see that vanishing of the linear variation with arbitrary H around this solution requires two conditions, for analytic and non-analytic parts separately. In 2D it is too restrictive, and sets both β-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.

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 I: 0 0 0 1 0 0 0 1‚ (21) and all its similarity transformations.
We have the interaction term where the prefactor of m 2 M 2 P l is omitted (or absorbed into the β-s). Also the β 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 with p i " e i`η´1 g˘" e i pI`Hq: where h i " e i pHq are elementary symmetric polynomials of the matrix H µ ν " η µα¨δ g αν . We now build the perturbation theory in the form of  (20) is obviously equal to zero. Nevertheless, the third polynomial can be directly solved for according to (24): Now we substitute the expansions (25) and (26) into the system of equations (22), (23) and keep only terms of the first order: The system is obviously degenerate. We have with one first order variable undetermined at the linear level. Taking this (27)  and, after taking it into account, the second one (23) finally allows to know the first order completely:´e 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 X would go with the reversed sign. It is worth to perform a consistency check here. Let λ i , i " 1, 2, 3 be eigenvalues of H.
One can easily check by Viet theorem that the three solutions of the cubic equation (29) are 1 2 p´λ 1`λ2`λ3 q, 1 2 p´λ 2`λ3`λ1 q, and 1 2 p´λ 3`λ1`λ2 q. These are indeed the three possible first order contributions to e 1 with eigenvalues˘?1`λ 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 At the same time, the second equation (23) contains e p2q 1 both quadratically and linearly. However, after taking the first equation (30) into account, one gets a linear equation for e p2q 1 which can be unequivocally solved. Explicit expressions are rather cumbersome, and we don't go into that in this paper 3 .
The first order contribution to the interaction Lagrangian is: If Minkowski space is to be a solution of equations of motion, then we need to impose δ p1q L int " 0. Since e p1q 1 is generically irrational for rational values of h i , we need the two parts, with h 1 and with e p1q 1 , to vanish independently: It is interesting to note that one can formally plug the Minkowski space with our choice (21) of a g´1η into the naive equation of motion (5), , and deduce the same condition (31) for β-s. Now we look at the second order Lagrangian for perturbations: We see that it contains the root e p1q 1 of the cubic equation (29), and therefore is not analytic in H. If we enumerate the eigenvalues of H such that e p1q 1 " 1 2 p´λ 1`λ2`λ3 q, we get: The three branches of solution correspond to three possible choices of two out of three eigenvalues for putting them into δ p2q L int .

Square roots in the physical dimension
In four dimensions the system of equations (12) would be slightly more complicated: 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: The remaining equations (32) -(34) give at the first order: The system is degenerate, and we can only find that leaving one variable undetermined.
At the second order we get a similar degeneracy, and find The degeneracy gets broken at the fourth order level producing a generic (in terms of independent parameters) quartic equation for e p1q 1 . Let us now have a look at the interaction Lagrangian: he first order contribution to which reads Again, the Minkowski space is a solution when δ p1q L int " 0. And we require the two parts, with h 1 and e p1q 1 , to vanish independently: As in 3D, one can formally plug the Minkowski space with (36) into the naive equations of and deduce the same condition (37) for β-s. Now we look at the second order Lagrangian for perturbations: If we enumerate the eigenvalues of H such that e p1q 1 " 1 2 p´λ 1`λ2`λ3`λ4 q, then we get: 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 ě 2. Then the structure of quadratic 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 X 9 I can be treated as a limiting case of a cosmological matrix diag pN ptq, aptq, . . . , aptqq 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).

4D: the second case
Another option available in four dimensions can be represented by the root a η´1gˇˇg 4 " 1; they correspond to vanishing determinant in equation (19). And the last polynomial (16) is The first order equations (13) -(15) are particularly simple and degenerate: At this point one might start worrying about the self-consistency. Wouldn't the perturbation theory break down at the second order when e p1q 1 " 0? Indeed, given some eigenvalues λ i , there should be six values of e p1q 1 of the form 1 2 pλ 1`λ2´λ3´λ4 q. Can a third order symmetric polynomial h 3´1 2 h 1 h 2`1 8 h 3 1 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 which gives a finite value for e p2q 1`e p2q 3 " Then, at the fourth order, the system of equations (13) 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: and, given also (40), those values of e p1q 1 are indeed the solutions of (42). Let us now have a look at the interaction Lagrangian: ith the first order expression being We look for Minkowski space to be a solution: δ p1q L int " 0, and demand Probably, at this point the reader would not be much surprised if we say that the naive equation of motion which entails the same condition (44) for β-s. The second order Lagrangian for perturbations is hich, given the identifications (40) and (43) above, is equivalent to 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 ď 2 and j ě 3 or vice versa. Then the structure of quadratic Lagrangian has two Fierz-Pauli- 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 g " diagt1, a 2 ptq, a 2 ptq, a 2 ptqu, f " diagtW 2 ptq, Y 2 ptq, Y 2 ptq, Y 2 ptqu with X " g´1f " tW 2 , Y 2 {a 2 , Y 2 {a 2 , Y 2 {a 2 u, the classical equations of motion are well defined only for square roots of a g´1f "˘t˘W, Y {a, Y {a, Y {au type since they do not make equal eigenvalues different. But we must be careful if choosing another square root: a g´1f " tW, Y {a,´Y {a,´Y {au 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

Discussion and Conclusions
It is amazingly difficult to promote the simple action of Fierz and Pauli to a full-fledged nonlinear theory of massive graviton. In the metric formulation one needs to compute the square root function of the matrix g´1f . 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.
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ˆn matrix X can be brought by a similarity transformation to the Jordan normal form which is a block-diagonal form SXS´1 " J k 1 pλ 1 q À J k 2 pλ 2 q À . . . with k i standing for the dimensionality of the i-th block, and the latter is given by with the complex number λ 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 " λ 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 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ˆn matrix X over complex numbers, one can always find an eigenvector 4 , and choosing it as the n-th element of the basis we make X n,k " 0 for k ă n. After that we do the same for the upper left pn´1qˆpn´1q corner of the new matrix, and so on. By induction, we bring X to an upper triangular form. Obviously, the numbers on the diagonal are eigenvalues, i.e. the roots of the characteristic polynomial of 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, λ 1 . Then, we have hereX is non-degenerate and X 0 is nilpotent. If we make a similarity transformation by matrixˆI C O I˙t hen the matrix X´λ 1 I, and also the matrix X, becomes block-diagonal ifXC´CX 0 " M which can be solved as C " ř iě0X´i´1 M X i 0 which is a well-defined finite sum sinceX 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 λ 1 with all eigenvalues equal to a given number λ 1 . Obviously, it suffices to study the nilpotent matrix X 0 " X λ 1´λ 1 I of a given size kˆ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 ď k such that X k 1 0 " 0. Then there is a vector e 0 for which X k 1´1 0 e 0 ‰ 0. This vector together with all vectors X i 0 e 0 generates a vector subspace V 0 in which the matrix X 0 acts as a Jordan block with λ " 0. If k 1 " k, we are done. Otherwise, let's take a suitable 5 vector e 1 from a linear complement to V 0 . We want to build a new cyclic subspace by vectors X i 0 e 1 . This strategy might fail if for some i we get 0 ‰ X i 0 e 1 P V 0 . Then it means that the vector X i 0 e 1 is a linear combination of vectors X j 0 e 0 with j ě i (for otherwise we have X k 1 0 e 1 ‰ 0). Therefore, there is someẽ P V 0 such that X i 0ẽ " X i 0 e 1 , and one can substitute e 1 by e 1 1 " e 1´ẽ 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 ř k i " k. The spectral theorem is proven. Therefore, pJ 3 p0qq 2 is equivalent to J 2 p0q À J 1 p0q under renumbering of the basis' elements. The Jordan block splits into two. Obviously, the splitting is related to the fact that f 1 pλq " 0. Note that it never occurs if f 1 pλq ‰ 0.
Given also that pJ 2 p0qq 2 " O, we see that a 3ˆ3 matrix A such that A 2 " J 3 p0q 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 A¨A " X, we will need to know the full set of solutions for commutation equation of the form AX " XA.
Let us work in the basis in which the given matrix X acquires the canonical Jordan form X " À i J k i pλ i q. This is a block diagonal form, and so let us consider the matrix A in a block form composed of k iˆkj matrices A i,j . Then the equation AX " XA takes the form A i,j J k j pλ j q " J k i pλ i qA i,j (no summation). If λ i ‰ λ j we get pλ i´λj qA i,j " J k i p0qA i,j´Ai,j J k j p0q. The matrices Jp0q in the right hand side are nilpotent. We can multiply this equation by λ i´λ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 λ i " λ j , then we have a simple equation A i,j J k j pλq " J k i pλqA i,j . The Jordan blocks have very simple structure, and one can easily check that if k i " k j (for example, i " jdiagonal blocks) then A i,j is an upper triangular matrix with components pA i,j q a,a`b " C b given by arbitrary numbers C b for b " 0, 1, . . . k´1. If k i ‰ k j then the solution is either O˙w hereÃ i,j is an upper triangular matrix of the same form as presented above and the size of minpk i , k j qˆminpk i , k j q.
If the matrix X was diagonalisable, then all the blocks are 1ˆ1-dimensional, and we simply have that the matrix element A ij is arbitrary if λ i " λ j , and A ij " 0 if λ i ‰ λ j . It is very easy to understand: in those subspaces where the matrix 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.
always T " S, while the general solution is that T S´1 commutes with À i J k i pλ i q. In other words, T is a product of S and an arbitrary matrix which commutes with À i J k i pλ i q. Of course, it gives a new square root only if T S´1 does not commute with À i a J k i pλ i q.
Now, we can see where appears a continuous freedom in square roots. If there are two Jordan blocks in X with equal eigenvalues, λ i " λ j , then a similarity transformation matrix S which leaves 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 ? λ i ‰ a λ j , then such transformation does change the square root ? X. This is how we get the continuous freedom corresponding to the manifold of similarity transformations which change ? X leaving X intact.

A.5 Cayley-Hamilton theorem
Let P pλq be the characteristic polynomial of X. The theorem states that P pXq " O. It trivially follows from our considerations above. Indeed, if there is a Jordan block J k i pλ i q in the spectral decomposition of X, then it means that λ i is a root of P of order not less than k i . Therefore, P pJ k i pλ i qq " O for every block in the spectral representation of X.