Numerical simulations of single and binary black holes in scalar-tensor theories:circumventing the no-hair theorem

# Numerical simulations of single and binary black holes in scalar-tensor theories: circumventing the no-hair theorem

## Abstract

Scalar-tensor theories are a compelling alternative to general relativity and one of the most accepted extensions of Einstein’s theory. Black holes in these theories have no hair, but could grow “wigs” supported by time-dependent boundary conditions or spatial gradients. Time-dependent or spatially varying fields lead in general to nontrivial black hole dynamics, with potentially interesting experimental consequences. We carry out a numerical investigation of the dynamics of single and binary black holes in the presence of scalar fields. In particular we study gravitational and scalar radiation from black-hole binaries in a constant scalar-field gradient, and we compare our numerical findings to analytical models. In the single black hole case we find that, after a short transient, the scalar field relaxes to static configurations, in agreement with perturbative calculations. Furthermore we predict analytically (and verify numerically) that accelerated black holes in a scalar-field gradient emit scalar radiation. For a quasicircular black-hole binary, our analytical and numerical calculations show that the dominant component of the scalar radiation is emitted at twice the binary’s orbital frequency.

## I Introduction

Scalar fields are ubiquitous in physics, either as a proxy for more complex interactions or as fundamental quantities in their own right. For example, one of the best studied modifications to general relativity is scalar-tensor gravity, in which space-time curvature couples to scalar fields (which are sufficiently light to be relevant for astrophysics and/or cosmology). The historical development of this theory goes back to the 1940s, and involves several research groups with different views on the physical interpretation of the scalar degree of freedom Goenner (2012); Brans (2005); Singh and Rai (1983). In recent times, interest in scalar-tensor gravity has been driven by theoretical attempts to unify gravity with quantum mechanics at high energies and solve the cosmological constant and hierarchy problems, as well as observations of the cosmic microwave background and the highly-anticipated direct detection of gravitational waves Fujii and Maeda (2003); Faraoni (2004); Berman (2007).

The simplest version of scalar-tensor gravity is Brans-Dicke theory Brans and Dicke (1961); Dicke (1962); Brans (1962), which consists of a single massless scalar , whose coupling to curvature is controlled by a dimensionless parameter . Generalizations include varying scalar-curvature couplings and a scalar potential (“Bergmann-Wagoner” theories Bergmann (1968); Wagoner (1970)) as well as the possibility of multiple interacting scalar fields: see e.g. Fujii and Maeda (2003); Damour and Esposito-Farese (1992); Chiba et al. (1997) for comprehensive treatments of the subject. These generalizations (commonly referred to as “scalar-tensor theories”) are compelling due to their simplicity, but (perhaps as a consequence) they are also very well constrained observationally. The classic book by Will Will (1993) presents an overview of the subject, and comprehensive reviews of the state of the art in experimental tests of gravitational theories can be found in Will (2005); Clifton et al. (2012).

In this paper we shall focus, for simplicity, on scalar-tensor theories involving a single scalar field. In general relativity, because of the conservation of total momentum for isolated systems, gravitational radiation is quadrupolar in nature. A possible smoking gun of scalar-tensor gravity is the existence of dipole radiation, essentially due to violations of the equivalence principle. Solar System experiments and observations of binary pulsar systems place strong constraints on the coupling functions of the theory Will and Zaglauer (1989); Damour and Esposito-Farese (1992); Will (1993); Damour and Esposito-Farese (1993, 1996a, 1996b, 1998); Will (2005); Damour (2007); Bhat et al. (2008); Lazaridis et al. (2009); Esposito-Farese (2011); Alsing et al. (2012); Freire et al. (2012), and other interesting constraints may come from the direct observation of gravitational radiation from binary systems in the near future Will (1994); Damour and Esposito-Farese (1998); Scharre and Will (2002); Will and Yunes (2004); Berti et al. (2005a, b); Yagi and Tanaka (2010); Yunes et al. (2012); Berti et al. (2012). Despite all of these observational constraints, striking and potentially observable astrophysical phenomena are still possible in these theories. Such phenomena include superradiant instabilities (see e.g. Arvanitaki and Dubovsky (2011); Dubovsky and Gorbenko (2011); Yoshino and Kodama (2012); Witek et al. (2012) for discussions of this phenomenon in the context of the “string axiverse” scenario Arvanitaki et al. (2010)) and the related possibility of floating orbits around rotating black holes (BHs) Cardoso et al. (2011).

### i.1 Classical no-hair theorems

Theoretical studies impose remarkable constraints and limitations on scalar-tensor theories. First of all, the famous BH no-scalar-hair theorems first proved by Hawking Hawking (1972), Thorne and Dykla Thorne and Dykla (1971), and Chase Chase (1970) state that stationary BH solutions in Brans-Dicke theory are the same as those in general relativity. These results have been generalized and expanded upon by many authors. For example, an extension to multiple scalars has been established by Heusler Heusler (1995), and an extension to Bergmann-Wagoner and theories has been established by Sotiriou and Faraoni Sotiriou and Faraoni (2012). The no-scalar-hair theorems have also been confirmed by numerical studies of gravitational collapse Scheel et al. (1995a, b); Shibata et al. (1994); Harada et al. (1997); Novak (1998); Kerimo and Kalligas (1998). More generally, it has been observed that the Kerr metric is a solution in a wide class of gravity theories Psaltis et al. (2008). However, the fact that stationary vacuum solutions of scalar-tensor theories agree with those of general relativity does not mean that the dynamics of BHs in these theories must be the same Scheel et al. (1995b); Kerimo (2000); Barausse and Sotiriou (2008); Berti (2013).

For a comprehensive discussion and literature survey of no-hair theorems, the reader is referred to the reviews of Bekenstein Bekenstein (1996) and Chruściel, Costa, and Heusler Chrusciel et al. (2012). Although the literature is vast, there are two basic assumptions lying at the heart of most no-hair theorems. The first is that of stationarity, whose necessity has been demonstrated by Jacobson’s “Miracle Hair Growth Formula,” a perturbative construction of a hairy BH with time-dependent scalar boundary conditions Jacobson (1999).

