Exact ghost-free bigravitational waves

We study the propagation of exact gravitational waves in the ghost-free bimetric theory. Our focus is on type N spacetimes compatible with the cosmological constants provided by the bigravity interaction potential, and particularly in the single class known by allowing at least a Killing symmetry: the AdS waves. They have the advantage of being represented by a generalized Kerr-Schild transformation from AdS spacetime. This means a notorious simplification in bigravity by allowing to straightforwardly compute any power of its interaction square root matrix, opening the door to explore physically meaningful exact configurations. For these exact gravitational waves the complex dynamical structure of bigravity decomposes into elementary exact massless or massive excitations propagating on AdS. We use a complexified formulation of the Euler-Darboux equations to provide for the first time the general solutions to the massive version of the Siklos equation which rules the resulting AdS-waves dynamics, using an integral representation originally due to Poisson. Inspired in this progress we tackle the subtle problem on how matter couples to bigravity and concretely if this occurs through a composite metric, which is hard to handle in a general setting. Surprisingly, the Kerr-Schild ansatz brings again a huge simplification in how the related energy-momentum tensors are calculated. This allows us to explicitly characterize AdS waves supported in one case by a massless free scalar field and by a wavefront-homogeneous Maxwell field in another. Considering the most general allowed Maxwell source instead is a highly nontrivial task, that we accomplish by exploiting again the complexified Euler-Darboux description and taking advantage of the classical Riemann method. In fact, this allow us to find at the end the most general configurations for any matter source.

VII. Matter coupling for AdS waves 9 A. Effective coupling to scalar fields 10

I. INTRODUCTION
Gravitational theories are characterized by the existence of wavy configurations, as in any other relativistic field theory; they are responsible for propagating the interaction and transferring energy across spacetime, which eventually can be measured independently of where they were generated. This was recently achieved by the LIGO and Virgo Collaborations with the direct detection of gravitational waves for the first time [1][2][3][4][5][6], which provided the missing crucial element in Einstein's legacy to understand the gravitational interaction as geometry and was awarded the 2017 Nobel Prize in Physics [7]. However, other open questions remain to be answered, such as the dark matter and dark energy problems, which might be tackled by modifying General Relativity. Therefore, upcoming experimental data require theoretical results capable of illuminating the way to proceed. In this sense, the full understanding of gravitational-wave solutions is promising, since they contain fingerprints to track down features of different theories, especially beyond the perturbative treatment.
In four dimensions, one of these modified theories is the massive gravity of de Rham, Gabadadze, and Toley (dRGT) [8]. It has received a lot of attention for being a consistent way to include nonlinear massive gravitons in a manner that is free of the so-called Boulware-Deser ghosts [9], which permeated many previous attempts. These efforts to include massive gravitons not only followed the theoretical insight of Fierz and Pauli [10]through a long endeavor to consistently incorporate selfinteractions for a massive spin-2 particle-but they also rely on the belief that such massive gravitons can help to explain some currently incomprehensible characteristics of the observed Universe.
The dRGT theory requires a reference metric in addition to the dynamical massive one in order to formally restore diffeomorphism invariance, which entails a precise nonderivative self-interaction potential built with the square root matrix of the product of both metrics. As a consistent extension, Hassan and Rosen [11] showed that the reference metric can be promoted to an independent dynamical one without spoiling the nice features of the original theory (i.e., it is generally covariant and ghostfree), which is the why the resulting theory is called bigravity. This is the framework we will work with, so will devote the following section to reviewing its main ingredients.
The first gravity theory propagating a ghost-free massive degree of freedom (d.o.f.) was Topologically Massive Gravity in three dimensions [12,13]. It relies in a parity-violating Chern-Simons description, which was superseded by later proposals as New Massive Gravity [14] or Zwei-Dreibein Gravity [15]. The latter are theories which in the perturbative regime recover the Fierz-Pauli three-dimensional limit; in this sense, they are analogous to dRGT and bigravity theories in lower dimensions. All of these three-dimensional theories support exact gravitational-wave configurations that exhibit the particular dynamics of each theory beyond the perturbative level, and consist of the pure gauge modes of standard (2 + 1) gravity plus nontrivial contributions reflecting their respective massive excitations [16][17][18][19][20] (see also the recent review [21]). It is expected to have similar revealing behavior in bigravity, i.e., that exact gravitational waves decompose the complex dynamical structure of the theory into elementary exact massless and massive excitations.
Exact gravitational waves over flat spacetime rigged by bigravity were previously investigated in Ref. [22] under the restriction that both metrics have the same profile. As we will see in the next section, the interaction terms inherently include cosmological constants for both metrics; hence, it is natural to search for wave solutions that can be interpreted as being propagated over (anti-)de Sitter [(A)dS] backgrounds. One such configuration is what is known in the literature as AdS waves, which were exhaustively studied in General Relativity by Siklos [23]. They can be represented by a generalized Kerr-Schild transformation from AdS spacetime. Fortunately, it has been proved by some of the authors that dealing with such transformations provides a notorious simplification in bigravity, allowing to straightforwardly compute the interaction square root matrix and any of its powers, independently of the seed metric [24]. Section III reviews the appearance of exact gravitational waves in General Relativity and ends by motivating the AdS-wave ansatz, including its residual symmetries and the reductions they induce in the solutions. Later, Sec. IV deals with the study of AdS waves in bigravity and how to decouple the resulting partial differential equations system for both independent wave profiles. It turns out that one of the decoupled profiles obey a massless Klein-Gordon equation on AdS, which is just the so-called Siklos equation [23]. The other satisfies a massive Klein-Gordon equation on AdS, with an effective mass proportional to the Fierz-Pauli one, defining a massive deformation of the Siklos equation. As a warm up, we start in Sec. V by analyzing sum-separable solutions of these equations to get a clear decomposition in the involved physical modes and gain intuition about the resulting configurations. Thereafter, we address the general setting in Sec. VI by regarding the Siklos operator as a complexified version of Euler-Darboux operators [25,26], which allows us to find the most general AdS waves rigged by the dynamics of bigravity. Another controversial issue of bigravity is the way to couple matter, and Sec. VII is devoted to this issue with a discussion on the effective metric and explicit calculations for both scalar and Maxwell fields. Their general solutions are also characterized by again exploiting the Euler-Darboux description and making use of the Riemann method, which in fact allows us to find the general solutions of the involved exact excitations for any matter source. Final remarks and perspectives of the present work appear in Sec. VIII. The case of pp-waves propagating in flat spacetime is additionally analyzed by completeness in Appendix A. Other appendices are devoted to detailed derivations which are essential to obtain the main results of the paper.

II. THE BIGRAVITY THEORY
Bigravity as formulated by Hassan and Rosen [11] is a four-dimensional ghost-free theory describing two dy-namical metric fields g µν and f µν interacting via a nonderivative potential, originally proposed by dRGT to consistently describe a single massive spin-2 field [8]. One of the metric fields is massive and the other is massless; hence, the theory consequently propagates a total of 5+2 d.o.f. The interaction is encoded through terms computed from a matrix defined by the following quadratic relation The action defining bigravity is where m is the graviton mass, R[g] and R[f ] are the scalar curvatures associated to the metrics g µν and f µν , respectively, κ g and κ f are their corresponding gravitational constants, and κ is in general a function of these constants. The interaction is specified by the potential where b k are the coupling constants of the theory and each term of the potential is defined as U 0 (γ) = where [γ n ] stands for the trace of the nth power of the square matrix γ defined in Eq. (1). Here, the first contribution plays the role of a cosmological constant for the metric g µν . Additionally, the last one can be compactly rewritten in the action (2) as a constant times √ −f , and thus it plays the role of a cosmological constant for the other metric f µν .
The action principle yields a set of two coupled Einstein field equations where G µ ν and G µ ν are the Einstein tensors for g µν and f µν , respectively. The interaction tensors are defined by varying the potential with respect to both metrics where In principle, all of the coupling constants b k are left free. But, as we already denote the graviton mass by m, in the consistent linear limit enjoyed by the theory the latter constant must correspond to the Fierz-Pauli mass in flat space. This imposes a constraint between the coupling constants [24] which we shall use in the present work.

