Schwarzschild quasi-normal modes of non-minimally coupled vector fields

We study perturbations of massive and massless vector fields on a Schwarzschild black-hole background, including a non-minimal coupling between the vector field and the curvature. The coupling is given by the Horndeski vector-tensor operator, which we show to be unique, also when the field is massive, provided that the vector has a vanishing background value. We determine the quasi-normal mode spectrum of the vector field, focusing on the fundamental mode of monopolar and dipolar perturbations of both even and odd parity, as a function of the mass of the field and the coupling constant controlling the non-minimal interaction. In the massless case, we also provide results for the first two overtones, showing in particular that the isospectrality between even and odd modes is broken by the non-minimal gravitational coupling. We also consider solutions to the mode equations corresponding to quasi-bound states and static configurations. Our results for quasi-bound states provide strong evidence for the stability of the spectrum, indicating the impossibility of a vectorization mechanism within our set-up. For static solutions, we analytically and numerically derive results for the electromagnetic susceptibilities (the spin-1 analogs of the tidal Love numbers), which we show to be non-zero in the presence of the non-minimal coupling.


Introduction
Black holes (BHs) are arguably among the most interesting objects in the universe. Their experimental agreement with the gravitational wave (GW) emission of binary systems [1] and with the imaging of the event horizon of a supermassive BH [2], promises to be but the first phase of a research era that is bound to culminate in a deeper understanding on the nature of BHs and gravity, as well as on numerous other related questions in astrophysics and fundamental particle physics.
As with most physical systems, a powerful way to probe BHs is by perturbing them and then see how they respond. While we cannot do this in the lab, such perturbed BHs are naturally produced by the merger of compact astrophysical objects. The details of how the post-merger BH relaxes toward equilibrium may in principle be measured through the GWs emitted during the process-the so-called ringdown phase. Quantitatively, the dynamics of the ringdown can be modeled by a superposition of quasi-normal modes (QNMs), whose characteristic frequencies are in one-to-one correspondence with the observable GW signal.
Although radiation in the form of GWs is a universal outcome of perturbing a BH, it need not be the only one. Indeed, a dramatic event such as a BH merger may reasonably be expected to excite other fields besides the metric, and these too will subsequently relax back to equilibrium via emission of the corresponding radiation. This radiation can again be described by QNMs, i.e. characterized in particular by dissipation caused by the presence of the BH horizon. More importantly, the matter fields' QNM spectra depend on the underlying spacetime, thus serving as an alternative probe of the BH. Furthermore, and crucially, the QNMs of a field encode physical information that is not directly available in the GW signal, namely about the coupling of the respective field with gravity and the equivalence principle.
In view of these considerations, we see at least two reasons that motivate the study of QNMs of matter fields in a BH background. The first concerns the fields themselves. As we have said, BH mergers are phenomena unlike anything we might achieve with Earth-based experiments. Thus, we may hope to make use of them as a way to test the existence of new particles, especially those that dominantly interact with the Standard-Model sector indirectly through gravity. The second reason which we have already alluded to regards the question of how fields couple to gravity. Establishing the existence of matter-gravity interactions beyond those dictated by the minimal coupling prescription is an exciting prospect that may in principle be achieved through the measurement of QNMs. In fact, as we will discuss later, QNMs offer a particularly clean signature of non-minimal gravitational interactions.
In this paper, we study QNMs of a massive vector field on a Schwarzschild BH background with a particular non-minimal coupling with gravity. Before describing our set-up in detail, let us briefly comment on the existing literature on the subject of vector-field QNMs in BH spacetimes. The study of massless, minimally coupled vector fields in four dimensions and with flat asymptotics dates back to the work of Chandrasekhar [27]. The Proca equation for a massive vector field and the corresponding QNM spectra have been investigated in [28][29][30][31] for a Schwarzschild(-AdS) BH and only recently in [32][33][34] for a Kerr BH.
Here we go beyond previous studies of spin-1 particles by considering the most general Lagrangian of a vector field A µ subject to the following assumptions: (i) The Lagrangian is quadratic in the vector field. This follows from our aim to investigate linear perturbations about vacuum solutions of general relativity (GR), specifically the Schwarzschild metric. Generically, this implies that the vector field must vanish at the background level, and therefore it is sufficient to focus on a quadratic theory for the purpose of deriving the QNM spectrum.
While this is true generically, we should remark that there exist vector-tensor theories that admit so-called "stealth" BH solutions, i.e. solutions that coincide with vacuum GR solutions in spite of having a non-trivial vector field background profile [35]. Our analysis therefore does not encompass this case.
A corollary of this premise is that metric and vector perturbations are decoupled at linear order. The QNM spectrum of GWs is thus exactly the one derived in GR [7,27] and so may be ignored.
(ii) The theory describes precisely five dynamical degrees of freedom, i.e. two in the metric and three in the vector field (or two in the case of a massless spin-1 field, which we will treat as a special case). In other words, we demand the absence of additional propagating modes associated to a loss of constraints or to higher-order equations of motion.
For a generic spacetime background, this Lagrangian extends the Proca theory by the addition of two non-minimal coupling operators. Unsurprisingly, these operators are precisely those obtained from linearizing the Lagrangian of the Generalized Proca theory of a self-interacting massive spin-1 field [36,37]. Our derivation thus serves as a proof of the uniqueness of Generalized Proca theory at the level of linear perturbations about the trivial state A µ = 0. For a Ricci-flat background the theory further simplifies, leaving only one non-minimal coupling operator, with F µν the vector field strength and G 6 a coupling constant. Our main objective in this paper is to numerically derive the QNM spectrum of vector field perturbations for a range of values of the parameter G 6 and the bare mass µ of the field. Interestingly, for a given BH mass, G 6 is restricted to a window of values given by where r g ≡ 2GM is the Schwarzschild radius (with G the Newton coupling and M the BH mass). This criterion follows from the requirement of stability (due to ghost and/or gradient instabilities) of the BH under perturbations of the vector field in the localized approximation, i.e. in the limit where the size of the perturbation is much shorter than the typical length scale characterizing the background variation [38,39]. Related to the question of stability of BH spacetimes under perturbations of generalized vector fields, one may ask if the criterion (1.2) is not only necessary but also sufficient for ensuring stability. While we plan to address this with exhaustivity in a dedicated work, here we provide evidence that this is indeed the case for a Schwarzschild BH. Our claim is based on the analysis of quasi-bound states of the vector field, that is solutions of the generalized Proca equation which decay at spatial infinity. Like QNMs, quasi-bound states have an associated spectrum of complex frequencies, which are of particular interest as they may be used to diagnose the presence of instabilities: A quasi-bound state frequency with positive imaginary part signals an exponentially growing mode and thus an unstable system, at least within the linearized regime. It is worth remarking that the same judgment cannot be made based on the QNM spectrum because, as we will review, the imaginary part of a QNM frequency must be negative by the definition of a QNM.
The principal result of this first study of quasi-bound states of a non-minimally coupled vector field is that the fundamental frequency mode for each degree of freedom of the field has a negative imaginary part within the numerically accessible part of the range given by Eq. (1.2). However, as the computational cost of our numerical routine increases as one approaches the bounds in (1.2), we are unable to numerically access values of the coupling G 6 arbitrarily close to the critical points. We partially address this shortcoming by providing an analytical argument, valid for a subset of the spectrum, which shows that quasi-bound states are stable whenever G 6 is within but arbitrarily close to the stability bounds.
We have mentioned that the QNM spectrum of matter fields may serve as a powerful tool to test the minimal coupling paradigm dictating the form of matter-gravity interactions. This question is of fundamental importance, so it behooves us to understand which signatures of a non-minimal coupling operator for a given field might be clean and robust enough so as to be potentially detectable. Perturbations of a massless field are arguably one such probe, since minimally coupled massless fields on BH spacetimes in GR are known to be very special, at least due to two properties: isospectrality of their QNM spectra and vanishing linear response coefficients in the static limit.
Isospectrality refers to the equivalence of the QNM spectra of parity-even and parityodd perturbations [27,40]. The property is featured by massless fields in four dimensions with flat or de Sitter asymptotics, at least for spins s = 0, 1, 2. 1 Isospectrality is however known to fail in higher dimensions [43], for asymptotically anti-de Sitter BHs [13,44], for massive fields [30,45], and for BHs in non-linear electrodynamics [46,47] or in the presence of higher-curvature corrections [16,18]. To our knowledge, the breaking of isospectrality due to non-minimal couplings has not been systematically addressed, although it is known to occur for certain couplings of scalar fields [21,[23][24][25][26]. Here we fill the gap of spin s = 1 by showing through numerical results that the parity-even and -odd spectra of a massless vector field with the non-minimal coupling of eq. (1.1) are indeed distinct.
Although QNMs are the main focus of our work, static perturbations are also interesting in that they define the static response coefficients associated to a given field. For massless spin-2 perturbations the response coefficients physically encode the tidal deformability of the BH and are known as Love numbers [48], see also [49][50][51]. For a massless spin-1 probe field they may be interpreted as the electromagnetic susceptibilities of the field in a BH background. It is a remarkable and well-known property that the static response coefficients of massless fields of spin s = 0, 1, 2 exactly vanish for four-dimensional BHs in GR [52][53][54][55][56][57][58]. The property is however absent in higher dimensions [59,60] as well as for BHs in beyond-GR theories [16,61,62]. In addition, and similarly to isospectrality, the vanishing of Love numbers and susceptibility coefficients is not expected to hold in the presence of non-minimal couplings, although again we are not aware of any exhaustive analyses (see [61] for results in some particular models). Here, we compute the electric and magnetic susceptibilities of dipolar perturbations of a massless vector field as functions of the coupling G 6 in eq. (1.1), and show that they are non-vanishing in agreement with expectations.
We now give an outline of the paper's contents: In Sec. 2, we describe our set-up, including (i) our uniqueness argument for the non-minimal coupling, (ii) the decomposition of the vector field in spherical harmonics, and (iii) the definition of QNMs according to the boundary conditions for the mode equations. In Sec. 3, we present our main results, namely the calculation of the QNM spectra for each mode of the vector field and for a range of values of the coupling G 6 and mass µ. The spectra of a massless field and the breaking of isospectrality are treated as a special case. In Sec. 4, we consider quasi-bound states. This provides evidence for the stability of the system under consideration beyond the local approximation. In Sec. 5, we consider static perturbations and, focusing on a massless field and dipolar modes, derive the electric and magnetic susceptibilities as functions of G 6 . We discuss our results and give some final remarks in Sec. 6. In Appendix A, we provide details of the numerical method used in our calculations.

Non-minimally coupled Proca field
Our study will focus on linear perturbations of a massive vector field A µ about a GR background solution. The background state of the vector field is the trivial one, A µ = 0, as per our definition of a GR solution, i.e. one with vanishing vector hair. It therefore suffices to focus on Lagrangians that are precisely quadratic in the field A µ , while the dependence on the metric tensor is in principle arbitrary. Note that, a priori, we make no restriction on the number of derivatives acting on A µ .
We will additionally require that the theory describe exactly five degrees of freedomtwo in the metric and three in the vector field-so as to avoid Ostrogradsky-type ghosts on all backgrounds. A sufficient condition to achieve this is to demand that the equations of motion of the Stückelberg formulation of the theory be of second order in derivatives. This condition is however not a priori necessary, as it may occur that the theory possess the correct number of constraints even in the presence of higher derivatives in the field equations [63]. We shall nevertheless disregard this possibility here and focus on the simpler set-up with second-order equations of motion.
Our claim is that the most general four-dimensional Lagrangian subject to these assumptions is given by where M Pl is the Planck mass, µ is the mass of the vector field, and G 4,X and G 6 are coupling constants. The notation chosen for the latter two coefficients is explained by the connection between the Lagrangian (2.1) and the Generalized Proca theory. As mentioned in the introduction, eq. (2.1) may be obtained upon linearizing the Generalized Proca Lagrangian about the trivial vector background A µ = 0. 2 In particular, the operators multiplying G 6 in the second line (which may be more compactly written in terms of the dual Riemann tensor) will be recognized as the unique extension, as demonstrated by Horndeski [64], of the standard Einstein-Maxwell theory, here restricted to quadratic order. In this paper we confine our attention to a background given by the Schwarzschild metric and neglect the backreaction of the vector field on the geometry. This assumption is valid at linear order in perturbation theory since, as we remarked, metric and vector fluctuations do not couple at this order. The generalized Proca equation for a Ricci-flat spacetime reduces to In this set-up, we are therefore left with two dimensionless parameters: µr g and g 6 ≡ G 6 /r 2 g (with r g the Schwarzschild radius). Observe that the Lorenz constraint, follows as a consequence of eq. (2.2) whenever µ = 0. In the massless case, we shall instead impose a different constraint as a gauge condition. As mentioned in the introduction, eq. (2.2) features pathological solutions (ghosts and/or gradient-unstable modes) unless the coefficient g 6 is confined to the range [38,39] While this result was obtained from an analysis of localized perturbations, these bounds on g 6 will be seen to translate into the statement that the mode functions of the vector field should be insensitive to additional poles appearing in the equations of motion. In terms of the Schwarzschild radial coordinate r, these poles are given by with r + ≡ g 1/3 6 r g and r − ≡ (−2g 6 ) 1/3 r g . Demanding that these poles be hidden inside the event horizon then yields (2.4). Thus, although this range was originally derived from different considerations, it has the important implication that the equations will allow for consistent QNM solutions, which at least generically would not be possible if one had poles in the physical domain r > r g . 3 2 The Generalized Proca Lagrangian contains the functions G4(X) and G6(X) (among others), with X ≡ − 1 2 A µ Aµ. Our coupling constants G4,X and G6 correspond respectively to G 4 (0) and G6(0), which are finite by our assumption that Aµ = 0 is a well-defined state. 3 QNMs are by definition everywhere regular and with fixed boundary conditions. The presence of a pole would impose an additional matching condition and thus an overdetermined system for the QNM frequency and the amplitude of the QNM function. Such a system will generically have no solution. The same remark, of course, also applies to quasi-bound states.

Uniqueness
Generalized Proca theory was constructed as the most general model which reproduces the (shift-symmetric) scalar Horndeski theory in the so-called decoupling limit where the longitudinal mode of the vector field becomes a dynamical scalar [65,66]. 4 As such, the theory is perfectly general, given the assumptions of its construction, on flat spacetime. The uniqueness of Generalized Proca is however not immediate when the coupling with gravity is taken into account, since the covariantization of the decoupling limit theory need not match, term by term, that of the full theory. In particular, one cannot a priori disregard non-minimal couplings to the curvature tensor beyond those obtained in [37] (see also [72]), as the latter were derived as "counterterms" to cancel the pathological operators that appear upon minimal covariantization. Here, we provide a sketch of the proof of the uniqueness of the Lagrangian (2.1); a detailed proof will be given in a dedicated work where we analyze the general problem without assuming linearity in the vector field.
To reiterate the problem, we seek the most general Lagrangian for a vector field A µ and metric tensor g µν subject to the assumptions of (i) general covariance, (ii) quadratic order in the vector field, and (iii) second-order equations of motion for all the fields in the Stückelberg formulation. Note that we make no assumption on the derivative order of the fields at the level of the Lagrangian.
In the Stückelberg formulation of the theory, the Lagrangian is a functional of the fields (g µν , A µ , φ) and is invariant under diffeomorphism and U (1) gauge symmetries. The latter property implies that the vector and Stückelberg scalar can only appear through the invariants F µν and D µ φ ≡ ∇ µ φ + µA µ (here µ is the mass of the Proca field). The main proposition is that these two building blocks cannot couple with each other in the Lagrangian. To establish this one notes that covariant derivatives of D µ φ may be chosen as fully symmetrized without loss of generality. Indeed, any mixed-symmetric or antisymmetric projection of ∇ µ 1 · · · ∇ µ n−1 D µn φ can be traded by F µν (and derivatives thereof) and/or curvature tensors contracted with fully symmetrized derivatives of D µ φ. Since derivatives of F µν cannot be made fully symmetric, it follows that they cannot be contracted with the tensor ∇ (µ 1 · · · ∇ µ n−1 D µn) φ. An exception to this is the divergence of the field strength, ∇ µ F µν , and its derivatives; for instance, ∇ µ F µν D ν φ is a valid operator that seemingly contradicts our claim. However, the divergence ∇ µ F µν may in principle be solved for algebraically from the vector field equation of motion, implying that any instance of this term in the Lagrangian may be eliminated through a field redefinition.
The Lagrangian is therefore "separable" in the building blocks F µν and D µ φ, which may then be analyzed independently. The operators involving only F µν and the metric constitute purely vector-tensor gauge invariant terms, hence they satisfy the assumptions of the Horndeski theorem for Einstein-Maxwell theory [64], with the known quadratic-order result 5 L ⊃ 1 4 Other prescriptions for constructing vector-tensor theories have been considered in the literature [67][68][69][70], leading to various extensions of Generalized Proca. See also [71] for an effective field theory approach. 5 The double-dual Riemann tensor is defined as R µνρσ ≡ 1 together with the standard Maxwell, Einstein-Hilbert and cosmological constant terms. The remaining operators in the Lagrangian must then all be expressible in terms of D µ φ and covariant derivatives of this invariant. Because the conditions we are imposing on the equations of motion must hold for all field configurations, they must hold, in particular, when A µ = 0. But in this case D µ φ → ∇ µ φ and we have precisely the assumptions of the Horndeski theorem for scalar-tensor theory [73] (with the extra condition that φ may not appear without derivative), with the known quadratic-order result together with the standard scalar kinetic term. In the general case with A µ = 0, we know that the scalar field derivative must appear "covariantized" in D µ φ, so that the result correctly reproduces the Generalized Proca term upon setting unitary gauge φ = 0. This concludes our derivation of the Lagrangian (2.1), independently of its relation with the non-linear Generalized Proca theory. The implication is that any consistent extension or alternative to Generalized Proca must reduce to (2.1) when expanded at quadratic order about the vacuum A µ = 0, provided the theory admits this state.

Decomposition in vector spherical harmonics
We consider the exterior of a Schwarzschild BH spacetime with line element where r g = 2GM , G is the Newton coupling, and M is the mass of the BH. Given the background symmetries, the equation of motion for the vector field is separable after expanding in spherical harmonics, where, in our convention, the vector spherical harmonics are defined as in terms of the standard scalar spherical harmonics Y m (θ, φ). Under a parity transformation, (θ → π − θ, φ → π + φ), the functions Z  and may be analyzed separately. The vector spherical harmonics satisfy the orthonormality condition where M µν Z = diag 1/f 2 , f 2 , ( + 1)/r 2 , ( + 1)/(r 2 sin 2 θ) , which is used to factor out the angular dependence in the equation of motion.
In the following, we suppress the supersripts and m in the mode functions u m i and denote partial derivatives with respect to t and r respectively with dots and primes. We also introduce the operator with r * being the tortoise coordinate defined by dr * = f −1 dr. In the next subsection we provide the mode equations for the dynamical degrees of freedom. The case of a massless vector field requires a separate analysis, which is done in the following subsection. Readers interested only in the relevant equations and results may consult Tab. 1 for reference.

Mode equations
In order to eliminate the non-dynamical variables we make use of the Lorenz constraint, eq. (2.3), which reduces tou upon substituting the expansion in eq. (2.9). Here, and in the following, we denote tderivatives with a dot and r-derivatives with a prime.
In the case of monopole ( = 0) perturbations, u 3 and u 4 are absent in the spherical harmonic expansion. Using constraint (2.16) we can further eliminate u 1 in favor of u 2 ≡ u M , with the resulting equation For axial perturbations with ≥ 1, it is convenient to define u − ≡ P 1/2 + u 4 , which produces (2.18) For the polar modes with ≥ 1 we again eliminate u 1 using the constraint (2.16), obtaining two coupled equations for the variables u 2 and u 3 , As anticipated in table 1, the monopole mode is only sensitive to the pole P − , axial modes are only sensitive to the pole P + , while polar modes with ≥ 1 are affected by both.
We remind the readers that the parameter g 6 (implicit in the above equations, cf. (2.5)) is restricted to lie in the stability range (2.4), so that the poles P ± never vanish in the physical domain r > r g . Nevertheless, the observation is pertinent as we shall be interested in exploring values of g 6 close to the bounds.

Massless case
When the bare mass µ vanishes, the Lagrangian (2.1) is gauge invariant and the identification of the dynamical degrees of freedom requires a separate analysis. We will use the gauge freedom to set u 1 = 0 as was done in Ref. [30]. Note that this is a complete gauge fixing for perturbations compactly supported in space and time. For = 0, both u 3 and u 4 are again absent, while the generalized Proca equation implies that u 2 = 0, indicating as expected that there is no dynamical monopole mode. For the higher multipoles with ≥ 1 we introduce in terms of which the parity-even part of the equation of motion can be cast as so that there is a single polar mode (for each , m) in the massless case. As for the axial mode, being gauge invariant, one can directly set µ = 0 in eq. (2.18) to obtain the corresponding equation.

Boundary conditions
We seek solutions to the mode equations in the frequency domain, where they assume the form for the modes u M , u − , u 2,3 and u 0 (cf. table 1), and remembering that the functional V couples both modes u 2,3 in the polar sector. The frequency ω is in general complex, assuming without loss of generality a positive real part. (If u solves the mode equation for some frequency with Re ω > 0, then u * solves the conjugate equation with Re ω < 0.) Presently, we further assume Im ω < 0, deferring a discussion of the opposite case to Sec. 4. The BH horizon serves as a causal boundary admitting only ingoing modes, hence the physical boundary condition is at the horizon, i.e. as r * → −∞.
At spatial infinity, r * → +∞, we have V(u, r) µ 2 u for every mode. The general asymptotic solution at spatial infinity is therefore By definition, QNM solutions correspond to purely outgoing waves at infinity, i.e. with c in = 0. 6 Having fixed boundary conditions at both the event horizon and at spatial infinity, we are left with an eigenvalue problem with a discrete set of QNM solutions characterized by a spectrum of frequencies {ω n } ∞ n=0 .

Quasi-normal modes: numerical results
Recall that our mode equations depend on the two parameters µ and g 6 . The standard Proca theory corresponds to g 6 = 0, whose QNM spectra on a Schwarzschild background were studied in Ref. [30]. Our main aim here is the extension of the analysis to non-zero values of g 6 within the stability range (2.4), sampling also over a range of mass values µ. We restrict our attention to the fundamental QNM frequency (n = 0) and lowest multipoles = 0, 1, except in the massless field case for which we present results also for the first and second overtones (n = 1, 2) of the dipole modes.
We numerically solve the mode equations using a spectral or collocation method with Chebyshev interpolation, using up to N = 80 collocation points to ensure converged results. In essence, the method turns a differential boundary-value problem into a non-linear eigenvalue problem with finite-dimensional matrix. A brief summary of the approach is given in 6 To see explicitly that e √ µ 2 −ω 2 r * is an outgoing wave, note that we choose the convention for the square root such that Re µ 2 − ω 2 > 0, which implies that sign(Im µ 2 − ω 2 ) = −sign(Im ω) = +1. Appendix A; the reader may find a succinct but more general exposition in [33,74], which also provides references to the relevant mathematical literature. Before proceeding, a word about terminology. The polar sector contains two degrees of freedom for each ≥ 1, hence two independent QNM spectra. We will refer to these modes as "scalar" and "vector", following [30]. The rationale behind these names is that, in the massless limit and with g 6 = 0, the polar mode equations match the form of the Regge-Wheeler (RW) equations for massless scalar and vector fields, in agreement with the Goldstone boson equivalence theorem. To see this explicitly, set µ = 0 and g 6 = 0, and introduce where if y 2 = 0, then y 3 is a pure gauge degree of freedom while y 2 satisfies the massless vector RW equation. These considerations can be generalized to the case with g 6 = 0, with the same conclusion: in the massless limit, the polar spectrum can be divided into two classes, one corresponding to a massless scalar field and another corresponding to a vector gauge field.
In Figs. 1-4, we present the results for the = 0, 1 fundamental (n = 0) QNMs, in case of non-vanishing mass µ = 0. The behaviour with 0 < µr g < 0.5 and −1/2 < g 6 < 1 is mapped out in terms of contour plots. To reveal pole-induced behaviour, we also plot the g 6behaviour at fixed exemplary µ. The chosen range for the vector field mass Îĳ is motivated by the fact that one expects interesting physical effects when the Compton wavelength of the field is comparable or larger than the size of the BH, i.e. µr g 1. This can also be understood more mathematically from the fact that the norm of the QNM frequency can typically be estimated as |ω| 2 ∼ V max , where V max is the height of the centrifugal potential barrier [75]. Now, for the QNM function to have the required wave-like behavior at spatial infinity, one also requires |ω| > µ. It follows that there will be no QNMs if µ 2 is greater than the height of the centrifugal barrier, i.e. if O(1) for the lower multipoles of most physical interest. 7 The individual results can be understood by revisiting Tab. 1. Each of the modes diverges at the critical values g 6 = −1/2 and/or g 6 = +1 iff the respective perturbation equation is affected by the corresponding pole P − (P + ). The = 0 monopole mode, cf. Fig. 1, is affected by P − only. The = 1 axial mode, cf. Fig. 2, is affected by P + only. Finally, the two coupled polar multipole modes are affected by one pole each: the scalar mode, cf. Fig. 3, by P − ; the vector mode, cf. Fig. 4, by P + . While not presenting respective results, we expect the same pole-induced behaviour to persist for all higher > 1 modes.
In case of vanishing mass µ = 0, the perturbations reduce to one axial and one polar mode only, cf. Sec. 2.4. For minimal coupling to the background metric, i.e., for g 6 = 0, the two QNM spectra are known to be isospectral, i.e., the axial and polar spectrum agree. Indeed, the respective perturbation equations (2.24) and (2.18) (with µ = 0) become redundant for g 6 = 0, i.e., for P ± = 1. For any g 6 = 0, isospectrality is broken. We verify this explicitly by presenting the n = 0, 1, 2 massless modes in Fig. 5.
As in the massive case, the observed behaviour close to g 6 = −1/2 and/or g 6 = 1 is determined by the poles in the respective perturbations equations. The axial massless mode, cf. eq. (2.18), is affected by P + only. The polar massless mode, cf. eq. (2.24), is affected by both poles. The onset of this pole-induced behaviour can be explicitly seen for  the n = 0 mode. For the n > 1 modes it becomes numerically challenging to resolve.
In fact, numerical convergence of the spectral methods worsens considerably with growing n. This can intuitively be understood as follows. The number of oscillations in r increases with n. The modes thus become more and more challenging to resolve with spectral methods based on a fixed number of collation points. The n > 0 results presented in Fig. 5, thus present the most challenging numerics of this work. Hence, we explicitly present convergence plots for exemplary points in App. A.
The behavior of the QNM spectrum for small values of g 6 is worth remarking. As one can glean from Fig. 5 (although we have also verified it from the numerical data), for each n the polar and axial QNM frequencies display a symmetry in their g 6 -dependence at linear order, exhibiting the same slope (within numerical precision) but with opposite sign. This feature may hint at the existence of an electromagnetic duality for small but non-zero values of G 6 . 8 We will encounter a similar phenomenon when we consider the electromagnetic susceptibilities in Sec. 5.
Finally, we comment on the observed crossing of the axial n = 1 and n = 2 imaginary parts close to g 6 = 1. We are not aware of other examples of such a crossing of imaginary parts. While we find no indication for convergence issues in the applied spectral methods, a confirmation of this result by independent numerical techniques would be welcome. . For any non-vanishing g 6 = 0 isospectrality is broken. We show the fundamental (n = 0) as well as first (n = 1) and second (n = 2) overtone in increasingly lighter shading. Where the curves end, spectral methods with N = 80 are found to be insufficient to ensure proper convergence. Exemplary convergence plots (for the points marked with triangles on the g 6 -axis) are presented in App. A.

Quasi-bound states
Quasi-bound state solutions to the mode equations are defined by boundary conditions corresponding to an ingoing wave at the event horizon and a vanishing amplitude at spatial infinity. The latter requirement selects c out = 0 in (2.27), oppositely to the case of QNMs. Importantly, the bound state behavior u ∼ e − √ µ 2 −ω 2 r * holds regardless of the sign of Im ω.
The last remark is apposite given our interest in establishing whether the theory admits unstable solutions with Im ω > 0, even if the coupling g 6 lies in the range (2.4) in which localized perturbations are stable. The latter condition is necessary for consistency, as it has been shown that localized modes must be either stable or else suffer from ghostor gradient-type instabilities, while tachyon-type unstable solutions cannot occur [39]. The caveat to this statement is that tachyonic solutions cannot be fully diagnosed in the localized approximation, which is by definition oblivious to modes of physical size comparable or larger than the length scales of the background. In other words, we would like to assess if global solutions could undergo instabilities in the regime where the theory is free from pathologies related to ghosts and negative-gradient modes.
Here we provide strong evidence that tachyonic quasi-bound state solutions for a massive vector field cannot occur on a Schwarzschild BH background. Our first argument in support of this claim is given by the numerical results, presented in Sec. 4.1, for the fundamental (n = 0) bound state solution for each of the lowest multipole modes of the field ( = 0, 1), sampling over a range of values of the parameter g 6 and a few values of the mass µ. While this certainly does not constitute a full proof, one naturally expects tachyon modes to appear for the lowest values of n and if they exist at all. 9 Indeed, tachyonic instabilities should, by definition, eventually disappear as the typical radial and angular wavelengths of the solutions (characterized respectively by n and ) become short enough.
A more critical loophole in this numerics-based argument is our inability to access values of g 6 arbitrarily close to the stability bounds (2.4), as our numerical routine becomes increasingly less efficient as we approach those values. This is important in view of the expectation (suggested also by the numerical results) that quasi-bound state frequencies will differ the most from their values in standard Proca theory precisely near the critical g 6 points. Fortunately we can patch this issue by means of an analytical proof which shows that the imaginary part of the frequency cannot be positive. This argument is also not a complete one, however, first because it does not apply to the polar modes with ≥ 1, and second because it assumes that g 6 lies sufficiently close to either of the critical points. These caveats notwithstanding, the argument is otherwise general, valid for any mass µ (with µ 2 > 0) and for all multipoles , m. We describe the argument and its application to the monopole and axial modes in Sec. 4.2.

Numerical results
As for the case of quasi-normal modes, cf. Sec. 3 and App. A, we numerically solve the quasibound state mode equations via spectral methods with Chebyshev interpolation. We restrict the analysis to the fundamental monopole ( = 0) and lowest multipole ( = 1) modes. We make sure that all of the following results are converged to at least 5% accuracy. (For most parameter values the accuracy is much higher, cf. App. A for exemplary convergence plots in the QNM case.) In Fig. 6, we summarize the behavior of the fundamental = 0, 1 quasi-bound states in the complex-frequency plane. 10 To do so, we show the quasi-bound states for −0.4 < g 6 < 0.9 and representative µ × r g = 1, 2/3, 1/2. With µ → 0, these all converge to Re[ω/µ] → 1 and 0 > Im[ω/µ] → 0. This holds for any constant g 6 , at least in the investigated range. For the axial mode, cf. upper-right panel in Fig. 6, we find Im[ω/µ] −→ 0, for all investigated values of µ. We find no indications for the onset of such scaling for the other modes, at least within the investigated range.
In Fig. 7, we exemplify the power-law behaviour that all modes exhibit as µ → 0. Here, we choose to present results at a representative value of g 6 = 1/2 only. The behaviour closely resembles the one previously found for g 6 = 0 [30]. 9 For instance, in the case of a massive spin-2 field on a BH background it is the n = 0, = 0 mode, and only this mode, which is unstable for a certain range of parameters [42,45]. 10 With some abuse of terminology, we refer to the polar = 1 modes as "scalar" and "vector" as we did for QNMs, although in reality quasi-bound states cease to exist in the massless limit and therefore the Goldstone boson equivalence limit is not meaningful.  Figure 6. Behaviour of the fundamental = 0, 1 bound states with changing µ and g 6 . For each mode, we present the parametric curves mapped out by −0.4 < g 6 < 0.9 (cf. color legend on the right) for three different values of µ × r g = 1, 2/3, 1/2. The dots indicate the respective value for g 6 = 0.
To summarize, we find no indication for the presence of an unstable mode (i.e., one with Im[ω] > 0). Whenever modes scale towards Im[ω] = 0, we have identified the respective power-law scaling. We view this as strong numerical evidence for the absence of unstable quasibound states.

Integral formula
Next we turn to the analytical proof of the fact that Im ω < 0 in our set-up. The method is essentially the one put forth in Ref. [76] in the context of asymptotically anti-de Sitter BHs (see also Ref. [44] for further applications). The interesting observation is that the argument also applies to bound state perturbations of asymptotically flat BHs, albeit with some differences.
We consider eq. (2.25) in the case of a single ODE, so that V(u, r) ≡ V (r)u. We introduce the redefined mode function v ≡ e iωr * u. The boundary conditions imply that v approaches a constant, v + , at the horizon and that it decays exponentially at spatial We multiply through by v * and integrate, Note that each term in this integral gives a finite result thanks to the exponential decay of v (while V is non-singular by assumption). The first term can be integrated by parts, noting that the boundary term vanishes, Taking the difference of this equation with its complex conjugate we get which can be plugged back in (4.3) to produce We see that if the potential function V were positive definite, then we would immediately infer that Im ω < 0 and conclude the proof. However V is not positive definite in the equations within our set-up. Nevertheless, we can prove that, for each mode, its contribution to the integral is indeed non-negative whenever g 6 is sufficiently close to the critical points, i.e., for values such that the poles P ± coincide with the event horizon. The effective potential of the monopole mode is given by and it is easy to see that V M is not positive definite for all values of µ 2 and g 6 . However, as we are interested in the case when g 6 lies near the bound g 6 = −1/2, we define ≡ g 6 + 1/2 and isolate the leading-order contribution to the integral (4.5) in an expansion in small , i.e., This integral is manifestly positive. Similarly, for the axial modes the effective potential reads which is also not positive definite for all µ, and g 6 . The relevant pole is now g 6 = 1, so we let ≡ 1 − g 6 and evaluate the integral to leading order in the limit of small , 9) and the result is likewise manifestly positive. This establishes that Im ω < 0 for quasi-bound state perturbations corresponding to monopole and axial modes. For the polar modes with ≥ 1 the argument does not readily apply since in this case one has to deal with a system of coupled equations and with additional terms proportional to derivatives of the mode functions, cf. eq. (2.20). Even though an analogue of eq. (4.5) can be straightforwardly derived, we have been unable to find a bound for the integral of the resulting effective potential. Nevertheless, we see no reason why polar perturbations should behave qualitatively different from the rest of the spectrum, and the numerical results certainly seem to confirm this. Moreover, as we remarked previously, the expectation is that unstable modes, if they exist, should manifest themselves at the lower end of the multipole ladder. Given our proof of the stability of monopole fluctuations, we take these combined results as strong evidence for the absence of instabilities in the whole quasi-bound state spectrum and the whole range of allowed values of the non-minimal coupling g 6 .