The second assumption is the truncation of the scalar-tensor action to second order in the derivative expansion. At this level, the most general action [Eq. (2) in the single-scalar case] is very simple and contains only three terms, whereas scalar-tensor gravity at the four-derivative level is too complicated to be studied in complete generality, and thus, attention is often restricted to particular models.

One such model is quadratically modified gravity, whose action contains all possible terms quadratic in the Riemann tensor, coupled to a scalar. In this theory, BHs have been studied perturbatively, and solutions with scalar hair have been found Yagi et al. (2012a); Yunes and Pretorius (2009); Konno et al. (2009); Yunes and Stein (2011); Pani et al. (2011). Moreover, in the special case of a scalar coupled to a topological invariant – namely, Einstein-dilaton-Gauss-Bonnet (EDGB) or Dynamical Chern-Simons (DCS) gravity – a no-hair theorem for neutron stars has been established Yagi et al. (2012b). The conclusion here is that spherically-symmetric neutron stars have vanishing scalar monopole moment, but higher-order scalar multipoles need not vanish. However, the presence of derivatives higher than second order in the field equations severely complicates the implementation of numerical simulations.

Although four-derivative actions generically lead to field equations with more than two derivatives, there are some noteworthy exceptions. One is the Einstein-Skyrme system, a nonlinear sigma model with target-space , containing a term in the action quartic in scalar derivatives, and admitting linearly-stable BH solutions with scalar hair which have been described numerically in Droz et al. (1991); Heusler et al. (1991, 1992); Shiiki and Sawado (2005). Another model with second-order field equations is the galileon Nicolis et al. (2009), which is related to both higher-dimensional Lovelock gravity Van Acoleyen and Van Doorsselaere (2011) and massive gravity Luty et al. (2003); de Rham et al. (2011). It satisfies Solar System constraints by means of the Vainshtein mechanism Vainshtein (1972), and a no-galileon-hair theorem for spherically-symmetric BHs has been recently established Hui and Nicolis (2012). Another interesting example is Bergmann-Wagoner scalar-tensor gravity coupled to non-linear electrodynamics, where the non-vanishing trace of the electromagnetic stress-energy tensor enters as a source into the scalar-field equation. This allows for stable asymptotically-flat BHs with scalar hair, which have been studied numerically in Stefanov et al. (2007a, 2008, b, 2011, 2009); Doneva et al. (2010).

### i.2 A generalized no-hair theorem

The “classical” no-hair theorems described in the preceeding paragraphs, which are statements about stationary vacuum space-times, have been extended to the context of a BH binary system. Employing the “generalized EIH” formalism developed by Eardley Eardley (1975), Will and Zaglauer Will and Zaglauer (1989) have shown that the leading-order post-Newtonian (PN) dynamics of a BH binary in Brans-Dicke theory is indistinguishable from that in general relativity. Recently, this result has been extended to general scalar-tensor theories in the extreme mass-ratio limit Yunes et al. (2012) and it has been shown to hold up to 2.5PN order for generic mass ratio Mirshekari and Will (2013). Thus, even a dynamical (vacuum) spacetime with two interacting BHs does not have scalar hair in the PN limit. We can regard this conclusion as a “generalized no-hair theorem.”

The generalized no-hair theorem relies on the assumption that the binary system is isolated, in the sense that cosmological and environmental effects (say, due to the galactic background surrounding the BHs) are neglected. More precisely, it is assumed that: (1) there is no matter in the system, (2) the scalar field has zero potential, (3) the scalar-tensor action is truncated to second order in the derivative expansion, and (4) the metric is asymptotically flat (in all conformal frames) and the scalar is asymptotically constant.

The vacuum assumption (1) can be relaxed either by considering BHs in astrophysical environments or by considering compact stars, which are affected by the well-known spontaneous scalarization phenomenon (cf. Barausse et al. (2012) for recent numerical studies). Other recent numerical studies created a scalar-field “bubble” around the binary by using a nonvanishing potential, i.e. relaxing assumption (2). They found that the scalar-field bubble is rapidly accreted by the BHs, modifying the binary dynamics Healy et al. (2011). As for assumption (3), compact binary dynamics in quadratically modified gravity has been studied analytically in a perturbative framework, where one takes the point of view that the model should be considered as an effective low-energy theory Yagi et al. (2012b, c, 2013). Whether these theories are well posed for numerical evolutions is currently a matter of debate. We will not consider this problem in the present paper, but it is an interesting topic for future research.

Relaxing either assumption (3) or assumption (4) introduces a new length (or time) scale in the BH binary dynamics. In the case of assumption (3) this scale is determined by the Compton wavelength of a heavy particle, whose square enters into coefficients of four-derivative terms in the action. In the case of assumption (4), the new scale is determined by cosmological and/or galactic effects. A priori, it is not obvious which of these effects is dominant, and thus it is worthwhile to explore both possibilities. To our knowledge, the relaxation of assumption (4) has not been explored in the literature, and it is the main focus of our paper.

Asymptotic flatness of the metric is only an approximation to the dynamics of an astrophysical binary. Observations show that the Universe is expanding on timescales which are very large, but nevertheless finite with respect to astrophysical BH binary evolution. As shown in Jacobson (1999); Horbatsch and Burgess (2011), imposing time-varying boundary conditions endows the BHs in a binary with scalar charge, and therefore the binary can emit dipole scalar radiation. Furthermore, many cosmological models consider the existence of background scalars which can be anchored on matter Sin (1994); Hu et al. (2000); Matos and Guzman (2000); Schunck (1998); Sahni and Wang (2000); Matos et al. (2000); Arbey et al. (2001, 2003); Boehmer and Harko (2007); Alcubierre et al. (2002). In this case, one can for instance conceive of a BH binary evolving in the background of a nearly-static but nonuniform scalar field anchored on the galactic matter. The characteristic lengthscale of such a scalar-field profile would be much larger than the binary separation, and therefore it would have the same effect on the dynamical evolution of the binary system as the enforcement of boundary conditions which are not asymptotically flat. Finally, BH dynamics in the background of scalar fields could also be relevant to understanding accretion inside hypothetical supermassive boson stars, where huge scalar-field gradients are expected Macedo et al. (2013).