III. EXACT GRAVITATIONAL WAVES: THE CASE OF ADS WAVES
Gravitational waves are most commonly understood as small perturbations to a background spacetime in the form g µν = g (0) µν + h µν , where the quadratic and higher contributions of h µν are neglected, which linearizes the intrinsically nonlinear Einstein field equations and forces the perturbations to satisfy the standard wave equation. This is the kind of gravitational waves recently detected by the LIGO and Virgo observatories [1][2][3][4][5][6]. Remarkably, in spite of its inherent nonlinearity, General Relativity also allows the existence of exact gravitational waves. These are solutions of the Einstein equations for which the linearization leading to the wave equation does not rely on any approximation, but rather emerges from a different mechanism usually involving the existence of a principal null direction [27]. When such a vector field associated to the Weyl tensor presents multiplicity, this defines the so-called algebraically special spacetimes. In the case where this multiplicity is maximal (fourfold), the spacetimes are classified as type N and characterize the exact gravitational waves. If these null vector fields are additionally geodesic (optical rays), 1 they are classified in terms of the irreducible contributions to the projection of its covariant derivative on the two-dimensional spatial sections orthogonal to them. When the antisymmetric part of this projection vanishes (nontwisting rays) and the traceless contribution of the symmetrical part is also zero (shear-free rays), there are only two possibilities: the trace also vanishes (nonexpanding rays), or it is nontrivial (expanding rays). In vacuum, both cases are described 1 It is not necessary to assume this and other properties of the null congruences under certain conditions summarized in the celebrated Goldberg-Sachs theorem [28] and its generalizations [27,29]. by the Kundt [30] and Robinson-Trautman [31] classes of exact gravitational waves, respectively. The nonexpanding Kundt class contains a subcase where the multiple principal null direction is additionally a Killing field and therefore a covariantly constant vector, which describes plane-fronted gravitational waves with parallel rays, or in short pp-waves. This was one of the early exact examples of gravitational waves [32] and probably the most studied. Their dynamics under bigravity is explored in Appendix A.
As emphasized in the previous section, the theory we focus on in this paper naturally presents a pair of cosmological constants. Hence, it is more appropriate to study the kind of exact gravitational waves that can be propagated under these circumstances. In General Relativity, the generalization of the Kundt and Robinson-Trautman waves in the presence of a cosmological constant was realized by the CINVESTAV group in the early 1980's [33][34][35], (see also Ref. [36] and the review [37]). Here again, the generalized Kundt class has a subcase where the nonexpanding ray becomes a Killing vector and which was exhaustively studied by Siklos [23]. The symmetrical Siklos spacetimes only exist for a negative cosmological constant and are defined by the metric where the null Killing field is ∂ v . In General Relativity the gravitational profile F satisfies the wave equation on AdS space written in Poincaré coordinates where u and v play the role of retarded and advanced times, respectively. This permits the interpretation that these solutions are exact gravitational waves propagating on an AdS background [38] or, in short AdS waves [18]. The linearization that gives rise to the wave equation arises in this case because the metric can be written as a generalized Kerr-Schild transformation where the vector field is proportional to the Killing vector ∂ v and retains its null and geodesic properties on AdS. Satisfying such properties on any seed metric assures that the mixed components of the transformed Ricci tensor depend linearly on the profile F and their derivatives [27]. This explains the emergence of the wave operator. The line element (8) where f = f(u), P = P (u), and B 0 = B 0 (u) are arbitrary functions of the retarded time, a dot denotes a derivative with respect to u and the coefficients of the quadratic and linear terms in the wavefront coordinates at the profile transformation are determined from the above functions by These transformations determine the residual symmetries of the AdS waves (8), and they were discussed originally by Siklos in Ref. [23] and were extended to any dimension in Ref. [39]. They can be exploited in the following way (see Ref. [18]): if the solution profile contains a quadratic term in x and y with the same coefficient, a linear term in x, and/or a zero-order term, one can choose the functions in the transformation in order that B 2 , B 1 , and/or B 0 coincide with the coefficients of these terms, respectively. These selections entail differential equations for f and P through Eq. (12b), whose solutions define precise local coordinate transformations that eliminate the involved contributions in the transformed profile. This feature will be of great help, as we will see in Sec. V and later.