Electromagnetic susceptibilities
Static response coefficients characterize the change of a system under an external timeindependent field. For a gravitational field, these coefficients correspond to the tidal Love numbers, which are in principle directly measurable through gravitational wave observations, e.g. of binary systems. For a U (1) gauge field the response coefficients are the electric and magnetic susceptibilities defining the polarizability of the object (in analogy with electromagnetism, although the field of course need not be the Standard Model photon).
As mentioned in the introduction, a remarkable property of four-dimensional BHs in GR is that they do not polarize under the effects of a Maxwell-type field. Yet the expectation is that this attribute will be broken in more general set-ups, in particular if the external U (1) field contains additional interactions, either with itself or with the spacetime metric.
Within our set-up of a Schwarzschild BH and in linear response theory, we have seen that it is only the Horndeski non-minimal coupling operator, eq. (1.1), which can contribute to beyond-GR effects without introducing additional degrees of freedom. The question is then whether the electromagnetic susceptibilities are indeed non-vanishing when the coupling G 6 is non-zero. Here we confirm that they are non-vanishing, focusing for simplicity in the case of dipolar perturbations.

Boundary expansion
We consider the mode equations in the gauge invariant case, i.e. eq. (2.24) for the polar or "electric" field and eq. (2.18) (with µ = 0) for the axial or "magnetic" field, setting ω = 0 as we are interested in the static limit.
For each equation, only one linear combination of the two independent solutions is regular at the event horizon. Demanding regularity thus fixes one integration constant, while the other remains arbitrary, simply setting the overall amplitude of the mode function. Then, modulo this overall constant, the solution at spatial infinity is fully determined, and is in general given by a sum of two modes, one which grows and one which decays with the radius r, i.e., The leading coefficient of the growing mode, c ext , is interpreted as the strength of the applied field, while that of the decaying mode, c resp , gives the corresponding response of the system. Their ratio, k ≡ c resp c ext , defines the linear susceptibility of the system for the given external field. There are two remarks to keep in mind about the structure of the boundary expansion in eq. (5.1), both related to the fact that the expansion is a Frobenius series. The first is that the series multiplying r +1 in the growing mode may in general contain logarithmic terms. However, in four dimensions these are always subleading and do not affect the definition in (5.2). 11 The second observation is that, because is an integer, the split between the growing and decaying modes is potentially ambiguous as they contain the same powers of r after some order [77,78]. Various ways to deal with this issue have been proposed in the literature [51,59,79,80], although in our case it will suffice to simply define the growing mode such that it does not contain the power r − (which is not to say that the series terminates, since all subsequent powers may a priori be present). This prescription makes the susceptibility k unambiguous and is physically justified by the fact that k so defined is an observable enjoying the property we seek: it vanishes in the absence of non-minimal coupling but is otherwise non-zero, as we now show.
For illustration, let us focus on the dipole modes ( = 1). Interestingly, we find that there is no logarithmic term in this case. However, unlike in the ordinary Maxwell set-up, the series for the growing mode does not terminate. Explicitly, for the first few terms we find 4c resp − g 6 (10 − 11g 6 )c ext 8 3(2c resp + 5g 6 c ext ) 10 4c resp + g 6 (10 − g 6 )c ext 8 (5.4) respectively for the electric and magnetic components.
Although the mode equations do not seem to admit an exact solution, they may straightforwardly be solved iteratively by expanding in powers of the coupling g 6 . 12 After selecting the regular solution in each case, as explained above, we are then able to infer k. We find, up to order O(g 4 6 ), respectively for the electric and magnetic susceptibilities. The expressions in (5.5) confirm the vanishing of the susceptibility in standard Maxwell theory, i.e. with g 6 = 0. For small but non-zero g 6 we find instead the expected dependence k = O(g 6 ). We observe that the electric and magnetic susceptibilities are equal, up to a sign, at linear order in g 6 . We recall that the same phenomenon was observed for the QNM spectrum in the massless (i.e., gauge-invariant) case, cf. Fig. 5. Also remarkable is that the electric susceptibility does not appear to receive non-linear corrections. These results are intriguing and clearly beg for a deeper physical understanding. We hope to come back to this question in future work.