Scalar-field gradients can therefore allow us to circumvent the generalized no-hair theorem, i.e., to have a spacetime which contains only BHs and still emits scalar radiation. Indeed, as we shall discuss below [cf. Eq. (28), Section III.2.2], in the presence of a spatially-varying scalar-field profile , a nonrotating BH of mass with world-line would have a scalar charge1

 Q(t) = 4M2d→x(t)dt⋅→∇φ(→x(t)) = 8πσM2d→x(t)dt⋅^z,

where in the second line we have assumed that the scalar-field gradient is directed along the -axis, and we have parametrized its magnitude by a real parameter . If the BH is accelerated, or if the scalar gradient is nonuniform, the scalar charge would evolve in time, yielding scalar radiation. As we show in Appendix D, for a stellar-mass BH () moving near the galactic center a typical scalar-field gradient is . For a supermassive BH with a typical gradient could be as large as , comparable in order of magnitude to the numerical simulations presented in this paper.

### i.3 Executive summary and plan of the paper

The main goal of this work is to explore the consequences of the presence of a scalar-field gradient, which is equivalent to imposing nontrivial boundary conditions on the dynamics of a BH binary, and to verify numerically whether, as suggested by Eq. (I.2), such a setup can allow scalar radiation from a BH binary system in scalar-tensor theory. Here we present an executive summary of our main results and an outline of the paper.

In Section II we lay out our theoretical framework by introducing generic scalar-tensor theories and presenting the relations that allow us to transform between the Jordan frame (where physical quantities should be computed) and the Einstein frame (where we will perform our calculations). In particular, we show how gravitational radiation in the Jordan frame can be computed from a knowledge of the Newman-Penrose scalars in the Einstein frame.

In Section III we introduce analytical approximations for scalar fields in the background of single and binary BH spacetimes. These approximations are useful to validate (and provide insight into) our numerical simulations. In fact, numerical evolutions of initial data corresponding to a single black hole in a scalar gradient show that, after a short transient, the scalar field relaxes to the static configurations predicted by these perturbative calculations. For a quasicircular black-hole binary, in Section III we show analytically that the dipole component of the scalar radiation is emitted at twice the binary’s orbital frequency. This prediction is validated by our numerical simulations, which also show that the dipole component dominates the scalar emission.

In Section IV we present the details of our numerical implementation. The results of our simulations are discussed and compared with analytical results in Section V. In Section VI we summarize our findings and point out possible directions for future work.

To improve readability, in the Appendices we collect technical material that illustrates various important points of our analysis. Appendix A shows that a BH moving with constant velocity in a uniform scalar-field gradient does not emit scalar radiation. Appendix B (which is complementary to Section III.2.2) collects some lengthy formulas illustrating the structure of gravitational radiation from a BH binary moving in a scalar-field gradient. In Appendix C we provide explicit expressions for the evolution equations used in our numerical code. In Appendix D we estimate the order of magnitude of the scalar-field gradients expected in scalar-field models of dark matter.

## Ii Theoretical Framework

We focus on general single-scalar-tensor theories in vacuum with vanishing scalar potential, and at most two derivatives in the action. These theories are equivalent to Einstein’s theory extended to include a minimally coupled scalar field with vanishing potential. This statement (which will be clarified below) has a nontrivial consequence: the addition of minimally coupled scalars to Einstein’s gravity allows one to study a multitude of scalar-tensor theories at once. For this reason our simple framework offers an opportunity to take a glimpse at a rather broad spectrum of physics beyond Einstein’s theory.

Our starting point is the action of a general scalar-tensor theory for a single scalar field , written as

 S=∫d4x√−g16πG(F(ϕ)R−8πGZ(ϕ)gμν∂μϕ∂νϕ−U(ϕ)), (2)

where is the Ricci scalar associated to the metric , and and are arbitrary functions (see e.g. Fujii and Maeda (2003) and references therein). This form of the action corresponds to the choice of the so-called “Jordan frame”, where the scalar field is nonminimally coupled with gravity and all other matter fields obey the equivalence principle. The dynamics of matter fields would be described by an additional term on the right-hand side, that we set equal to zero because we are interested in BHs in vacuum.

Our numerical evolutions are more easily performed in the so-called Einstein-frame representation, that is related to the Jordan-frame representation by a conformal rescaling of the metric. In the Einstein frame, the scalar field is minimally coupled with the metric tensor and it affects the scalar-matter coupling in the matter action . Working in the Einstein frame is convenient because we focus on pure BH spacetimes, i.e., we set . Since we are interested in the effect of boundary conditions, we shall assume for simplicity (as in Horbatsch and Burgess (2011), and at variance with Healy et al. (2011)) that the effect of the scalar-field potential is negligible: .

### ii.1 From Jordan to Einstein and back

With , the explicit transformations that recast the previous action in the Einstein frame are Damour and Esposito-Farese (1996b)

 gEμν = F(ϕ)gμν, (3) φ(ϕ) = ∫dϕ[32F′(ϕ)2F(ϕ)2+8πGZ(ϕ)F(ϕ)]1/2, (4) A(φ) = F−1/2(ϕ). (5)

The Einstein-frame action is then

 S=116πG∫[RE−gμνE∂μφ∂νφ]√−gEd4x, (6)

and it leads to the following equations of motion:

 GEμν = ∂μφ∂νφ−12gEμνgρσE∂ρφ∂σφ, (7) □Eφ = 0, (8)