IV. AdS WAVES IN BIGRAVITY
In this work we undertake the task of studying the dynamics of AdS waves in bigravity. Hence, we will take both metrics as a generalized Kerr-Schild ansatz of the form (10), additionally allowing the presence of a global conformal factor in the second metric which provides the freedom to use two different AdS radii, whose ratio is precisely C = ℓ f /ℓ g . In the context of bigravity, not only does the Kerr-Schild ansatz provide the well-known linearization leading to the exact wavy behavior of General Relativity described in the previous section, but (as was noticed first in Ref. [24] by some of the authors) the null character of the vector (11) in the generalized transformation also provides a nilpotent contribution to the interaction square where ∆ S is the operator defining the Siklos equation, defined below in Eq. (22a), given by with ∆ L being the standard Laplacian operator. With the help of the profile redefinitions 3 the system (19) is decoupled and becomes The first equation for F corresponds to the Siklos equation, since he was the first to describe the dynamics of AdS waves [23]. Up to a factor, this is just the wave operator (d'Alembertian) evaluated on the AdS metric (9) written in Poincaré coordinates. In other words, the profile F describes an exact massless excitation. The second equation for H is the massive version of the Siklos one since Eq. (22b), up to a factor, is just the massive Klein-Gordon equation evaluated on the AdS spacetime (9) with a mass given in terms of the Fierz-Pauli one as followsm Correspondingly, the profile H characterizes an exact massive excitation. Finally, it is important to know how the exact decoupled excitations are defined modulo diffeomorphisms in order to properly identify their physically relevant contributions. Since we are using the same coordinates to write both metrics (f and g), it is easy to check that the new decoupled profiles change under the residual symmetries (12) of AdS waves in the following waỹ This means that only the massless profile inherits the characteristic indeterminacy of the AdS waves and the massive one remains essentially the same modulo a trivial scaling.
In the following we shall characterize these decoupled exact excitations, first by inspecting their sum-separable sector in order to identify the principal modes ruling the dynamics, and later by unveiling their full space of solutions using Euler-Darboux operators. From the obtained behaviors the original AdS-wave profiles can be reconstructed by inverting the redefinitions (21) as It is very illustrative to start studying the wavefrontcoordinates sum-separable solutions to the decoupled exact excitations (22). This procedure provides a quick integration and, more importantly, reveals the d.o.f. propagated by the theory remaining on these configurations which helps to gain intuition on the exhibited physical modes as well as the special points of the parameter space of the theory.

A. Massless profilesm = 0
The closest case to General Relativity corresponds to the vanishing of the effective mass (22c),m = 0. A particularly interesting possibility is having zero-mass modes without requiring the flat-space graviton mass m to vanish, but rather imposing P 0 = 0. This corresponds to a constraint on the coupling constants of the theory that reduces the coupled differential system (19) to the pair of decoupled Siklos equations As explained before, by looking for a clear decomposition in the prevailing modes of the theory we will search for solutions that are sum separable with respect to the wavefront coordinates, i.e., F 1 (u, y, x) = X 1 (u, x) + Y 1 (u, y) and F 2 (u, y, x) = X 2 (u, x) + Y 2 (u, y). As a consequence of this process, each Siklos equation is separated into ordinary Euler equations for the wavefront coordinates. The linearly independent solutions are power laws and together lead to the more general separable solutions in the massless case where the f 's and h's are arbitrary functions of the retarded time u. We may be tempted to use the residual symmetries (12) to remove all the terms that have an unphysical meaning in standard gravity [18]; interestingly, this can be achieved for only one of the metrics, while in general the other profile will keep all of terms We note that F 1 corresponds to the physical General Relativity mode, which in this context describes the wellknown Kaigorodov spacetime [42]. Instead, F 2 includes additional contributions that cannot be dropped and are fingerprints of the massive d.o.f. of the present theory.

B. Massive profilesm = 0
Let us return to the most general case in which the effective mass is not trivial and the AdS waves are described by the decoupled exact excitations (21) obeying the system (22). Once again, in order to exhibit the decomposition in the prevailing modes we will search for decoupled solutions that are sum separable with respect to the wavefront coordinates, i.e., F (u, y, x) = X 1 (u, x) + Y 1 (u, y) and H (u, y, x) = X 2 (u, x) + Y 2 (u, y). The separation in terms of ordinary Euler equations for the wavefront coordinates now results in the solution where are the roots of the characteristic polynomial determining the linearly independent power-law solutions of the ordinary Euler equation for the massive modes (since their separable x dependence must be trivial), and corresponds to the well-known Breitenlohner-Freedman bound, which is the lowest value the square of the mass of a stable scalar field can take on an AdS background [43]. Since the profile F remains massless we already exploited the residual symmetries (23) to get rid of the unphysical terms (which are like those of the previous subsection) and reduce its solution to Eq. (28a). Hence, the profiles of both metrics that are free from irrelevant contributions are constructed using the inversions (24) and are given by The saturation of the Breitenlohner-Freedman bound, m 2 = m 2 BF , leads us to a logarithmic profile due to the multiplicity in the powers of the solutions to the involved ordinary Euler operator, ρ + = ρ − = 3/2, where we have used the residual symmetry (23) to get rid of the unphysical terms. Again, the profile functions can be explicitly written through Eqs. (24) as As we saw before, one advantage of looking for sumseparable solutions is that it allows to implement the whole residual symmetry [Eq. (12) or (23)] to get rid of the nonphysical terms. This clearly exhibits a mode coming from General Relativity as well as two additional physical massive modes that are characteristic of these kind of theories; these are the prevailing modes respecting the symmetries of the system. For completeness, the same approach is applied to the pp-wave problem in Appendix A, providing new results which extend those already reported for massive gravity [22].

A. General exact massless excitations
In order to find the solution to the system (22) in a general setting, it is quite convenient to incorporate the wavefront coordinates in a complex variable z = x + iy together with its complex conjugatez, in terms of which the Siklos operator (20) can be cast into the form The right-hand side is just a complexified version of the Euler-Darboux operator E α,β , which we review in Appendix C (see also Refs. [25,26]), in the case where the parameters of the operator (C1) take the values α = −1 = β. The complex variables imply that the Euler-Darboux operator E α,β is no longer hyperbolic, but rather elliptic. In other words, the Siklos equation (22a) is just a complexified subclass of the Euler-Darboux differential equations [23]. 4 Consequently, the general solution of the Siklos equation can be straightforwardly read from the related expression (C7) when the parameters of the Euler-Darboux operator are negative integers and using the reality condition on the solution. Concretely, for the Siklos case, α = −1 = β we end with for an arbitrary complex function ω(u, z), holomorphic on the complex wavefront coordinate z. This is just the general solution first reported by Siklos in Ref. [23]. The connection between the Siklos equation and the Euler-Darboux operators was pointed out by Siklos himself [23]. But he did not exploit this fact since he arrived at the solution in a different and clever way by using the following third-order identity, satisfied for any function f , that allows another representation of the Siklos operator in terms of the Laplacian This allows to use the well-known fact that the real part of an arbitrary holomorphic function is harmonic and the general solution to the Laplace equation can always be represented in this way. In this sense, it is not strictly necessary to invoke the alternative derivation we present first. However, as we shall see in the following subsection, the former view is essential to study the general solutions of the massive case for which no previous results are known. Before that, we need to reanalyze the manifestation of residual symmetries since only the massless modes are practically sensitive to them, see Eq. (23). Such symmetries are encoded in the holomorphic function ω(u, z), implying that it must necessarily change under the transformation (12) according tõ Hence, any quadratic, linear, and constant holomorphic dependences on the wavefront coordinates having real coefficients depending on retarded time u can be eliminated from ω(u, z) with the help of residual symmetries.

B. General exact massive excitations
The equation defining the exact massive modes (22b) can be written in complex wavefront variables as (36) where the effective mass is given by Eq. (22c). This massive generalization is no longer described by an Euler-Darboux equation, but fortunately (as we review at the end of Appendix C) it is connected to the extension (C14) of Euler-Darboux operators containing the original massless version and whose behavior can be determined again in terms of the standard Euler-Darboux description. This is achieved by means of the redefinition where is any of the roots defined in (28c), which remarkably allows us to rewrite the massive equation as As a consequence, the nonderivative massive contribution is dropped at the cost of changing the kinetic parameters to α = β = ρ − 1, which of course for a generic value of mass gives nothing more than the general form of an Euler-Poisson-Darboux equation. Notice the special selection of the constant coefficient in front of the right-hand side of redefinition (37a); in contrast to the treatment in Appendix C where all the variables in Eq. (C15) are real, here we need to guarantee that the final result is real. It is necessary to remark that (as is proven in Appendix C) we can equivalently represent the massive configuration in terms of one root or the other, which is why we refer to them generically as ρ in Eq. (37b). In Appendix C we address how to find the general solution to these equations in terms of a superposition constructed by means of an integral representation. However, before addressing the general problem with the methods described in Appendix C, we find it illustrative to first study a particular class for which the profile can be written without the use of integrals; this case involves a discretization of the mass above the Breitenlohner-Freedman bound (28d).

Special case with discrete mass
We start by assuming that the parameters of the Euler-Poisson-Darboux equation (38) take integer values, that is, ρ − 1 = n with n ∈ Z; inserting this into the definition of the exponent (37b) we obtain the following restriction for the effective masŝ i.e., the related configurations involve only discrete values of mass with a gap above the Breitenlohner-Freedman bound (28d) for the first value, making all of them physically acceptable. Additionally, each of these discrete values presents a double degeneration allowing two solutions to share the same mass value: the first characterized by a positive integer and the second by a nonpositive one. These solutions must satisfy Eq. (38), which becomes where ς = ς(u, z) is a complex function that depends arbitrarily on its arguments, but it is holomorphic in the complex wavefront coordinate. Notice that we have made appropriate choices for the coefficients in order to have a manifestly real profile. For more general values of the mass we need a different approach, explained below.

BF
Let us consider now a generic value for the mass in (37b). This requires to integrate the Euler-Poisson-Darboux equation (38) for a generic value of its parameter. The involved solution can be written from the general solutions (C11) to the Euler-Darboux equations we review in Appendix C. Such solution is built as a superposition of two linearly independent particular solutions which are valid when their parameters take generic values such that α + β = 1. For massive configurations α+β = 2(ρ−1) and the integral representation (C11) applies for ρ = 3/2, i.e., the Breitenlohner-Freedman bound (28d) cannot be saturated in Eq. (37b). Therefore, any exact massive mode above this bound,m 2 > m 2 BF , is given by where ϕ and ψ are two arbitrary real functions and ρ is given by Eq. (37b). We remark that consistently the above expression is real. The only effect of taking the complex conjugate H is that the complex argument of the arbitrary functions changes by i.e. the interpolation parameter is reversed from t to 1−t. Using 1 − t as the new integration parameter and the fact that the rest of each integral is invariant under this reversing, the integral representation (42) remains intact and it is concluded that As the last word, a corroboration of the fact that the solution (42) is not the most general one when the Breitenlohner-Freedman bound ρ = 3/2 (m 2 = m 2 BF ) is approached resides in the fact that in such a limit we end up with a solution that possesses a single arbitrary function,φ = ϕ + ψ, instead of two as it should be for a second-order equation. In the following we explain how to appropriately saturate this celebrated bound.

Saturating the Breitenlohner-Freedman bound
For α + β = 2(ρ − 1) = 1 the two solutions from which the superposition (42) is built is no longer linearly independent, see Appendix C. This occurs for ρ = 3/2 or from (37b) when the Breitenlohner-Freedman bound is saturatedm However, intriguingly, the saturated solution can be obtained from the generic one through a nontrivial limit exhibited in general for Euler-Darboux equations in Appendix C. We start with the observation that the generic solution (42) can be rewritten as where the arbitrary functions are properly redefined (see Appendix C). This expression is more appropriate for exploring the saturation of the bound without losing generality. In fact, the most general solution for the Breitenlohner-Freedman mass (43) is just obtained by taking the limitm 2 → m 2 BF in the previous expression This kind of logarithmic behavior emerges when massive configurations approach the Breitenlohner-Freedman bound and is a characteristic of many massive theories [16][17][18][19][20].

VII. MATTER COUPLING FOR ADS WAVES
A widely discussed topic in bigravity is how matter should be coupled to the gravities. Since there is no experimental feedback, there are many possibilities that seem (at least theoretically) consistent. One proposal to democratically couple matter without reintroducing the Boulware-Deser ghost is through the construction of an effective metric [44] A remarkable feature of such a metric is that it is symmetric under the simultaneous exchange of the metrics g ↔ f and couplings α ↔ β for the usual case of interest where the vierbein and metric formalism coincide. Hence, the resulting matter coupling will be symmetric with respect to both metrics as it is, in fact, the vacuum bigravity (2) itself. The full theory is then described by the action where the first contribution stands for the bigravity action (2) and L M is the matter Lagrangian built with the effective metric (46).
The Einstein field equations now must include the contribution of the matter sources where the energy-momentum tensors are defined by In the second equalities we apply the chain rule to rewrite these tensors in terms of the standard energy-momentum tensor with respect to the effective metric, T E µν , which can be calculated as usual. However, the variation of the effective metric requires the knowledge of the variation of the square-root γ matrix, which has a cumbersome structure [45]. Some efforts to avoid such variation were made in Ref. [46] by contracting the Einstein equations (48) with the inverse of an appropriate Jacobian. Notably, this is another difficulty that can be circumvented using the linearizing properties of generalized Kerr-Schild transformations, which turn the computation of such variation into a very easy task. Due to the form of the γ matrix (14) for generalized Kerr-Schild transformations as (13), the effective metric (46) is reduced to the form whose inverse is This gives straightforwardly and, consequently, the energy-momentum tensors contributing to each set of Einstein equations only differ from the canonical one (calculated from the effective metric) by constant factors Another important consequence of the Kerr-Schild ansatz for the present context is that, after fixing the effective cosmological constants as in the vacuum (18), the left-hand sides of both Einstein equations (48) only have contributions along the null ray k µ , see Eqs. (16) and (17). This forces to any matter supporting the AdS waves to behave as pure radiation (a null pressureless fluid). The consequences of the resulting pure radiation constraints have been explored in standard gravity for scalar fields in Refs. [18,47].
Regarding the equations of motion for the matter fields, since they come from matter variation of action (47), they necessarily have the standard form but written in terms of the effective metric. We shall now consider two concrete examples of matter fields coupled to AdS waves in order to test the previous considerations. We start with a massless free scalar field, and follow with the study of the Maxwell field.

A. Effective coupling to scalar fields
Let us first consider the effective coupling to a massless free scalar field-the simplest matter one can think of. The Lagrangian only consists in the kinetic term, but constructed with the effective metric This gives the standard energy-momentum tensor together with the wave equation associated to the effective metric, which rules the scalar field dynamics The scalar field is easily integrated from the emerging pure radiation constraints, i.e., the vanishing of all the components of the energy-momentum tensor except the one along the null ray, as is imposed by the Einstein equations [18,47]. Considering the combination the result is that the scalar field is an arbitrary function of the retarded time which automatically satisfies the wave equation (56). This result is exactly the same even if one attempts to promote the scalar field to be self-interacting by adding a potential to the Lagrangian (54). The outcome is that no potential is compatible with supporting an AdS wave unless a non-minimal coupling to the involved metric is also incorporated [18,47]. Hence, considering only minimal coupling to the effective metric, the more general situation is just that of a massless free scalar field. So far, we are left only with a contribution along the null ray in both Einstein equations which gives rise to two inhomogeneous differential equations for the profiles F 1 and F 2 , that can be decoupled as in the vacuum (59b) Because we are dealing with an inhomogeneous linear system, its most general solution is built by superposing the general solution to the homogeneous (vacuum) version of the equations with any particular solution of the inhomogeneous (with sources) ones The homogeneous version corresponds to the vacuum system (22), whose general solution was studied in detail in the last section and is given by the Siklos solution (33) for the massless case F h and by the integral representation (42) for the massive one H h . In order to incorporate the inhomogeneous contributions, due to the simple form of the scalar inhomogeneities in Eqs. (59) it is enough to look for particular solutions which are independent of the x coordinate, then, the above equations become inhomogeneous ordinary Euler equations for the y coordinate. The solutions are easily found and they are proportional to the power exhibited at the inhomogeneity, except when that power resonates with one of the vacuum power-law modes (28). 5 The result is the following 5 We emphasize this is a genuine resonance phenomenon. Notice we can change to a different coordinate y = ℓ exp(t/ℓ) for which the ordinary Euler equations become linear equations with constant coefficients, where the resonance phenomenon is normally defined. There, the power-law solutions become exponentials and the power exponents become mass-dependent frequencies. The resonance is associated to the precise values of the parameters specifying the system and is independent of the variables chosen to describe it.
where the scalar source inhomogeneity enters in resonance with the vacuum modes when the mass becomes What happens in this case is that one of the vacuum powers (28c) becomes ρ + = 2 and is equated by the scalar source value. We stress that although the resulting scalar resonant mass is negative, it describes physically admissible configurations above the Breitenlohner-Freedman bound With respect to the inhomogeneous contributions (61), nothing special occurs for the rest of the masses, even for the Breitenlohner-Freedman bound. In fact, the logarithmic behavior atm 2 = m 2 sr of (61b) has a resonant origin and is different from the one appearing at the vacuum for the multiplicity that occurs when the Breitenlohner-Freedman bound is saturated.
Superposing these inhomogeneous contributions with those already studied for the vacuum according to Eq. (60), we obtain the most general form in which a scalar field can support bigravity AdS waves. In what follows we study the related, more complex problem for a Maxwell field.

B. Effective coupling to Maxwell fields
The following is another natural scenario in which bigravity could interact with matter. We are interested now in how electromagnetic radiation fields bend the AdS waves. For this purpose, consider the Maxwell Lagrangian constructed with the effective metric, where the electromagnetic strength is given as usual in terms of the vector potential, F µν = 2∂ [µ A ν] . The standard electromagnetic energy-momentum tensor resulting from varying with respect to the effective metric is Maxwell equations in terms of the effective metric are now obtained taking the variation through the vector potential After a nontrivial gauge-fixing procedure described in Appendix B, it is shown that the most general Maxwell potential A µ supporting AdS waves, i.e., that is compatible with the pure radiations constraints, is proportional to the null rays k µ The Maxwell equations (64) reduce to the harmonic equation whose general solution is the real part of a general holomorphic function of the complex wavefront coordinate where, for later convenience, we choose to write the holomorphic function as the derivative of other holomorphic function a(u, z). The only non-vanishing components of the Faraday strength tensor are and the effective energy-momentum tensor (63) acquires the form of a pure radiation field [see (B3)] which contributes to the equations for the AdS-wave profiles in the form of inhomogeneities coming from the Maxwell field. The resulting equations can be decoupled using the same combinations (21), just as in the vacuum case. Thus, we get the inhomogeneous second order system The general solution is again represented by the superposition (60) with the same homogeneous contributions F h and H h given by the vacuum configurations (33) and (42), respectively. The inhomogeneous contributions F i and H i are due now to the more complex Maxwell source, and understanding their behavior in the more general setting requires a different treatment.
An interesting point to be noticed is that for both studied sources, if the gravitational constants are tuned by the condition the massive profile H becomes decoupled from the sources and behaves exactly as the vacuum configurations studied in Sec. VI. In order to fully understand the Maxwell contributions to bigravity AdS waves, we shall proceed by gradually increasing the degree of difficulty. We start by analyzing the massless sector,m 2 = 0, which resembles that of General Relativity with a Maxwell source originally studied by Siklos in [23].

Massless sector
The inhomogeneous Siklos solution for the Maxwell source can be straightforwardly read from the following fourth-order identity relating the Siklos and Laplacian operators for any function f , see Ref. [23], Applying the identity to the massless sector (m 2 = 0) of the system (70) we obtain that the inhomogeneous contributions are given by Since the profile F describes a massless mode, their inhomogeneous contribution (73a) remains the same even for more general values of the mass. Hence, the most general massless configuration supported by electromagnetism is given by the superposition profile (60) together with the Faraday strength following from the vector potential (65), which here take the explicit forms [23] We know that the massless profiles are determined modulo residual symmetries according to the transformation (23). This imposes that pairs of holomorphic functions-determining the vacuum and electromagnetic contributions-that are related under the transformation (12) as all represent the same configuration. Here, C 0 = C 0 (u) and C 1 = C 1 (u) are additional arbitrary complex functions of the retarded time which parametrize the indetermination of the gauge field and which naturally induce a generalization of the residual symmetry already known for the vacuum (35). This transformation was unveiled by Siklos in Ref. [23], without incorporating the diffeomorphic part. As he emphasized, it is difficult to obtain due to the quadratic contribution of the electromagnetic holomorphic dependence in the massless geometric profile (74a). This exhausts the understanding of the massless excitations in presence of Maxwell sources. The situation is different for the profile H which is massive in nature; the inhomogeneous contribution (73b) represents only a particular case and the general solution for generic values of the mass with Maxwell sources is a more difficult problem. There is no obvious generalization of the identity (72) which would eventually allow to also derive a local expression for the massive solution. However, the inhomogeneous solutions of any linear partial differential equation generically allows an integral representation. Due to the connection of the Siklos operator with the Euler-Darboux ones, we exploit this fact in subsection VII B 3 to use a complexified version of the Riemann method for hyperbolic equations, which we briefly review in Appendix D. But first, we attack the particular case of a Maxwell field whose strength is homogeneous in the wavefront coordinates which ends up possessing an extra symmetry; in this way the solution allows a local representation.

Massive sector: wavefront-homogeneous Maxwell source
Now we consider a particular example inspired by the fact that the separable vacuum solutions studied in subsection. V B, modulo the use of residual symmetries, are AdS waves allowing as an additional Killing vector the spatial translations ∂ x . Hence, we assume now that the Maxwell field is compatible with such symmetry and consequently ∂ x F ux = 0 = ∂ x F uy . Taking into account the general form of the Faraday strength (68), these conditions imply where D 0 = D 0 (u), D 1 = D 1 (u) and D 2 = D 2 (u) are arbitrary complex functions of the retarded time. Additionally, D 0 and D 1 can be eliminated by using the residual transformation (75). The resulting strength (74b) is in general homogeneous in all of the wavefront coordinates. Keeping in mind that the full electromagnetic contribution to the inhomogeneity of the equations is now independent of the wavefront coordinate x, it is enough to consider particular solution constructed in the same way, i.e. F i = F i (u, y) and H i = H i (u, y). With this, the system (70) becomes a pair of ordinary Euler equations Consequently, for a Maxwell field homogeneous in the wavefront coordinates the inhomogeneous contributions in the decoupled AdSwave profiles are given by where once again the inhomogeneity due to the electromagnetic source enters in resonance with one of the vacuum power-law modes (28) when the mass takes the value The electromagnetic resonance is produced because the power-law dictated by the electromagnetic source equals the vacuum power ρ + in Eq. (28c), which takes the value 4 at the above mass. The full solution is obtained by superposing these inhomogeneous contributions with the vacuum ones following Eq. (60), which describes all AdS waves of bigravity supported by electromagnetic fields homogeneous in the wavefront coordinates. The most general solutions with no additional restrictions on the Maxwell sources are addressed in the next subsection.

Massive sector: general Maxwell source
As has been pointed throughout the paper, the Siklos operator and its massive generalization can be expressed as Euler-Poisson-Darboux operators using as independent variables the complex wavefront coordinate z = x+ iy and its complex conjugatez, together with the redefinition (37) for the massive profile. This allows to describe the behavior of massive excitations in the presence of general Maxwell sources (70b) by means of an inhomogeneous Euler-Poisson-Darboux equation (79) The particular solutions of all inhomogeneous Euler-Darboux equations (C20) can be obtained using the Riemann method reviewed in Appendix D. This allows to express the particular inhomogeneous solution to massive excitations sourced by general Maxwell fields according to Eqs. (37), (D6c) and (D9) as where 2 F 1 describes the standard hypergeometric function.
We have used the double triangular integral representation (D10) with negative-inclined hypotenuse together with the election (D12) for the solution, but the double representation with positive-inclined hypotenuse (D13) is equally adequate. This is where these double representations become relevant; despite its complexity, they significantly ease the construction of a real elliptic solution from the hyperbolic one. In order to elucidate this far from obvious fact, we emphasize first the following points: the nontrivial dependence of the exact massive excitation in the complex wavefront coordinates emerges by evaluating the solution to the standard hyperbolic inhomogeneous Euler-Darboux equation, as defined by the double triangular integral (D10), in the following way The inhomogeneity in Eq. (79) is a real function of the originally real wavefront coordinates, which means that when it is written in terms of the dummy variables inside the integral (80) obeys that Finally, the Riemann-Green function corresponding to the Euler-Darboux equation (D9) is real for real arguments and additionally symmetrical under the simultaneous interchange of each pairs of variables, i.e.
We are now in a position to prove the reality of the so-lution (80); taking its complex conjugate we obtain where the first and second equalities are self-explanatory, while in the third one we use the interchanging properties (82) and (83) and later we reparametrize the dummy variables by (ξ ′ , η ′ ) → (ξ ′′ = η ′ , η ′′ = ξ ′ ). In terms of the new variables both integrals are interchanged, which justifies the fourth equality. This proof also works following the same steps if one uses the integral representation for a triangle with positive-inclined hypotenuse (D13). Another unavoidable exercise is to reconsider from the present perspective the massless scenariom 2 = 0 (ρ = 0). This implies analyzing the relation between the integral representation for the corresponding particular inhomogeneous solution obtained from the Riemann method and the local one already exhibited in Eq. (73b) that was originally provided by Siklos [23]. In fact, for ρ = 0 the integral expression (80) allows a double integration by parts. We find it more easy to do the calculation using the analog triangular integral representation with positive-inclined hypotenuse (D13), which gives where in the second equality we used the previously proven fact that the second integral is the complex conjugate of the first. Additionally, the integrand can be expressed as the mixed derivative of the following function whose proper evaluation gives the Siklos inhomogeneous local solution as the first contribution of the third equality The remaining contribution defined as H h is just a solution to the homogeneous equation, which finally proves the compatibility between the Riemann and Siklos particular inhomogeneous solutions in the massless scenario. Even though the result is more easy to obtain with the above representation, it is equally valid if the integration is performed over the triangle with a negative slope.

C. Effective coupling to any matter source
The Riemann method works so efficiently and in such a general fashion, that one can go even one step further. Similar to the previously described AdS waves supported by a Maxwell field, one can work in a more general setting where the spacetime ripples of bigravity are caused by any kind of matter. The first step is to consistently solve the related pure radiation constraints, which by the proportionalities (53) are all expressed in terms of the energy-momentum tensor associated to the effective metric. Using Bianchi identities this also entails solving the related field equations that are naturally built on the effective background. If a nontrivial answer emerges from this process, it only remains to evaluate the surviving contribution to the energy-momentum tensors along the retarded time and write the single Einstein equation of each set. Considering as before the decomposition (37), the dynamics of the decoupled exact excitations are rigged now by the pair of inhomogeneous Euler-Poisson-Darboux equations Thereby, in this more generic instance, they would be given by the following expressions This completes the exhaustive scan of the AdS-wave configurations in bigravity.

VIII. CONCLUSIONS
In the preceding work we have tackled the problem of searching for exact gravitational waves in the context of the ghost-free bimetric theory. Due to the inherent presence of cosmological constants in this theory, we focused on the Kundt class that is compatible with them, and in particular on the single example where the involved nonexpanding null ray becomes a Killing symmetry: the AdS waves. These configurations are additionally characterized by admitting a Kerr-Schild representation. From previous results [24], we are aware that this is an exceptional simplifying tool in the Hassan-Rosen theory not only for exactly linearizing the Einstein equations as in General Relativity, but also for the usually involved task of computing the interaction between the metrics. The massive nature of bigravity manifests itself in a very clear way, since the resulting wave equations for the profiles are explicitly coupled by the mass term. Even though they can be decoupled through redefinitions, one of the new profiles is inherently massive. We learned how the physical d.o.f.-compatible with the symmetries of the problem-are carried by the waves by looking for configurations that are sum separable with respect to the coordinates defining the wavefront. In this respect, another consequence of the presence of two metrics comes into play even before turning on the mass. It is well known that spaces belonging to the Kundt class possess residual symmetries [27], meaning that some of the metric potentials can be gauged away with appropriate transformations. In bigravity, we can only perform such transformations on one of the metrics, implying that the d.o.f. that would be wiped out in standard gravity remain present here, and actually propagate physical modes of the theory.
Regarding the AdS-wave solutions, we can distinguish two cases in every setting we studied. The massive profile depends on an effective mass that can be turned off even if the flat-space Fierz-Pauli mass remains finite. This is achieved by introducing a fine-tuning between the couplings, which produces two copies of exact massless profiles. The general solution to this branch in General Relativity was already reported by Siklos [23], whereas here we explored the consequences of its observation in the sense that the equation ruling the behavior of the wave profiles (now known as the Siklos equation) turns out to be a particular complexified case of the Euler-Darboux equations. This proves to be useful to understand the more complex scenarios described below. The more general situation involving no fine-tuning corresponds again to a massless profile, but it also includes a massive one. The massive behavior is described by an extension of the previous Euler-Darboux operators. One of the properties of the extended Euler-Darboux equations consists in the possibility of eliminating nonderivative contributions through a redefinition, returning to the standard Euler-Darboux equation. As consequence, the massive part can be removed at the cost of generalizing the parameters of the kinetic contributions. Once again, this allowed us to provide the most general solution to the problem: a closed local form where the Euler-Darboux parameters are integers, which is valid for a discrete family of mass values, and an integral representation (originally due to Poisson) for general mass values.
Since the wave profiles obey the (massive) wave equation over the AdS spacetime, the well-known criterion of stability on AdS D of respecting the Breitenlohner-Freedman boundm 2 ≥ −(D − 1) 2 /4ℓ 2 , appears repeatedly on the massive sector of the problem. With this in mind, it is interesting to inspect the separable solutions. Form 2 = −9/4ℓ 2 , we have profile solutions in powers of the conformal AdS coordinate y, which include the standard Kaigorodov massless contribution [42] that cubically decays to infinity (located at y = 0 in the AdS background in which the modes are propagated) plus two other massive modes. For −9/4ℓ 2 <m 2 < 0, the massless Kaigorodov contribution becomes subleading at infinity with respect to the decay of both massive contributions. However, form 2 ≥ 0 the leading massive contribution no longer decays at infinity. In contrast, the other massive term decays and becomes subleading even with respect to the massless Kaigorodov one. Curiously enough, for the precise valuem 2 = −9/4ℓ 2 (corresponding to the critical mass of the Breitenlohner-Freedman bound) there is a degeneracy in the massive solutions, leading to logarithmic AdS waves-a situation also observed in 2 + 1 massive gravities at critical points where they are supposed to be dual to logarithmic conformal field theories [16][17][18][19][20]. This situation was replicated in the most general solution for which the integral representation of the exact massive excitation also acquires a logarithmic dependence when the Breitenlohner-Freedman bound is saturated.
The natural step to follow was considering the wave dynamics in the presence of matter; however, this question should be carefully examined in bigravity. The problem of matter couplings has been widely discussed and the proposed alternatives should always be concerned with not awaking again the undesired d.o.f. In this work we tested the approach of introducing an effective metric-a composite of the two original ones-constructed in such a way that the whole theory is symmetric when replacing one original metric by the other up to coupling redefinitions. In general, calculating the two involved energymomentum tensors through this effective metric is not an easy task, but we found that even in this regard the Kerr-Schild ansatz provides another significant simplification. By definition, the effective metric contains a mixing term that is the product of one of the metrics with the coupling square root matrix; this term is the origin of the difficulties when performing metric variations to deduce the energy-momentum tensors. Interestingly enough, under the adoption of Kerr-Schild forms this interacting term is linearized and the effective metric results in a linear superposition of both metrics. This leaves a straightforward calculation in the variational problem and both energy-momentum tensors become proportional to the standard one written in terms of the effective metric.
After unraveling how to calculate the energymomentum tensors, another relevant well-known issue when AdS waves are supported by matter sources is that their Kerr-Schild structure is also replicated at the level of the Einstein equations. In fact, their geometric lefthand side has only a contribution along the ansatz null ray, which forces the involved matter to behave effectively as pure radiation, i.e., a null dust. This entails solving the resulting pure radiation constraints for the given source, and only when the result is nontrivial can the AdS waves be supported by this kind of matter [18,47]. The first source analyzed was a free massless scalar field. This situation is the simplest one, since the outcome of solving the pure radiations constraints is that the most general scalar field supporting an AdS wave is just an arbitrary function of the retarded time. The resulting inhomogeneity added to the Siklos equation is a simple quadratic dependence on the wavefront coordinate y. The corresponding particular inhomogeneous part of the solutions obeys exactly the same dependence since it can be assumed homogeneous in the other wavefront coordinate, which reduces the Siklos equation to an Euler ordinary equation. There is an exception when the mass value is such that the scalar source enters in resonance with behavior characteristic of the vacuum; this situation gives rise to logarithmic modes accompanied by the quadratic resonance power.
The second type of matter source studied-the Maxwell field-is quite richer. The most general selfgravitating vector potential that can be constructed that is compatible with the pure radiation constraints has a single component along the retarded time direction, which must be the real part of an arbitrary holomorphic function of the complex coordinate formed with the wavefront coordinates. This function also contains an undetermined dependence on the retarded time. The resulting inhomogeneities added to the Siklos equation are rather more involved now, since they possess a quadratic dependence on the complex derivatives of this holomorphic function. The most general solution to the inhomogeneous part for the massless cases was already presented by Siklos himself [23]. The tough part comes when dealing with the sourced effectively massive equation. A preliminary simple case is to consider a homogeneous Maxwell strength. This reduced the problem of finding a particular inhomogeneous solution to solving again an ordinary Euler equation, but this time with a quartic inhomogeneity in the conformal coordinate y. The resulting solutions share the same y dependence for all masses, except for an electromagnetic resonant mass value for which the companion coefficient develops the characteristic logarithmic behavior of resonance phenomena. Describing the massive AdS-wave dynamics under the most general Maxwell source is a highly nontrivial task. Here, once again, we exploit the connection with the hyperbolic Euler-Darboux equations, whose inhomogeneous versions can be solved by means of the Riemann method. It consists in providing the general solution to a given initial value problem by solving a related characteristic boundary value problem, giving the so-called Riemann-Green function, which acts as the kernel in an integral representation of the inhomogeneous solution. For the Euler-Darboux equations it can be justified that their Riemann-Green functions are determined in terms of hypergeometric functions [26]. Adapting this approach to our elliptic complexified paradigm allows us to provide the general solution for the exact massive (massless) wave profiles not only when they are supported by a Maxwell field, but also in the case when the source is generalized to any matter consistent with the pure radiation constraints. Whether the associated integral representation can be given a closed local form or not depends on the specific value of the mass. For example, for zero mass it can be shown that the particular Riemann solution consistently reduces to the Siklos one modulo homogeneous solutions. A curious situation (common to all scenarios) is that the matter can be decoupled from the massive sector by means of a fine-tuning in the parameters defining the effective metric.
To complete our work, we additionally studied the situation in which the exact gravitational waves are propagated over flat spacetime. This involves supplementary constraints on the bigravity coupling constants in order to get rid of their naturally defined effective cosmological constants. Within the Kundt class valid under this circumstance, the exceptional cases allowing isometries are the well-known pp-waves. Again, these exact waves were decomposed into a massless excitation and a massive one. We analyzed in detail the configurations that are sum separable with respect to the wavefront coordinates. For the zero-mass case, we obtained two decoupled massless profiles representing linearly polarized plane waves (standard in General Relativity), plus linear and homogeneous terms in the wavefront coordinates which are usually gauged from pp-waves. However, here the residual symmetries allow to remove the extra terms of only one of the profiles; they remain untouched in the other and now describe the propagation of the genuine physical d.o.f. of bigravity. In the properly massive case, the decoupled solutions describe on the one hand the linearly polarized plane wave of General Relativity for the massless mode, and on the other hand a superposition of Yukawa exponential decays and growths in each wavefront direction, corresponding to massive modes on flat spacetime.
Finally, we emphasize that not only are the methods introduced in this paper to deal with exact gravitational waves in the presence of a cosmological constant original in the context of bigravity, but also that there have been no similar studies even in General Relativity. We believe that these methods are useful beyond this setting, since they teach us how to find the general solution to the Klein-Gordon equation with sources on the AdS background, for configurations which are only restricted to be invariant under a single light-cone translation. It is not too risky to conjecture that these techniques could have potential applications, for example, in the AdS/CFT context. As was briefly reviewed in Sec. III, the pp-waves are characterized within all of the exact gravitational waves by allowing a covariantly constant null vector field (parallel rays) and the requirement of propagating plane wave-fronts, which is summarized in the line element [32] where x = (x, y) denotes a Euclidean vector. Consequently, the multiple principal null direction is at the same time a Killing field ∂ v and a gradient Hence, the pp-waves are the simplest nontrivial example of a Kerr-Schild transformation from Minkowski flat spacetime and again the null and geodesic character of k warrants the linearization, allowing the profile F to obey the wave equation in General Relativity. Thus, they are interpreted as exact gravitational waves propagating on the flat background.
The pp-wave metric (A1) is form invariant under the family of transformations [27] u = λ (u + u 0 ) , where u 0 , λ, and the matrix Λ ↔ ∈ SO(2) are constants, P = P (u) and B 0 = B 0 (u) are arbitrary functions of the retarded time, a dot denotes derivative with respect to u, and the coefficients of the linear terms in the wavefront coordinates at the profile transformation are determined from the above functions by Similar to the case of AdS waves, if the solution profile F contains linear terms in the wavefront coordinates x and/or a term which is only function of the retarded time u, such contributions can be eliminated by appropriate diffeomorphisms; see Ref. [17]. The starting point to explore the behavior of pp-waves in bigravity is a pair of Kerr-Schild ansatz, where the second metric is allowed to contain a global conformal factor As was previously argued in Sec. IV, the interaction square-root matrix has the form (14) for any Kerr-Schild ansatz and the interaction terms reduce to the same structure as the Einstein equations (16). The fact that the wavefronts of the pp-waves are planes implies there is no longer a diagonal contribution in their Einstein tensors as those of the tensors (17); consequently, the effective cosmological constants (18) do not appear in this case, which imposes the constraints The only nontrivial contributions are those of the offdiagonal components along the null ray This system becomes decoupled for the same combinations F and H decoupling the AdS waves (21), which yields the following equations defining again exact massless and massive excitations propagating now in flat spacetime The effective mass appearing in the later Klein-Gordon equation is the same as that already defined in Eq. (22c). These decoupled profiles are determined up to the residual symmetries of the pp-waves (A4), since they change according tõ Similar to the case of AdS waves, only the massless profile inherits the indeterminacy and the massive one is almost preserved modulo a trivial scaling. The decoupled exact massless profile satisfies the wave equation which, due to the Killing vector ∂ v , becomes just the harmonic equation (A8a) with general solution where G(u, z) is any function depending arbitrarily on the retarded time, but that is holomorphic in the complex wavefront coordinate z = x + iy. A particularly interesting case is the fully massless onem = 0 (P 0 = 0). Looking back to Eq. (A7), both of the original profiles are harmonic and thus are described as in Eq. (A10). Even with this quite general solution, it is useful to explore the particular cases that are sum separable with respect to the wavefront coordinates since they also bring here a clear decomposition in the prevailing modes of the theory after removing the unphysical contributions by means of the pp-waves residual symmetry (A4). Thus, we will search for the most general solutions of the form F 1 (u, x) = X 1 (u, x) + Y 1 (u, y) and F 2 (u, x) = X 2 (u, x) + Y 2 (u, y). Such solutions are where all functions of the retarded time are arbitrary and we have eliminated the linear and homogeneous terms in the wavefront coordinates appearing in the first profile by using the residual symmetry (A4). The quadratic contributions to both profiles are just the sum-separable part of the well-known plane waves of General Relativity [27], which corresponds to the linearly polarized ones [48].
Applying the same analysis to the general massive casê m = 0, the sum-separable solutions of Eq. (A8) are given by where e ±m x = (e ±mx , e ±my ) denote Euclidean vectors and the residual symmetry (A9) allows us to get rid of the unphysical contributions in the massless profile.
Here, the massless profile again corresponds to the linearly polarized plane-wave contribution of General Relativity, and the massive one describes a superposition of Yukawa exponential decays and growths in each wavefront direction, corresponding to standard massive modes in flat spacetime. From the inversions (24), one gets the solutions for the original profiles Let us note that it is not straightforward to obtain thê m = 0 solution from the last outcome. We end here the revision of bigravity pp-waves.
Appendix B: Most general vector potential supporting AdS waves As emphasized in Sec. VII, the fact that Einstein tensors have the structure (17) for AdS waves (8) fixes the effective cosmological constants coming from the interaction contributions as in Eq. (18), which in turn implies that both Einstein equations (48) establish the vanishing of all of the components of the energy-momentum tensors (53) that are not exclusively in the direction of the null ray k µ . In particular, for a Maxwell field the following constraints must be satisfied These imply that the involved strength components vanish and that the energy-momentum tensor (63) acquires the form of a pure radiation field along the null ray Let us start by analyzing the consequences of the first three pure radiation constraints (B2) with x a = v, which can be integrated as whereÂ a are integration functions which are independent of the null ray parameter v. Consequently, modulo the gauge transformation we end with a gauge wherẽ Implementing now the last pure radiation constraint (B2) in the above gauge lead us toÃ The gauge (B7) is still preserved by residual gauge transformations, from which we choosẽ which allows a further fixing after using Eq. (B9) Hence, satisfying the pure radiation constraints entails the existence of a gauge where the vector potential can be expressed as It remains to explore the repercussions of this gauge fixing on the Maxwell equations (64), which become This can be understood as an inhomogeneous linear partial differential equation for the retarded time component A u , where the inhomogeneity is determine by the spatial component A y . The most general solution is a superposition of the kind (60) for the potential A u = A h u + A i u . The contribution A h u must be harmonic since it represents the general solution to the related homogeneous equation, and can be written as the real part of a general holomorphic function in the complex wavefront coordi- Here, for further applications it is convenient to write the holomorphic function as the derivative of another holomorphic function a(u, z). Regarding the inhomogeneous contribution A i u , since only a particular solution to the inhomogeneous equation is needed, we can assume that it does not depend on the spatial coordinate x, i.e., y). Thus, the inhomogeneous equation (B13) is reduced to and is integrated as The simplest particular solution is obtained by choosing K 1 = 0 = K 2 . Consequently, the most general solution to the Maxwell equation (B13) is general is separable into left and right movers encoded in the arbitrary functions f and g. The general solutions for other integer values of the parameters can also be easily expressed. In order to accomplish this, another important feature to realize is that the different solutions are related not only algebraically [as in Eq. (C2)], but also through differentiation. For example, if we know a solution of E α,β (u) = 0, their derivatives are also solutions, in this case to E α+m−1,β+n−1 (u) = 0. In particular, by evaluating α = 1 = β, we have that (C6) Now, by again evaluating α = 1 = β and redefining the integers as m ′ = n − 1 and n ′ = m − 1, we obtain Hence, the general solutions to the Euler-Darboux equations where the parameters take nonpositive integer values are once more generated from Eq. (C3) by differentiation. As is explained in Sec. VI, the Euler-Poisson-Darboux subcase u(−1, −1) is just the solution of the Siklos equation characterizing the AdS waves of General Relativity in presence of a negative cosmological constant [23].
We now consider the case where the parameters α and β appearing in Eq. (C1) are real numbers. The first step is to look for particular separable solutions of the form By inserting this into the Euler-Darboux equation (C1), we obtain the following solution where a is the separation constant. Moreover, we can use the relation (C2), when α + β = 1, to construct another particular solution of E α,β (u) = 0 as follows The previous particular solutions are the seeds to construct the general solution to the linear Euler-Darboux equation (C1) as a superposition where ϕ and ψ are arbitrary functions and in the second equality a new interpolation parameter has been introduced, a = ξ + (η − ξ)t. This integral representation for the general solution was provided first for α = β by Poisson [51], and then generalized to α = β by Appell [54]. The case α + β = 1 is not covered by the expression (C11); we cannot use the relation (C2) to construct a second independent particular solution. However, interestingly enough, we can use the final form (C11) to derive the general solution of this case as a nontrivial limit. First, we notice that the generic solution (C11) can be rewritten in the following form where, without losing generality, we are considering that the arbitrary functions could be defined from the beginning asφ = ϕ + ψ,ψ = (α + β − 1)ψ.
The general solution to the Euler-Darboux equation when α + β = 1 is precisely the limit of Eq. (C12) when the sum α + β approaches unity; therefore which is just the general solution reported in Ref. [25] without justification. This kind of solutions has turned out to be relevant in General Relativity [specifically u(1/2, 1/2)] in the construction of the so-called Weyl class describing static axisymmetric vacuum spacetimes [55], or in their analogues after a double Wick rotation characterizing colliding plane gravitational waves with collinear polarization [56]; see Refs. [27,48] for more recent reviews on these subjects.
In what follows we emphasize there are extensions of the Euler-Darboux equations [52] which can be recast into the standard form. In the previous case, this is done by just letting Indeed, u must obey the standard Euler-Darboux equation provided that ρ is a solution to the quadratic equation In this regard, it is important to make the following remark. In general, the latter polynomial has two roots, say, ρ ± . Therefore, from Eq. (C15) one may be tempted to write the general solution to the extended Euler-Darboux equation (C14) as a superposition built from the contributions of both roots. However, it is enough to consider a single root since, using relation (C2), it is easy to check that where d = (α+ β − 1) 2 − 4λ is just the discriminant of the quadratic equation (C17). Hence, the general solution to the extension (C14) is the same when written with one root or the other For α = −1 = β, these are precisely the solutions of the massive generalization of the Siklos equation we report in Sec. VI and where the third parameter λ is determined by the mass. More relativistic examples for which the extended operator (C14) becomes relevant were reviewed in Ref. [57], where almost all aspects covered in this appendix and the following one were also reviewed from a similar perspective. Finally, we also need to address the problem of how to solve the inhomogeneous Euler-Darboux equations since it is important in this paper to understand the behavior of AdS waves with matter fields as sources. The inhomogeneous solutions of linear hyperbolic equations such as these can be obtained by using the Riemann method. In fact, Riemann devised his method by studying the Euler-Poisson-Darboux equations [58]. Hence, for completeness, we briefly review this method and how it must be applied to the Euler-Darboux case in the last appendix.