Numerical results
Having understood analytically the polarizability properties of a BH in the approximation of small g 6 , we now turn to the exact results derived numerically using a shooting method. We have computed the solutions of the mode equations in the vicinity of the BH horizon by expanding in powers of (r − r g ), up to order four. We then evaluate the function at some small (r − r g ), which is used as initial condition to integrate numerically up to some large radius. The result is then matched to the boundary series discussed in the previous subsection, which we expand up to order r −8 , so as to obtain c ext and c resp , and hence the susceptibility k. We have checked that the results are robust against changes in the initial and matching radii as well as in the order at which we terminate the series ansatze. 12 We recall that g6 enters in the equations through the combinations P± = 1 − r 3 ± /r 3 (cf. (2.5)), where r± < rg < r and r 3 ± ∝ g6. Therefore, for any r in the physical domain, the mode equations are indeed analytic at g6 = 0 and the expansion in Taylor series is justified. The results for the electric and magnetic susceptibilities are shown in Fig. 8, plotted as functions of the coupling g 6 within the stability range eq. (2.4). The plots also show the comparison with the approximate analytical behaviors in eq. (5.5), which are indeed in perfect agreement with the numerical results. For the electric susceptibility we confirm the interesting outcome that the linear truncation in eq. (5.5) appears to be exact, within our numerical precision, even as we get very close to the critical values of g 6 (we can reliably compute k for |g 6 − g crit 6 | 10 −3 ). In contrast, the magnetic susceptibility shows a clear departure from the polynomial approximation for sizable values of g 6 , as one would generically expect 13 .