where the label denotes quantities built out of the Einstein-frame metric . This will be the starting point of our analysis. It is important to stress that, even though we are formally investigating a minimally coupled theory, our results are in principle applicable to a wide range of scalar-tensor theories. By working in the Einstein frame we can focus on quantities that depend on the intrinsic properties of the binary, rather than quantities that would be measured by gravitational-wave detectors. The latter may be obtained for any specific theory using the transformation from the Einstein frame to the Jordan frame.

### ii.2 Gravitational waves in the Einstein and Jordan frames

In the Jordan frame, we describe perturbations of the metric and of the scalar field as follows:

 gμν=g(0)μν+hμν,ϕ=ϕ(0)+ϕ(1). (9)

In the Einstein frame this corresponds to , , and from Eqs. (3)-(5) one gets

 hμν=F(ϕ(0))−1(hEμν−g(0)μνF′(ϕ(0))ϕ(1)), (10) ϕ(1)=[32F′(ϕ(0))2F(ϕ(0))2+8πGZ(ϕ(0))F(ϕ(0))]−1/2φ(1). (11)

In general, gravitational waves in scalar-tensor theories have 3 degrees of freedom Eardley et al. (1973a, b). In the Einstein frame they correspond to the two transverse-traceless components of the metric perturbation, plus the scalar field. The calculation of these quantities does not present difficulties. The physical degrees of freedom, which should be computed in the Jordan frame, can be read off from Eqs. (10)-(11). An alternative procedure is presented in Ref. Barausse et al. (2012), which shows the transformation of the corresponding Newman-Penrose quantities, with similar results. In this work we will present results for the scalar field and for the curvature scalar in the Einstein frame, which is directly related to the two Einstein-frame polarization states and via

 ΨE4=¨hE+−i¨hE×, (12)

where dots denote time derivatives. Quantities in the Jordan frame can be found using Eqs. (10)-(11), once one specifies the underlying theory. From here onwards, having established the relation between the Einstein and Jordan frames, we shall work exclusively in the Einstein frame, dropping the label “E” from all quantities. Unless specified otherwise, we will use geometrical units and set . Note however that here is a bare gravitational constant, and it is different from the quantity measured by a Cavendish experiment.

## Iii Scalar fields in single and binary black-hole backgrounds: analytical approximations

In this Section we introduce approximate analytical solutions that describe single and binary BHs in a scalar field gradient. These solutions will be useful below, either as code checks or for the interpretation of our numerical results.

### iii.1 Single black holes: linearized analytical solutions

Let us assume that the scalar-field gradient is of such low amplitude that the scalar can be treated as a perturbative effect on the spacetime metric. Under this assumption we can neglect terms quadratic in the scalar field, and therefore the field equations reduce to

 Rμν =0, (13) □φ =0. (14)

We will drop this perturbative approximation in Section IV, where the field equations (7) and (8) will be solved numerically.

Equation (13) is of course identical to Einstein’s equations in vacuum. We will consider the Schwarzschild (and later the Kerr) metrics as background BH solutions, and we will solve the Klein-Gordon equation (14) on these backgrounds. For the reasons explained in the introduction we are interested in numerical evolutions in a scalar-field gradient. Therefore we will consider background scalar-field solutions generated by distant, fixed, infinite homogeneous planes with constant surface scalar-charge density . To our knowledge, Press Press (1972) was the first to study a closely related problem (his setup differs from ours in that he considered a spherical shell of scalar charge as the source of the scalar field). Here we recast some of his results in a form suitable for comparison with our numerical setup.

#### Spherically symmetric black-hole background

Let us first consider the Schwarzschild solution in isotropic coordinates:

 ds2 = −(1−M2^r)2(1+M2^r)2dt2+(1+M2^r)4(dx2+dy2+dz2) (15) = −(1−M2^r)2(1+M2^r)2dt2+(1+M2^r)4 ×[d^r2+^r2(dθ2+sin2θdϕ2))].

Here is the isotropic radius, which is related to the areal radius by .

The scalar-field equation (14) on this background admits a simple axisymmetric solution sourced by distant planes of constant scalar-charge density, i.e.:

 φext = 2πσ(r−M)cosθ=2πσ(^r+M24^r)cosθ (16) = 2πσz(1+M24^r2),

where is the direction orthogonal to the charged plane. If no BH is present, Eq. (16) reduces to the field generated by homogeneously charged infinite planes, with constant gradient . For large the constant-gradient behavior applies also to the case where the background is a Schwarzschild BH.

We shall take as initial condition and test our numerical framework by checking that, after a transient, the scalar-field profile settles to the analytical solution , up to corrections of second order in the scalar field.

In Appendix A we derive a solution describing a BH moving with small constant velocity in a direction orthogonal to the charged planes, and show that it does not emit scalar waves. Indeed, as we will show analytically in Section III.2 and numerically in Section IV, the BH must have nonvanishing acceleration in order to generate scalar radiation.

#### Rotating black-hole background

In the case of a rotating BH with a rotation axis that is not orthogonal to the charged planes, axisymmetry is lost; however, a simple solution can still be found Press (1972). Let us choose our coordinates so that the BH spins along the -axis, at an angle with respect to the direction orthogonal to the charge-carrying planes. Then it can be shown that a solution is

 φext = 2πσ(r−M)(cosγcosθ+sinγsinθcosϕfa) (17) = 2πσ(r−M)[zrcosγ+xrfasinγ],

where ,

 fa=1r−MR[δΓ[2−ia/δ]eaπ2δPia/δ1(r−Mδ)], (18)

denotes the -function and is an associated Legendre function of the first kind. At large distances , and it can be seen that the charged plane generates a uniform gradient in the -plane. Figure 1 shows contour plots of in the plane, for different values of and of the angle . Notice how the field lines are distorted and frame-dragged close to the horizon.

