Numerical simulations of single and binary black holes in scalartensor theories: circumventing the nohair theorem
Abstract
Scalartensor 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 timedependent boundary conditions or spatial gradients. Timedependent 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 blackhole binaries in a constant scalarfield 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 scalarfield gradient emit scalar radiation. For a quasicircular blackhole 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 scalartensor gravity, in which spacetime 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 scalartensor 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 highlyanticipated direct detection of gravitational waves Fujii and Maeda (2003); Faraoni (2004); Berman (2007).
The simplest version of scalartensor gravity is BransDicke 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 scalarcurvature couplings and a scalar potential (“BergmannWagoner” theories Bergmann (1968); Wagoner (1970)) as well as the possibility of multiple interacting scalar fields: see e.g. Fujii and Maeda (2003); Damour and EspositoFarese (1992); Chiba et al. (1997) for comprehensive treatments of the subject. These generalizations (commonly referred to as “scalartensor 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 scalartensor 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 scalartensor 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 EspositoFarese (1992); Will (1993); Damour and EspositoFarese (1993, 1996a, 1996b, 1998); Will (2005); Damour (2007); Bhat et al. (2008); Lazaridis et al. (2009); EspositoFarese (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 EspositoFarese (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 nohair theorems
Theoretical studies impose remarkable constraints and limitations on scalartensor theories. First of all, the famous BH noscalarhair theorems first proved by Hawking Hawking (1972), Thorne and Dykla Thorne and Dykla (1971), and Chase Chase (1970) state that stationary BH solutions in BransDicke 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 BergmannWagoner and theories has been established by Sotiriou and Faraoni Sotiriou and Faraoni (2012). The noscalarhair 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 scalartensor 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 nohair 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 nohair 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 timedependent scalar boundary conditions Jacobson (1999).
The second assumption is the truncation of the scalartensor action to second order in the derivative expansion. At this level, the most general action [Eq. (2) in the singlescalar case] is very simple and contains only three terms, whereas scalartensor gravity at the fourderivative 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, EinsteindilatonGaussBonnet (EDGB) or Dynamical ChernSimons (DCS) gravity – a nohair theorem for neutron stars has been established Yagi et al. (2012b). The conclusion here is that sphericallysymmetric neutron stars have vanishing scalar monopole moment, but higherorder 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 fourderivative actions generically lead to field equations with more than two derivatives, there are some noteworthy exceptions. One is the EinsteinSkyrme system, a nonlinear sigma model with targetspace , containing a term in the action quartic in scalar derivatives, and admitting linearlystable 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 secondorder field equations is the galileon Nicolis et al. (2009), which is related to both higherdimensional 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 nogalileonhair theorem for sphericallysymmetric BHs has been recently established Hui and Nicolis (2012). Another interesting example is BergmannWagoner scalartensor gravity coupled to nonlinear electrodynamics, where the nonvanishing trace of the electromagnetic stressenergy tensor enters as a source into the scalarfield equation. This allows for stable asymptoticallyflat 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 nohair theorem
The “classical” nohair theorems described in the preceeding paragraphs, which are statements about stationary vacuum spacetimes, 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 leadingorder postNewtonian (PN) dynamics of a BH binary in BransDicke theory is indistinguishable from that in general relativity. Recently, this result has been extended to general scalartensor theories in the extreme massratio 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 nohair theorem.”
The generalized nohair 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 scalartensor 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 wellknown spontaneous scalarization phenomenon (cf. Barausse et al. (2012) for recent numerical studies). Other recent numerical studies created a scalarfield “bubble” around the binary by using a nonvanishing potential, i.e. relaxing assumption (2). They found that the scalarfield 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 lowenergy 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 fourderivative 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 timevarying 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 nearlystatic but nonuniform scalar field anchored on the galactic matter. The characteristic lengthscale of such a scalarfield 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 scalarfield gradients are expected Macedo et al. (2013).
Scalarfield gradients can therefore allow us to circumvent the
generalized nohair 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 spatiallyvarying
scalarfield profile , a nonrotating BH of mass
with worldline would have a scalar
charge
where in the second line we have assumed that the scalarfield 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 stellarmass BH () moving near the galactic center a typical scalarfield 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 scalarfield 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 scalartensor 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 scalartensor 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 NewmanPenrose 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 blackhole 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 scalarfield 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 scalarfield 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 scalarfield gradients expected in scalarfield models of dark matter.
Ii Theoretical Framework
We focus on general singlescalartensor 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 scalartensor 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 scalartensor theory for a single scalar field , written as
(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 socalled “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 righthand side, that we set equal to zero because we are interested in BHs in vacuum.
Our numerical evolutions are more easily performed in the socalled Einsteinframe representation, that is related to the Jordanframe 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 scalarmatter 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 scalarfield 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 EspositoFarese (1996b)
(3)  
(4)  
(5) 
The Einsteinframe action is then
(6) 
and it leads to the following equations of motion:
(7)  
(8) 
where the label denotes quantities built out of the Einsteinframe 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 scalartensor 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 gravitationalwave 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:
(9) 
In the Einstein frame this corresponds to , , and from Eqs. (3)(5) one gets
(10)  
(11) 
In general, gravitational waves in scalartensor theories have 3 degrees of freedom Eardley et al. (1973a, b). In the Einstein frame they correspond to the two transversetraceless 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 NewmanPenrose 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 Einsteinframe polarization states and via
(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 blackhole 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 scalarfield 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
(13)  
(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 KleinGordon equation (14) on these backgrounds. For the reasons explained in the introduction we are interested in numerical evolutions in a scalarfield gradient. Therefore we will consider background scalarfield solutions generated by distant, fixed, infinite homogeneous planes with constant surface scalarcharge 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 blackhole background
Let us first consider the Schwarzschild solution in isotropic coordinates:
(15)  
Here is the isotropic radius, which is related to the areal radius by .
The scalarfield equation (14) on this background admits a simple axisymmetric solution sourced by distant planes of constant scalarcharge density, i.e.:
(16)  
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 constantgradient 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 scalarfield 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 blackhole 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 chargecarrying planes. Then it can be shown that a solution is
(17)  
where ,
(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 framedragged close to the horizon.
iii.2 Blackhole binaries: analytical approximation for quasicircular inspirals
Let us now turn to the more complicated case of a quasicircular 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 scalarfield background. This “largescale field” acts as a sort of reservoir: when the local field begins falling into the BH, there is an ingoing flux which restores the scalarfield 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 scalarfield gradient, so we require the field to behave in the following way:
(19) 
where is the areal radius, the constant gradient corresponds to the external field , and the second term on the righthand side is the solution of the homogeneous equation describing outgoing (approximately spherical) scalar waves. We thus get
(20) 
Since the boundary conditions are defined at large distances, , and we can write
(21) 
Multipole expansion of the scalar field
The angular dependence of the scalar field can be described through a multipole expansion of the form
where is the function appearing on the righthand 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 nearworldline boundary conditions that must be imposed in order to find an approximate solution are more complicated.
In the special case of a comparablemass BH binary system with small sizetoseparation ratio (or, equivalently, small orbital velocity) the problem of imposing correct nearworldline boundary conditions is substantially simplified when one employs a “pointparticle” 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
(23) 
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 pointparticle 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 pointparticle theory, and these sources automatically enforce the correct boundary conditions at the worldlines for both and . For instance, to leading nonrelativistic order, the scalarfield equation has the explicit form
(24)  
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
(25)  
(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
(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
(28)  
where is the velocity of body , and in the second line the full scalar field has been replaced by the zerothorder field .
Let us specialize to quasicircular orbits in the plane, so that the trajectories take the simple form
(29)  
(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 leadingorder 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
(31)  
(32) 
where and are the (constant) radius and angular frequency of the zerothorder orbit, respectively, and
(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:
and it vanishes in the equalmass limit;

dipole radiation is emitted at twice the orbital frequency, and more precisely
(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 timedependent scalar boundary conditions, rather than a gradient, and it found that monopole radiation is absent, while dipole radiation vanishes in the equalmass 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 scalarfield background [in our “rotated” polar coordinates, see Eqs. (74) below]. The source will produce a modulation in the background field of the form
(36) 
Expanding in circular harmonics, and
(37) 
which implies that the multipolar components of the field will have the following dependence:
(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 scalartensor 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 higherdimensional 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 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 threedimensional surfaces . Introducing the normal to the surface and the projector
(39) 
we write the fourdimensional metric in the form (; ):
(40) 
where are the lapse and the shift, respectively, and
(41) 
We shall denote by the covariant derivative on , i.e., the covariant derivative with respect to the threedimensional metric .
The Lie derivative with respect to , then, is . Defining the extrinsic curvature
(42) 
we have the following evolution equations for the (threedimensional) metric:
(43) 
Furthermore, we define the scalar curvature as
(44) 
so that
(45) 
Comparing with the definitions of the variables and in Salgado et al. (2008), we find that
(46) 
Therefore,
(47) 
The constraint equations (2.15) and (2.16) of Salgado et al. (2008) read
(48)  
(49) 
The evolution equation (2.17) of Salgado et al. (2008) can be written as
(50) 
Note that since , this expression coincides with Eq. (2.23) of Zilhao et al. (2010). Taking the trace of Eq. (50) we have
(51) 
and using Eq. (48) we find
(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
(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 stressenergy tensor, which enables us to compare the scalarfield terms with those in Alcubierre (2008). From (7) we have
(54) 
Since , , , defining, as on page 87 of Alcubierre (2008)
(55)  
(56) 
we get
(57)  
(58)  
(59)  
(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))
(61) 
We have
(62) 
and the trace of this equation (since ) yields
(63) 
Then, , and
(64) 
therefore Eq. (2.5.6) of Alcubierre (2008) coincides with our Eq. (50).
iv.2 BaumgarteShapiroShibataNakamura formalism
Our evolution equations use the BaumgarteShapiroShibataNakamura (BSSN) formalism, in which the dynamical variables are , defined as follows:
(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 scalarfield 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 blackhole formation
Our setup consists of a constant scalarfield 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
(66) 
where is the Minkowski metric. A scalar field with constant gradient is such that , and therefore we get the equation
(67) 
Imposing regularity at , the solution of this equation is
(68) 
An apparent horizon exists for . This condition is equivalent to finding the roots of
(69) 
The smallest root of this equation is at , i.e.
(70) 
in good agreement with the previous backoftheenvelope 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 scalarfield gradient.
v.1 Single blackhole 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 sphericalharmonic expansion
(72) 
where
(73) 
Since the binary moves in the plane, we use noncanonical rotated coordinates
(74)  
For small scalarfield gradients, the backreaction 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
(75) 
where is a constant. Because , the dominant nonvanishing component at late times (up to quadratic corrections in the field) should be given by
(76) 
In this expression the areal radius is effectively timedependent, 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 earlyintermediate time results are physically consistent: for both values of the gradient the scalarfield distribution is dominated by the dipolar component, and it is in good agreement with analytical predictions.
In summary, singleBH simulations in a scalarfield 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 blackhole evolutions: scalar and gravitational radiation
In this Section we discuss our numerical evolutions of BH binaries in a scalarfield gradient. In this initial study we focus on evolutions of nonspinning, unequalmass binaries with a gradient and mass ratio , because unequalmass binaries display two interesting features which would be absent by construction in the equalmass case: centerofmass recoil [see Eqs. (92)(94)] and monopole scalar radiation [see Eq. (95)].
Gravitational waveforms, as characterized by the NewmanPenrose scalar , are shown in Fig. 5. For such low values of the scalarfield 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 scalarfield 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 scalarfield
profile, this constant “background” value of the scalar shows up as
a large imaginary component
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 scalarwave 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
(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 strongfield 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