Discussion
The aim of this paper was to initiate the study of global solutions for massive vector fields non-minimally coupled to gravity in the linear approximation about GR backgrounds. We focused on the simplest but physically important case of a Schwarzschild BH background, and restricted our attention to a single non-minimal coupling operator, namely the Horndeski term given in eq. (1.1). In spite of the simplicity of the model under consideration, we showed that the set-up is in fact unique, in the sense that any vector-tensor Lagrangian must reduce to (2.1) upon linearization of the vector field about the vacuum A µ = 0, assuming the theory describes 3 + 2 dynamical degrees of freedom.
Our principal result is the outcome of the numerical calculation of the fundamental QNM frequency for the lowest multipole modes of the vector field, i.e. monopole (Fig.  1), axial dipole (Fig. 2), and polar scalar and vector dipoles (Figs. 3, 4). We explored a physically motivated range of values for the Proca mass µ, as well as the full range for the (normalized) non-minimal coupling parameter g 6 allowed by stability. However, our results exclude values very close to the bounds, eq. (2.4), where our numerics become unreliable. 13 From the truncated Taylor series for kM in eq. (5.5) one may also construct a Padé approximant in order to get a better fit of the numerical results. We thank Hector Silva for pointing this out to us.
It would be desirable to gain a better grasp on the behavior of QNMs when g 6 is at or arbitrarily close the critical values. This is not merely an academic question, since we recall that g 6 ≡ G 6 /r 2 g depends on the BH mass, so for any given non-zero coupling G 6 there will be a BH mass value such that either of the bounds is saturated. Of course, whether such a BH mass is physical, and whether the theory at that scale still makes sense, is a different question.
In the case where the vector field is massless the set-up simplifies considerably thanks to the U (1) gauge symmetry of the theory, and one is left with a single mode (for each , m) in each of the polar and axial sectors, i.e. the analogs of the electric and magnetic fields, allowing us also to compute the first two overtone QNM frequencies (n = 1, 2) as functions of g 6 . One interesting, although perhaps not unexpected, conclusion is that the isospectrality between polar and axial QNMs is broken by the non-minimal coupling, cf. Another set of valuable observables in the gauge-invariant setting is given by the electromagnetic susceptibilities, corresponding to the linear response of the BH to a static external field. While for a minimally coupled U (1) field BHs in GR do not polarize, as recalled in the introduction, our results demonstrate that this property ceases to hold in the presence of the Horndeski non-minimal coupling that we studied. We have shown this here explicitly for the dipole modes, cf. Fig. 8, for which we also provided some analytical understanding of the dependence of the susceptibility coefficients at linear order in the parameter g 6 , cf. eq. (5.5). We plan to undertake a more general analysis in a dedicated work.
The question on the stability of astrophysically relevant GR backgrounds under fluctuations of generalized vector fields motivated us also to study quasi-bound state solutions within our set-up. Unlike QNMs in asymptotically flat spacetimes, quasi-bound state frequencies may in principle develop a positive imaginary part, signaling a tachyon-type instability. Whether this indeed can occur for vector fields is an important issue because a tachyonic destabilization is a possible mechanism to generate compact astrophysical objects with vector hair starting from a hairless initial state-a phenomenon known as vectorization [81][82][83][84]. The no-go result of [39], together with the more general analyses in [85,86], cast doubt on vectorization as a viable mechanism, as they showed that localized perturbations must be either stable or else grow through wrong-sign kinetic or gradient operators. Our present results further supplement this claim by demonstrating that global bound-state solutions on a Schwarzschild BH background likewise do not exhibit tachyonic growth. We warn the reader that our argument does have some potential loopholes, as we explained at length in Sec. 4.2, although they are not expected to be critical. As an incidental outcome of our analysis, we also showed how an integral formula for the imaginary part of the quasibound state frequency due to Horowitz and Hubeny may be applied to asymptotically flat spacetimes. We hope to revisit this problem in a more general setting, e.g. by including matter and spin, in future investigations.