### iii.2 Black-hole binaries: analytical approximation for quasi-circular inspirals

Let us now turn to the more complicated case of a quasi-circular BH binary evolving in an external scalar field with constant gradient. We begin our discussion by addressing the delicate problem of specifying boundary conditions for this system, which will be used to perform the numerical simulations discussed in Sec. IV.

#### Boundary conditions

In our numerical setup the boundary conditions are imposed by fictitious faraway charges, which are not modeled in our numerical simulations. These charges are assumed to lie outside the numerical grid, and mimic an external profile due (say) to a galactic scalar-field background. This “large-scale field” acts as a sort of reservoir: when the local field begins falling into the BH, there is an ingoing flux which restores the scalar-field gradient. The presence of this field outside the boundaries of our numerical simulations is implemented through the boundary conditions. At large distances we want to allow for outgoing waves while imposing the existence of the scalar-field gradient, so we require the field to behave in the following way:

 φ=φext+Φ(t−r,θ,ϕ)r, (19)

where is the areal radius, the constant gradient corresponds to the external field , and the second term on the right-hand side is the solution of the homogeneous equation describing outgoing (approximately spherical) scalar waves. We thus get

 ∂∂r(rφ)+∂∂t(rφ)=∂∂r(rφext). (20)

Since the boundary conditions are defined at large distances, , and we can write

 ∂∂r(rφ)+∂∂t(rφ)=4πσrcosθ. (21)

#### Multipole expansion of the scalar field

The angular dependence of the scalar field can be described through a multipole expansion of the form

 Φ(t−r,θ,ϕ)=M+ni˙Di+12ninj¨Qij+⋯,

where is the function appearing on the right-hand side of Eq. (19), dots denote derivatives with respect to the retarded null coordinate and is the radial unit vector, which depends only on the angles . The calligraphic symbols , and denote the monopolar, dipolar and quadrupolar components of , respectively.

The full relativistic scalar equation to be solved is . In order to obtain a solution which describes the physics that we are interested in – namely, a BH binary in a scalar gradient – it is essential to impose correct boundary conditions, both at null infinity, and in the vicinity of the worldlines of the singularities of the two BHs. The boundary conditions at null infinity have been discussed above; the near-worldline boundary conditions that must be imposed in order to find an approximate solution are more complicated.

In the special case of a comparable-mass BH binary system with small size-to-separation ratio (or, equivalently, small orbital velocity) the problem of imposing correct near-worldline boundary conditions is substantially simplified when one employs a “point-particle” effective field theory, in which length scales smaller than the Schwarzschild radii are integrated out. In this effective field theory the matter action has the form

where is an index that runs over the bodies ( for a binary system), is the worldline of body , is the proper differential arclength along , and is the “effective point-particle Lagrangian” of body A. For a structureless particle of mass , one has .

The matter action (23) gives rise to sources in the field equations of the effective point-particle theory, and these sources automatically enforce the correct boundary conditions at the worldlines for both and . For instance, to leading nonrelativistic order, the scalar-field equation has the explicit form

 □fφ = 4πρφ (24) = 4π2∑A=1QAδ(3)(→x−→zA(t)),

where is the D’Alembert operator in flat space, is an explicit parametrization of the worldline , and are the scalar charges of the BHs.

Moreover, to leading nonrelativistic order, one finds that the multipole moments entering into Eq. (III.2.2) are given by

 M = ∫d3xρφ=Q1+Q2, (25) →D = ∫d3x→xρφ=Q1→z1+Q2→z2, (26)

and so on.

In general, calculating the scalar charges is a difficult problem. A simplification takes place if we assume that the BH masses have the same order of magnitude , and that

 σ≪√aM3∼1Mv, (27)

where is the typical orbital separation, and is the typical orbital velocity. Henceforth, terms of order will be dropped. Then the BH scalar charges may be found by Jacobson’s formula Jacobson (1999), which for Schwarzschild BHs yields

 QA(t) = 4M2A[∂φ(t,→zA(t))∂t+→vA(t)⋅→∇φ(t,→zA(t))] (28) = 8πσM2A→vA(t)⋅^z,

where is the velocity of body , and in the second line the full scalar field has been replaced by the zeroth-order field .

Let us specialize to quasi-circular orbits in the -plane, so that the trajectories take the simple form

 →zrel(t) = →z1(t)−→z2(t) (29) = a(t)[^ycosχ(t)+^zsinχ(t)], M→zCM(t) = M1→z1(t)+M2→z2(t), (30)

where is the orbital radius, is the orbital phase, is the angular frequency of the orbit, and is the center of mass of the binary system. The quantities , , and are all small (of order ), and vanish in the absence of radiation reaction. An explicit expression for their leading-order time evolution may be found by solving the 2.5PN equations of motion (given in Section 9 of Blanchet (2006)), while dropping conservative corrections. Carrying out this calculation yields

 a(t) = a(0)(1−256t5τq)1/4≃a(0)(1−64t5τq), (31) ˙χ(t) = ˙χ(0)(1−256t5τq)−3/8≃˙χ(0)(1+96t5τq), (32)

where and are the (constant) radius and angular frequency of the zeroth-order orbit, respectively, and

 τq=[M1M2(M1+M2)[a(0)]4]−1 (33)

is the time scale over which the quadrupole tensor radiation shrinks the orbit.

With an explicit description of the orbit in hand, the scalar charges and multipole moments (25)-(26) may be calculated (see Appendix B for the explicit expressions). In this way, we find that:

• monopole radiation is emitted at the orbital frequency:

 M = 8πσM[M(M21+M22)(→vCM⋅^z) + M1M2(M1−M2)(˙asinχ+a˙χcosχ)],

and it vanishes in the equal-mass limit;

• dipole radiation is emitted at twice the orbital frequency, and more precisely

 ˙→D=˙→DCM+˙→Drel,DC+˙→Drel,osc, (35)