Appendix D: The Riemann Method
Here we follow Ref. [26] but adopting a notation compatible with the previous appendix. First of all, let us consider a linear second-order differential operator with continuously differentiable coefficients. Multiplying the above by a smooth function R and differentiating by parts, the following identity is obtained where the functions determining the divergence on the right-hand side are not unique since we can add to them ∂ η Θ and −∂ ξ Θ respectively. However, the linear operator is unique and defines the adjoint of L. Consequently, the adjoint of L * is L. When additionally L * = L it is said that the operator L is self-adjoint. By integrating the identity (D2a) in a well-defined domain Ω whose boundary is a regular closed curve ∂Ω with anticlockwise orientation, it follows from Green's theorem that One of Riemann's significant contributions consists in an approach to solve the Cauchy problem of an inhomogeneous hyperbolic equation using the above Green integral identity [58]. For a hyperbolic operator there exist characteristic coordinates that allow to get rid of all of the second-order derivatives except the crossed one; as a result, the associated inhomogeneous equation can be rewritten as where the inhomogeneity f is a continuously differentiable functions, as are the coefficients D, E and K. Solving the Cauchy problem means finding the unique solution to this equation with given values of u, ∂ ξ u and ∂ η u on some initial curve. Let P = (ξ 0 , η 0 ) be the point in the future (or past) of the initial-value curve where it is required to know the solution. The first step of the Riemann method is to choose the regular boundary ∂Ω of the identity (D3) in such a way that it contains both the initial-value curve and the future (past) point. In the characteristic plane (ξ, η) the initial-values curve must be a duly-inclined regular arc, and we can define its intersections with the characteristics η = η 0 and ξ = ξ 0 passing through P as the points Q = (ξ 1 , η 0 ) and R = (ξ 0 , η 1 ), respectively. Hence, the regular boundary ∂Ω can be the deformed triangle with vertices P QR. The second step of the method involves picking out the function R in the identity (D3) as a homogeneous solution to the adjoint operator L * (R) = 0, satisfying the additional conditions R(ξ 0 , η 0 ; ξ 0 , η 0 ) = 1.
This defines a characteristic boundary value problem whose solution R(ξ, η; ξ 0 , η 0 ) is known as the Riemann-Green function. 6 The knowledge of this function for a hyperbolic linear partial differential equation allows to construct the unique solution which is compatible with the given initial conditions. Concretely, using the defining properties of the Riemann-Green function (D5) in the identity (D3) after integrating at the deformed triangle with vertices P QR and taking into account the definitions (D2b)-(D2c), it is possible to isolate the value of the inhomogeneous solution u at the point P = (ξ 0 , η 0 ), i.e.
u(ξ 0 , η 0 ) = u h (ξ 0 , η 0 ) + u i (ξ 0 , η 0 ), where the first contribution solves the homogeneous equation respecting the initial conditions and is given by u h (ξ 0 , η 0 ) = 1 2 u(ξ 0 , η 1 )R(ξ 0 , η 1 ; ξ 0 , η 0 ) + 1 2 u(ξ 1 , η 0 )R(ξ 1 , η 0 ; ξ 0 , η 0 ) whereas the second contribution is a particular solution to the inhomogeneous equation An important remark is that in deriving the solution (D6) it has been implicitly assumed that the initial value arc QR is negatively inclined, with the future point P on its right, which is the more natural situation for characteristic coordinates. However, using a different notation it is also possible to have a positively inclined initial value arc with the future point on its left. In this case, the deformed triangle P QR is clockwise oriented and as a consequence the solution is again expressed as Eqs. (D6), except that the expression for the inhomogeneous solution (D6c) changes sign. Returning to our main objective, if we know the general homogeneous solution to Eq. (D4), the first contribution (D6b) is already considered there. Hence, we can rest in the results of the Riemann approach and just keep the expression (D6c)-or its negative counterpart in the positively inclined caseto obtain a particular solution of the inhomogeneous hyperbolic equation, as long as it is possible to build the corresponding Riemann-Green function. For our case of interest [the Euler-Darboux equation (C20)], the first step is to find the corresponding adjoint operator that results in a particular case of the extended Euler-Darboux operators (C14), specifically It only remains to impose the boundary conditions (D5b)-(D5d) after using, for example, the Poisson integral representation (C11) for the general solution of the Euler-Darboux equations. The resulting Riemann-Green function acquires the following form R(ξ, η; ξ 0 , η 0 ) = (ξ − η) α+β (ξ 0 − η) α (ξ − η 0 ) β 2 F 1 (α, β; 1; χ), (D9a) which can be derived in different ways [26,52,57], where 2 F 1 stands for the standard hypergeometric function.
Here, we provide an alternative derivation that allows us to prove at the same time that Eq. (D6c) is in fact a particular solution of the inhomogeneous Euler-Darboux equation (C20) if the above expression is considered as the corresponding Riemann-Green function. Our main motivation is to have a concrete expression available for the particular solution of the complexified version we study in subsection VII B 3, instead of providing a general solution to the initial value problem. Hence, we consider as a concrete initial value arc the straight line passing by the origin with slope −1; this is equivalent to choosing (ξ 1 , η 1 ) = (−η 0 , −ξ 0 ). In order to make our derivation more transparent and also to be compatible with the main text, we change the notation henceforth by relabeling the point equipped with the Riemann solution using now standard coordinates, i.e., P = (ξ, η). Consequently, the other two vertices of the triangle forming the integration region Ω are now given by Q = (−η, η) and R = (ξ, −ξ). Additionally, we use primes to denote the dummy coordinates in the integration. There are two elementary ways to sweep the triangle surface Ω in the integration (D6c): in the first, for each horizontal coordinate of the triangle, we sum on the corresponding vertical interval, and in the second we do exactly the opposite. We shall use a superposition of both: they are equivalent in the real domain covered in this appendix, but this is not the case if they are finally evaluated in complex limits, as we pretend to do in the main text. The advantage of this representation is that it naturally gives a real result in the complexified context. Taking the above into account and inspired by Eqs. (D6c) and (D9), we propose as a particular solution to the inhomogeneous Euler-Darboux equation (C20) the following expression where a priori no other assumption is made on F (χ) except that it is a function of the variable (D9b) written in the new notation. Applying now the Euler-Darboux operator (C1) to the above expression, we obtain the following identity after a careful evaluation E α,β (u i ) = F (0)f (ξ, η) where we denote the standard hypergeometric operator byĤ a,b;c (F ) ≡ χ(χ − 1)F ′′ + [(a + b + 1)χ − c]F ′ + abF.
(D11b) It is obvious now that if we choose the hitherto arbitrary function as the hypergeometric F (χ) = 2 F 1 (α, β; 1; χ), the integral contributions in the identity (D11) vanish. Considering now that 2 F 1 (α, β; 1; 0) = 1, it is proved that Eq. (D6c) [or, more specifically, Eq. (D10)] is in fact a particular solution of the inhomogeneous Euler-Darboux equation (C20) when the Riemann-Green function (D9) is considered. We emphasize that there are particular cases in the parameter space for which it is even possible to integrate the expression (D10) and arrive at a closed local form; see subsection VII B 3. Finally, for later use we write down the counterpart to Eq. (D10) when a positive-inclined initial arc is considered. The concrete initial value arc taken is the straight line passing by the origin with slope 1, which is equivalent to choosing (ξ 1 , η 1 ) = (η 0 , ξ 0 ) in the notation of the beginning of the appendix. In the notation of the last part, this means integrating on the triangle with vertices P = (ξ, η), Q = (η, η), and R = (ξ, ξ), giving The analog of the identity (D11) also follows in this case, which brings us to the same conclusion (D12). With this, we end our thorough review of results of the Euler-Darboux equations which are indispensable to justify the more important findings of the main text.