Acknowledgments
We would like to thank Luca Santoni, Hector Silva and Shuang-Yong Zhou for useful conversations and comments. The work of SGS and JZ at Imperial College London was supported by the European Union's Horizon 2020 Research Council grant 724659 MassiveCosmo ERC-2016-COG. The work of AH at Imperial College London was supported by the Royal Society International Newton Fellowship NIF\R1\191008. AH also acknowledges that the work leading to this publication was supported by the PRIME programme of the German Academic Exchange Service (DAAD) with funds from the German Federal Ministry of Education and Research (BMBF). SGS thanks the Peng Huanwu Center for Fundamental Theory at USTC for generous hospitality. JZ is also supported by a scientific research starting grant No. 118900M061 from the University of Chinese Academy of Sciences.

A.1 Method
In this appendix we describe the numerical method used to compute QNMs and quasi-bound states in this paper.
In Sec. 2.2 we have decomposed the Proca equation into its angular and radial components. The QNMs and quasi-bound states can be found by solving the non-linear eigenvalue problem for the radial equations (2.17), (2.18), (2.19), (2.20) and (2.24), supplemented with the appropriate boundary conditions discussed in Sec. 2.5. For the mode equations to be amenable to our numerical routine, it is useful to factor out the wave behavior at the boundaries. This is achieved by redefining with the choice of + sign for QNMs and − sign for quasi-bound states, and where B(r) is a regular function of r (also implicitly of ω) that tends to constant values as r * → ±∞. In the case of axial perturbations, as well as monopole and massless polar perturbations, the radial equation can be written as a single second order differential equation for the function B(r). In order to compute the eigenfrequencies, we first approximate the differential equations with finite-dimensional matrix equations using a collocation method with Chebyshev interpolation. We firstly introduce Note that one can choose other mappings ξ(r), the only requirement being that the singularities introduced by the non-minimal coupling terms are far enough from the domain of ξ such that they do not dramatically affect the convergence. Now we expand B(ξ) in terms of a set of cardinal polynomials p k (ξ), where p k (ξ) is defined by p k (ξ n ) = δ nk , and ξ n are the Chebyshev nodes, ξ n ≡ cos π(2n + 1) 2N + 2 , with n = 0, 1, . . . , N . where M nk (ω) ≡ p k (ξ n ) + C 1 (ω, ξ n )p k (ξ n ) + C 2 (ω, ξ n )δ nk . (A.7) The derivative matrices p k (ξ n ) and p k (ξ n ) can be computed using the second barycentric form [87,88], explicitly p k (ζ n ) = w k /wn ζn−ζ k n = k − k =n p k (ζ n ) n = k , (A.8) p k (ζ n ) = 2p k (ξ n )p n (ξ n ) − 2p k (ξn) ξn−ξ k n = k − k =n p k (ζ n ) n = k . (A.9) With a good initial guess on the eigenfrequency, in our case e.g. the eigenfrequency of the standard Proca field [30], one can solve Eq. (A.6) for ω and the set B(ξ k ).
For the massive polar modes, we instead have two coupled differential equations, say for B 2 (r) and B 3 (r), after we make the ansatz (A.1) for u 2 and u 3 . The procedure is nevertheless the same, i.e. we approximate the two differential equations with a set of algebraic equations of the form (A.6), now with and with the matrix M nk (ω) being enlarged accordingly.

A.2 Accuracy checks
For all the presented figures we have performed accuracy checks (i) at the points closest to the respective poles at g 6 = −1/2 and g 6 = 1 as well as (ii) at other random points. For the fundamental modes, these converence tests agree very well with expectations from an analytical error estimate discussed, for instance, in [33,App. C.4]. In particular, we find exponential convergence with growing number N of Chebyshev nodes. Moreover, we can also see that the convergence properties worsen with closeness to singularities in the complex plane.
We also find that higher modes (n > 1) beyond the fundamental (n = 1) are increasingly difficult to obtain because convergence worsens significantly. While we do not provide an analytical argument, we expect that the underlying reason is an increase in the number of oscillations (in the radial coordinate r) with growing n. The more oscillatory the behaviour, the harder it becomes to resolve these oscillations with fixed number of nodes N .

A.3 Values of QNM frequencies
We provide the numerical values of the frequencies in Tables. 2-5 for readers to make comparison.