where the “CM” term is emitted at the orbital frequency, the “DC” component is nonoscillatory, and the component oscillates at twice the orbital frequency: cf. Eq.(100).

The physical problem addressed here differs from the situation investigated in Horbatsch and Burgess (2011). That study considered time-dependent scalar boundary conditions, rather than a gradient, and it found that monopole radiation is absent, while dipole radiation vanishes in the equal-mass limit. One may expect dipole scalar radiation to be emitted at the orbital frequency, rather than twice the orbital frequency. The reason why this expectation is erroneous in our case is that we have a background field with a gradient directed along the orbital plane, which combines with the oscillatory component sourced by the orbital motion.

A simple toy model can provide us with a complementary and perhaps more intuitive way to justify the expectation that dipole radiation must be emitted at twice the orbital frequency. Let us consider a rotating source with frequency on a scalar-field background [in our “rotated” polar coordinates, see Eqs. (74) below]. The source will produce a modulation in the background field of the form

 φ=φext[1+f(ϕ−Ωt)]. (36)

Expanding in circular harmonics, and

 φ=2πσrsinθsinϕ(1+∑mfmeim(ϕ−Ωt)), (37)

which implies that the multipolar components of the field will have the following dependence:

 φlm∼(e−i(m+1)Ωt+e−i(m−1)Ωt)+constant. (38)

Therefore the contribution should oscillate with frequency , the contribution with frequency , the contribution with frequencies and , and so on. As we will show below, this behavior is consistent with our numerical simulations.

## Iv Numerical implementation

Our numerical implementation of scalar-tensor theory closely parallels Salgado et al. (2008), but borrowing notation and conventions from Zilhao et al. (2010). The physical system studied in Zilhao et al. (2010) was quite different, since that paper considered higher-dimensional Einstein gravity in vacuum. However the equations can be cast as a system involving a scalar field coupled to gravity via dimensional reduction, so they are formally similar to the system considered here, as we show below.

### iv.1 3+1 decomposition

As a preliminary step for our numerical implementation, we perform a decomposition of the spacetime (see Alcubierre (2008) and references therein). Let us consider a slicing of the spacetime in a set of three-dimensional surfaces . Introducing the normal to the surface and the projector

 γμν=gμν+nμnν, (39)

we write the four-dimensional metric in the form (; ):

 ds2 =gμνdxμdxν =−α2dt2+γij(dxi+βidt)(dxj+βjdt), (40)

where are the lapse and the shift, respectively, and

 ∂t=αn+β. (41)

We shall denote by the covariant derivative on , i.e., the covariant derivative with respect to the three-dimensional metric .

The Lie derivative with respect to , then, is . Defining the extrinsic curvature

 Kij≡−12Lnγij, (42)

we have the following evolution equations for the (three-dimensional) metric:

 (∂t−Lβ)γij=−2αKij. (43)

Furthermore, we define the scalar curvature as

 Kφ=−12Lnφ (44)

so that

 (∂t−Lβ)φ=−2αKφ. (45)

Comparing with the definitions of the variables and in Salgado et al. (2008), we find that

 Qi =Diφ√8πG=γiμ∂μφ√8πG, Π =Lnφ√8πG. (46)

Therefore,

 Kφ=−12√8πGΠ. (47)

The constraint equations (2.15) and (2.16) of Salgado et al. (2008) read

 (3)R+K2−KijKij =Π2+Q2f=4K2φ+DiφDiφ, (48) DlKl i−DiK =−ΠQif=2KφDiφ. (49)

The evolution equation (2.17) of Salgado et al. (2008) can be written as

 (∂t−Lβ)Kij=(∂t−Lβ)γikKk j =−DiDjα+α((3)Rij+KKij−DiφDjφ−2KikKk j). (50)

Note that since , this expression coincides with Eq. (2.23) of Zilhao et al. (2010). Taking the trace of Eq. (50) we have

 ∂tK−βl∂lK+DiDiα−α((3)R+K2−DiφDiφ), (51)

and using Eq. (48) we find

 (∂t−Lβ)K=−DiDiα+αKijKij+4αK2φ, (52)

that should be compared to Eq. (2.19) of Salgado et al. (2008). The evolution equation (2.18) of Salgado et al. (2008) can be written as

 1α(∂t−Lβ)Kφ=−12√8πGLnΠ =KφK−12αDiφDiα−12DiDiφ. (53)

All terms of this expression appear, with the same coefficients, in Eq. (2.37) of Zilhao et al. (2010). The additional terms in that equation which are not present here are due to the more complicated dynamics of the scalar field arising from dimensional reduction.

It can be useful to write the equations also in terms of the stress-energy tensor, which enables us to compare the scalar-field terms with those in Alcubierre (2008). From (7) we have

 8πGTμν=∂μφ∂νφ−12gμν∂αφ∂αφ. (54)

Since , , , defining, as on page 87 of Alcubierre (2008)

 ρ =nμnνTμν, (55) jα =−γαμnνTμν, (56)

we get

 8πG ρ =nμnν∂μφ∂νφ−12gμν∂αφ∂αφ (57) =2K2φ+12DiφDiφ, (58) 8πG ji =−8πGγiμnνTμν (59) =−Diφnμ∂μφ=2KφDiφ, (60)

where we have used the fact that (see Alcubierre (2008)), thus . Therefore, Eqs. (2.4.6) and (2.4.9) coincide with our constraint equations (48) and (49). Furthermore, we can compute the quantity (cf. page 89 of Alcubierre (2008))

 Sμν=γμαγνβTαβ. (61)

We have

 8πGSμν=DμφDνφ−12γμνDiφDiφ+2γμνK2φ, (62)

and the trace of this equation (since ) yields

 8πGS=−12DiφDiφ+6K2φ. (63)

Then, , and

 4πG[(S−ρ)γij−2Sij]=−DiφDjφ, (64)

therefore Eq. (2.5.6) of Alcubierre (2008) coincides with our Eq. (50).

### iv.2 Baumgarte-Shapiro-Shibata-Nakamura formalism

Our evolution equations use the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism, in which the dynamical variables are , defined as follows:

 γij=χ−1~γij(with γij=χ~γij), χ=(detγij)−1/3, ~Aij=χ(Kij−13γijK), Γk=γijΓkij=χ~Γk+12~γkj∂iχ. (65)

Alternative notations replace our variable by a variable defined as . Here we will use as our dynamical quantity.

The Einstein equations in the BSSN formulation in the presence of a scalar field, which appears through the quantities , are implemented in an extended version of the Cactus Cac () based Lean code Sperhake (2007); Sperhake et al. (2010) in the form given in Eqs. (107)-(116) of Appendix C. Mesh refinement for our simulations is provided by Carpet Schnetter et al. (2004), horizon diagnostics by AHFinderDirect Thornburg (1996, 2004) and BH binary initial data satisfying Eq. (13) by a spectral solver Ansorg et al. (2004) provided through Cactus as the TwoPunctures thorn.

### iv.3 Initial data

Our initial data sets consist of either a single BH or a BH binary evolving in a background scalar field with a nonvanishing gradient. The background scalar field is generated by distant sources, which are kept fixed in our time evolutions. As a simple way to enforce this scenario, imagine that the background scalar field is generated by infinite homogeneous charged planes with surface density . As we saw in Section III, if a BH is present in the spacetime and the scalar field is small enough to be treated in the linear approximation the metric is unaffected, but the equilibrium solution for the scalar field changes. This setup is quite similar in spirit to the one adopted by Palenzuela and collaborators Palenzuela et al. (2010, 2009), the main difference being that they dealt with electromagnetic fields, and that their external magnetic field is generated by a current loop far away from the system.

Under these assumptions, the solution we found in Eq. (16) should be a good approximation to the initial data describing a single, nonrotating BH in a scalar-field gradient. Except for regions close to the BH singularity, the term in parentheses in Eq. (16) can furthermore be neglected. Therefore we initialize the scalar field using the simplified expression . Our numerical simulations of single BHs show indeed that the initial data (16) and its approximate version yield (after a brief transient, which we exclude from our analysis below) virtually identical evolutions of the scalar field and of the spacetime metric.

### iv.4 Upper bounds on the field gradient from the threshold of black-hole formation

Our setup consists of a constant scalar-field gradient at large distances. Because the energy density in any classical field theory is proportional to the square of the field gradient, one expects a roughly constant energy density . The total mass in a region of linear dimension then scales like , and therefore . Therefore we expect that the initial data will contain a horizon for . This condition imposes a nontrivial constraint on the size of our numerical grid. Here we provide a more formal argument supporting this conclusion.

We focus on conformally flat backgrounds with , . Then the momentum constraint is identically satisfied, while the Hamiltonian constraint yields

 0=△ψ+18ψηij∂iφ∂jφ, (66)

where is the Minkowski metric. A scalar field with constant gradient is such that , and therefore we get the equation

 △ψ=−π2σ22ψ. (67)

Imposing regularity at , the solution of this equation is

 ψ=Asin(πσr/√2)r. (68)

An apparent horizon exists for . This condition is equivalent to finding the roots of

 xcotx=1/2,withx≡πσr√2. (69)

The smallest root of this equation is at , i.e.

 σ≃1.16556√2πr≃0.52469r, (70)

in good agreement with the previous back-of-the-envelope estimate.

Our boundary conditions are enforced at distance and , respectively, for single and binary BH simulations. This provides a formal upper bound of on the magnitude of the gradients that we can simulate. However, as we will see below, even smaller gradients can generate exponentially growing instabilities (indicating collapse) in our numerical evolutions.

## V Numerical results

In this Section we first discuss our results for single BH evolutions, verifying that for small values of the gradient they are in good agreement with the analytical predictions of Section III.1. Then we study gravitational and scalar radiation from BH binaries in a scalar-field gradient.

### v.1 Single black-hole evolutions

In order to compare our numerical results with the analytical predictions of Section III, it is important to minimize coordinate effects. Let us denote by the radial coordinate used in our numerical simulation, which coincides with the isotropic coordinate at , i.e. . This need not be true at later times, since the gauge can dynamically change during the evolution. In order to monitor the scalar field as a function of time and check that it eventually settles to a configuration close to the analytical solution of Eq. (16), we need to perform integrations over spheres at given values of . Furthermore, we compute the areal radii at these locations by performing spherical integrations of the metric components, as discussed in Ref. Witek et al. (2010).

The spherical-harmonic expansion2 of the scalar field reads

 φ(t,r,θ,ϕ)=∑lmφlm(t,r)Ylm(θ,ϕ), (72)

where

 φlm(t,r)=∫dΩφ(t,r,θ,ϕ)Y∗lm(θ,ϕ). (73)

Since the binary moves in the -plane, we use noncanonical rotated coordinates

 x = rcosθ, (74) y = rsinθcosϕ, z = rsinθsinϕ.

For small scalar-field gradients, the back-reaction on the metric is very small and we expect to recover the stationary solution found in Section III.1, i.e., we expect that after an initial transient

 φ(t→∞,r,θ,ϕ)=φext(r,θ)=φ0+2πσ(r−M)cosθ. (75)

where is a constant. Because , the dominant nonvanishing component at late times (up to quadratic corrections in the field) should be given by

 φ10 =√4π32πσ(r−M). (76)

In this expression the areal radius is effectively time-dependent, and it is computed dynamically during the evolution as described above.

The numerical and analytical predictions (solid and dashed lines, respectively) are compared in Figs. 2 and 3. Let us focus first on Fig. 2, which refers to . In this case the numerically extracted dipole mode asymptotes quickly to the analytical prediction. The inset shows the percentage difference between the analytical and numerical value of the scalar field at late times as a function of the extraction radius: the agreement is remarkable at large radii but it gets progressively worse as we get closer to the BH, most likely because of gauge effects (we can exclude that the deviations are due to nonlinear effects, because these would scale with , whereas the disagreement seems to be independent of ). To summarize, the analytical and numerical predictions agree within a few percent for gradients .

The situation changes significantly for larger values of . In Fig. 3 we overplot the numerical and analytical components for and for a larger gradient, . As expected, when rescaled by their respective gradients the dipolar components are essentially the same, and they also converge to the value predicted by the analytical solution. However for this convergence can only be seen at intermediate times, whereas at late times the mode develops an instability: the figure shows that the roughly constant “baseline” value of the field seems to be superimposed to an exponentially growing oscillation that develops on timescales , both when we evaluate the field numerically (continuous black line) and when we use the areal radii to compute the analytical solution (dashed black line).

The existence of this instability is confirmed by Fig. 4, where we look at higher multipoles of the field on a logarithmic scale for (left panel) and for (right panel). The right panel of this plot shows that for the larger gradient an exponentially growing instability with growth time is present also in the subdominant , and multipoles. The left panel illustrates that when the instability (if it is present at all) develops on much longer timescales , so it does not affect our numerical simulations. Notice that early-intermediate time results are physically consistent: for both values of the gradient the scalar-field distribution is dominated by the dipolar component, and it is in good agreement with analytical predictions.

In summary, single-BH simulations in a scalar-field gradient show that our numerical evolutions are stable and reliable as long as the gradient is not too large. This conclusion is compatible with the arguments presented in Section IV.4 above.

### v.2 Binary black-hole evolutions: scalar and gravitational radiation

In this Section we discuss our numerical evolutions of BH binaries in a scalar-field gradient. In this initial study we focus on evolutions of nonspinning, unequal-mass binaries with a gradient and mass ratio , because unequal-mass binaries display two interesting features which would be absent by construction in the equal-mass case: center-of-mass recoil [see Eqs. (92)-(94)] and monopole scalar radiation [see Eq. (95)].

Gravitational waveforms, as characterized by the Newman-Penrose scalar , are shown in Fig. 5. For such low values of the scalar-field gradient, the impact on gravitational radiation emission is hardly noticeable.

The emission of scalar radiation is much more interesting. The scalar field acquires a nontrivial profile due to the dynamics of the orbiting BH binary. The scalar radiation has a crucial dependence on the binary setup, and more specifically on the angle between the orbital angular momentum of the binary and the direction of the scalar-field gradient. If this angle is zero, then effectively the individual BHs do not traverse any field gradient, and the scalar profile is expected to be trivial. Our numerical results confirm this expectation: the output of these simulations is indistinguishable from vacuum evolutions in pure general relativity Berti et al. (2007); Sperhake et al. (2010).

On the other hand, the induced scalar radiation should be maximized when the orbital angular momentum is perpendicular to the field gradient, so we now focus on this case. Our results are summarized in Figs. 5 and 6.

Because the binary evolves on the background of a dipolar scalar-field profile, this constant “background” value of the scalar shows up as a large imaginary component3 of the scalar-field modes, which is apparent in the left panel of Fig. 6 (in fact, we had to rescale the imaginary component by a factor in order to show this on this plot).

The real part of displays interesting dynamics (the imaginary component also has similar dynamics, but this is partially masked by the large background dipolar field, so the analysis of the real part turns out to be numerically “cleaner.”) At any extraction radius is initially zero, as the binary is simply traversing a constant scalar field. As the binary evolves, we expect to see scalar radiation crossing the extraction surface and producing a nonvanishing scalar profile. This is indeed observed in Fig. 6, where we show that has a “wavy” pattern at any fixed extraction radius. The scalar-wave nature of this pattern is well illustrated by the right panel of Fig. 6. There we take the time derivative of at the largest and smallest extraction radii, scaling the amplitude of the signal by the ratio of the extraction radii (as expected for a wave scaling like ), and we observe that the signal is indeed dipole scalar radiation emitted at twice the orbital frequency of the binary, consistent with the predictions of Sec. III.2.2. Furthermore, in Eq. (106) of Appendix B we show that the amplitude of the mode is consistent in order of magnitude with analytical predictions. We also observe a monopole component whose amplitude is comparable to the amplitude of , consistently with analytical predictions.

To understand the merger signal, when the two BHs collide and relax to a final nearly stationary state, it is useful to remember that, in vacuum, the merger of a BH binary with mass ratio produces a Kerr BH with spin Berti et al. (2007). Thus the lowest ringdown frequencies are expected to be, from perturbative calculations, for a scalar field, gravitational mode Berti et al. (2006, 2009). We find that indeed rings down with , in good agreement with perturbative calculations. An analysis of yields a ringdown frequency (with errors ), which is roughly consistent with perturbative calculations of scalar () perturbations of Kerr BHs Berti et al. (2009). Our simulations also show that the ringdown phase, and indeed the entire scalar signal, scales with . We conclude that our results indeed represent linear effects, as opposed to nonlinear mode couplings.

Although not completely obvious, there is a small DC component in Fig. 6 (right panel), which we estimate to be

 |˙φDC11||˙φpeak11|∼0.2, (77)

at early times, where is the absolute value of the waveform at a local peak (maximum or minimum). This numerical estimate can be compared to the analytical prediction, Eqs. (102)-(103). We start by estimating the angular frequency through the waveform frequency, and we find . Using Kepler’s law, one can estimate the orbital separation, and these two ingredients together with relations (32) allow us to estimate the relevant ratio of the time derivatives of expressions (102)-(103). We find a ratio which is smaller by almost one order of magnitude. This discrepancy can probably be explained by numerical uncertainties and strong-field nonlinear effects.

It is apparent from Fig. 6 (left panel) that the modes display a “drift”: after all the dynamics has died away, the field does not return to zero. Our data implies that at late times for