Black hole binaries and light fields: Gravitational molecules

In a nonrelativistic setting, the Klein-Gordon equation around a single BH reduces to the Schrödinger equation Brito et al. (2015a); Arvanitaki and Dubovsky (2011); Baumann et al. (2019a). There is thus a quasibound state structure for scalar fields around a BH which is identical to the spectrum of the hydrogen atom. Let us now consider the nonrelativistic limit of a scalar field around a BHB. To lowest order in a post-Newtonian expansion, the geometry of a binary (including that of a BHB, if we are not interested in near-horizon phenomena) can be written in the form

where

is Newton’s potential. Here, $M_{i}(i=1,2)$ are the individual component masses, and $\vec{r}_{i}(t)$ are their position vector. There are higher-order terms which depend on the specifics of the system, and which become relevant for relativistic and strongly gravitating systems, but which do not affect the physics we are interested in.

Using standard nonrelativistic limit procedures, we define the complex field $\Psi(t,\vec{r})$ as

Note that the Klein-Gordon field $\phi$ is assumed to be real throughout this work. The field $\Psi$ is, in general, complex. Assuming that the binary components are widely separated and that the angular frequency is so small that time-dependent terms in the Newtonian potential can be neglected, the Klein-Gordon equation (1) reduces to the Schrödinger equation

where we neglect the subleading (for weakly gravitating, nonrelativistic systems) terms

We can also recover Eq. (5) via a Lagrangian approach Baumann et al. (2019b).

Equation (5) is written in the lab frame, $x^{\mu}=(t,r,\theta,\varphi)$. It will be useful to write it in the binary rest frame (corotating frame), $\bar{x}^{\mu}=(\bar{t},\bar{r},\bar{\theta},\bar{\varphi})$, which we can do with the usual coordinate transformation

where $\partial_{\bar{\varphi}}=-\bar{y}\partial_{\bar{x}}+\bar{x}\partial_{\bar{y}}$ spans the orbital plane, and $\Omega$ is the BHB orbital angular velocity. The frames are related through

Denoting $\bar{\Psi}(\bar{t},\bar{r},\bar{\theta},\bar{\varphi})\equiv\Psi(\bar{t},\bar{r},\bar{\theta},\bar{\varphi}+\Omega\bar{t})=\Psi(t,r,\theta,\varphi)$, and remembering that $\nabla^{2}=\bar{\nabla}^{2}$, we can then rewrite Eq. (5) in the corotating frame as

where

and the (time-independent) potential $V$ is given by

where $r_{1,2}=\sqrt{(\bar{x}\mp a)^{2}+\bar{y}^{2}+\bar{z}^{2}}$ is the distance to BH $i$, and $\bar{x}=\pm a$ are the positions of each BH.

Equation (8) can be treated perturbatively when $\Omega$ is small. Let us first consider the unperturbed system,

Since the potential $V$ is not time dependent, we can consider the energy eigenstate problem; writing

we then have

Introducing prolate spheroidal coordinates

where $2a$ is the separation between the two BHs on the $\bar{x}$ axis and $-1\leq\eta\leq 1$, $1\leq\xi$, and $0\leq\chi<2\pi$, Eq. (12) becomes separable. Using the ansatz

with $m_{\chi}=0,\pm 1,\pm 2,\ldots$, we then find

where we define

Here, $A$ is a separation constant. $\bar{E}$ and $A$ are labeled by three integers $m_{\xi}$, $m_{\eta}$, $m_{\chi}$, characterizing the properties of the solutions of the coupled system. We will focus on bound-state solutions, for which $\bar{E}<0$.

It is then easily seen that this gravitational problem is completely equivalent to the quantum-mechanical Schrödinger equation for the electronic energy of a one-electron heteronuclear diatomic molecule. In particular, if we identify $Z_{1}+Z_{2}=\mu\alpha$, $Z_{2}-Z_{1}=\mu\Delta\alpha$ then the equations above are mathematically equivalent to the quantum-mechanical problem (in atomic units), where the nuclei have atomic numbers $Z_{i}$ each and are separated by a fixed distance $D\equiv 2a$ Burrau (1927); Wilson (1928); Brown and Steiner (1966); Nickel (2011); Leaver (1986). In particular, this system also describes the ionized dihydrogen molecule Burrau (1927); Wilson (1928). We thus have a formal equivalence between two similar systems, that of a molecule governed by electromagnetism and a simple binary system in a Newtonian setting. We will see below that the inclusion of full general-relativistic effects alters this picture only slightly.

Equations (15) are of spheroidal type. The first is an “angular”-type scalar spheroidal equation Leaver (1986); Berti et al. (2006), and it is coupled to the second (radial) equation through the (unknown) energy $\bar{E}$ and separation constant $A$. We have solved this system to find the characteristic energies $\bar{E}$, with two different methods. We use direct integration of the ordinary differential equations, shooting to the energy; in addition, we used a high-accuracy continued fraction approach to solve the same problem Leaver (1986). Our numerical results for selected values of the separation are shown in Table 1.

At zero separation, we are effectively dealing with one single BH, for which the energy levels are known to high precision. In spherical polar coordinates $(\bar{r},\,\bar{\theta},\,\bar{\varphi})$ the eigenfunction is

with $\mathcal{Y}^{m}_{\ell}$ being the scalar spherical harmonics and $L_{n}^{2\ell+1}$ a generalized Laguerre polynomial (it’s a polynomial of order $n$ in its argument). Note that $M=M_{1}+M_{2}$ is the BHB mass. With these definitions and conventions, the energy eigenvalue is

up to ${\cal O}(\alpha^{3})$. Higher-order expansions in $\alpha$ can be calculated using well-known techniques and are shown in Ref. Baumann et al. (2019a), with a different state label. Notice that we follow the state labeling of Ref. Brito et al. (2015a). The scalar profile of mode $(n,\ell,m)$ decays spatially as $r^{\ell+n}e^{-\gamma/2}$ and the angular profile is dictated by the corresponding spherical harmonic. The spatial extent of the scalar configuration is defined by the exponential decay, and is of order ${\cal S}\sim 1/(M\mu^{2})$.

Note that the prolate coordinates [Eq. (13)] are also given by

Defining spherical coordinates such that

one finds, when $a\to 0$,

The relation between quantum numbers is then, in this limit Bates and Reid (1968)

For small separations $a\mu$, our results—shown in Table 1—for the fundamental ($n=0$) mode are compatible with (up to terms of order ${\cal O}(a^{3}\mu^{3})$)

This analytic expansion was derived by Bethe, and in subsequent analytical work for the dihydrogen molecule and generalizations thereof Bethe (1933); Brown and Steiner (1966); Nickel (2011); Bates and Reid (1968); Abramov and Slavyanov (2001); Damburg et al. (1984). Our own numerical results in the limit of small binary separation are in perfect agreement with such expression.

We note that for our numerical simulations shown in Sec. IV, we will use a frame rotated by 90 degrees, which therefore superposes different $m$ states. In addition, and more importantly, at finite separation $a$, we can no longer use Eq. (20): a given state in the $(m_{\xi}\,,m_{\eta}\,,m_{\chi})$ basis involves, generically, a superposition of all modes in a spherical harmonic decomposition, such as that of Sec. IV.

Particle analogies of wave equations often play an important role in several physical phenomena. In particular, a lot of BH physics described by wave equations have particle descriptions which helps us to understand the system from a different point of view. For example, BH quasinormal modes (QNMs) can be related to closed null orbits around the BH based on the WKB approximation, where the QNM frequencies can be estimated from the orbital periods Cardoso et al. (2009). Conversely, the existence of closed null orbits around relativistic objects may hint at the possibility of the existence of QNMs Bernard et al. (2019).

We will therefore study the particle analog of the system we have been considering of a massive scalar field around a BHB. As we will now see, this particle description can be derived from the WKB approximation of Schrödinger’s equation and can be applied to states around a general separable spacetime.

In quantum mechanics, it is well known that the first order of the WKB approximation corresponds to the Hamilton-Jacobi equation for a classical particle, and the energy spectrum can be related to the particle motion through the Bohr-Sommerfeld condition Weinberg (2015)¹¹1Strictly speaking, the Bohr-Sommerfeld condition is $\oint p_{k}dq_{k}=2\pi\left(n_{k}+\frac{1}{2}\right)$; however, since we are only interested in its large-$n_{k}$ limit, we have neglected the $\frac{1}{2}$ term in Eq. (22).

where $q_{k}$ and $p_{k}$ are the canonical coordinates and $n_{k}$ is an integer. The integral is performed over a closed classical orbit which is described by the corresponding Hamilton-Jacobi equation. Since this is derived from the WKB approximation, we can use this relation in the large-$n_{k}$ limit to estimate the frequency of the bound states of Sec. II.1.

Furthermore, from the classically allowed regions of the Hamilton-Jacobi equation, we can have an idea about the spatial profile of the wave function. Since, as shown in the previous subsection, our system can be described by Schrödinger’s equation, we can apply these arguments to get the relation between (truly) bound states and closed orbits of a massive particle.²²2These are actual bound states because we are using a Newtonian approximation, and horizons are not present.

Let us then consider the classical motion under Newton’s potential [Eq. (10)], which corresponds to the motion of a massive particle (of mass $\mu$) around a BHB, in the binary rest frame. The Lagrangian for the particle is

where $(\bar{x}(\bar{t}),\bar{y}(\bar{t}),\bar{z}(\bar{t}))$ are the particle’s coordinates. Using the coordinates $(\xi,\eta,\chi)$, the Hamilton-Jacobi equation reads

where $S(\bar{t},\xi,\eta,\chi)$ is Hamilton’s principal function. In these coordinates the Hamilton-Jacobi equation is separable; we thus write

and after substituting into Eq. (23), we obtain

where $C_{0}$ is a separation constant. Comparing with Eq. (15), we see that $2a^{2}\mu C_{0}$ corresponds to $A$. For classical motion, the right-hand side of Eq. (24) must be positive, and the parameter space where this happens corresponds to the classically allowed region. The topology of this region can change depending on the parameters; for simplicity, we will focus on equal-mass binaries ($\Delta\alpha=0$) and on bound states around the binary ($\bar{E}<0$) with $m_{\chi}=0$.

Considering first the $C_{0}>0$ case, we see from Eq. (24) that $\eta$ is unconstrained ($-1\leq\eta\leq 1$), and the allowed region for $\xi$ is determined from

When $|\bar{E}|-\frac{\alpha}{a}+C_{0}<0$, the allowed region for $\xi$ is $\xi\in[1,\xi_{+}]$, and when $|\bar{E}|-\frac{\alpha}{a}+C_{0}>0$, the allowed region is $\xi\in[\xi_{-},\xi_{+}]$. Here, $\xi_{\pm}$ are solutions of the equation $|\bar{E}|\xi^{2}-\frac{\alpha}{a}\xi+C_{0}=0$:

Therefore, when $C_{0}>0$, the particle motion is an orbit around the binary (see the first and second panels in Fig. 1).

We now focus on the $C_{0}<0$ case. The allowed region for $\eta$ is $\eta\in\left[-1,-\sqrt{\frac{|C_{0}|}{|\bar{E}|}}\right]\cup\left[\sqrt{\frac{|C_{0}|}{|\bar{E}|}},1\right]$ for $|C_{0}|<|\bar{E}|$. Existence of solutions to Eq. (25) then implies that $|\bar{E}|-\frac{\alpha}{a}-|C_{0}|<0$, and the corresponding allowed region for $\xi$ is then $[1,\xi_{+}]$ (see the third panel in Fig. 1). In this case, the motion is an orbit around each individual BH.

Since the WKB approximation is valid only in the high-frequency limit, we cannot easily apply this classical picture to the states from the previous subsection. We can, however, see that there are two distinct profiles for these states: one corresponding to configurations bound to each individual BH, and another corresponding to a configuration bound to the binary as a whole.

In order to compare the particle picture to the states of the previous subsection, let us discuss the spectrum of these states using the Bohr-Sommerfeld condition Eq. (22). From this condition we immediately see that $m_{\chi}$ must be the integer $n_{\chi}$ and we further get

where $n_{\xi}$ and $n_{\eta}$ are integers, and the integration limits $\eta_{\pm}$ and $\xi_{\pm}$ are the roots of the expressions inside each square root. Since Eqs. (27) cannot be integrated analytically let us consider the hydrogen atom limit ($a\to 0$) by fixing $\ell^{2}=2a^{2}\mu C_{0}$ and $r=a\xi$ in the equal-mass case ($\Delta\alpha=0$). We obtain the spectrum

These expressions are in good agreement with the spectrum of the hydrogen atom in the large-$n_{\eta}$ and large-$n_{\xi}$ limit, confirming the validity of this particle picture.

Having solved Eq. (8) to zeroth order in $\Omega$, let us now consider, perturbatively, the effect of rotation. The first-order correction in $\Omega$ to the unperturbed energy levels that we have just computed is given by the expectation value of the rotation operator $L_{\bar{z}}=-i\partial_{\bar{\varphi}}$ for the system in the unperturbed state. In the coordinates $(\xi,\eta,\chi)$ this operator takes the form

Since $\bar{z}\sim\cos\chi$ and $\bar{y}\sim\sin\chi$, we can easily see that the expectation value of this operator in an eigenstate [Eq. (14)] is zero, $\langle L_{\bar{z}}\rangle_{\bar{\psi}}=0$. This shows that rotation effects only manifest themselves at order $\Omega^{2}$ and will therefore be neglected.

For our numerical implementation we ignore the backreaction of the massive scalar field on the BHB spacetime and follow the approach described in Ref. Bernard et al. (2019). In this approach one builds an approximate BHB background metric using the construction outlined in detail in Mundim et al. Mundim et al. (2014) (see also Ref. Johnson-McDaniel et al. (2009) for the equivalent construction used in the context of generating BHB initial data). It is important to mention that for our present construction we turn off the emission of gravitational radiation and therefore always consider binaries with constant separation. We do not turn off scalar radiation; the Klein-Gordon equation is evolved in full generality.

Our approach then consists in using this approximate BHB metric construction and solving the Klein-Gordon equation (1) in this spacetime. We emphasize that this metric, even though it is time-dependent, is prescribed—that is, it is not time-evolved. Our task is then to numerically solve Eq. (1) on a time-dependent background. To do so we write the equation in a first-order form by introducing

where $N$ and $\beta^{i}$ are the lapse function and shift vector, respectively. The resulting system is numerically evolved with a method of lines approach.

For our numerical evolutions, we use the Einstein Toolkit infrastructure Löffler et al. (2012); Zilhão and Löffler (2013); Babiuc-Hamilton et al. (2019) with Carpet Schnetter et al. (2004); Carpet for mesh refinement capabilities and the multipatch infrastructure Llama Pollney et al. (2011). The scalar field equations are evolved in time by adapting the ScalarEvolve code available in Ref. Witek et al. (2020), which was first used and described in Ref. Cunha et al. (2017). Since the background metric uses harmonic coordinates, in our evolutions we further excise the BH interior with the procedure outlined in Ref. Assumpcao et al. (2018). The overall infrastructure and evolution are essentially the same those used and tested in Ref. Bernard et al. (2019). We employ fourth-order-accurate finite-differencing stencils to approximate spatial derivatives, but there are lower-order elements in the code—in particular the so-called prolongation operation is only second-order accurate in time. Our results are compatible with a convergence order between orders 2 and 3, which is consistent with the overall setup. See Appendix B for further details.

We will evolve two different types of scalar field initial data, which will be referred to as “nonspinning” and “spinning” initial data. The first of these consists of a momentarily static Gaussian profile given by

where $A$ and $\sigma$ denote the amplitude and width of the Gaussian pulse, respectively.

Secondly, we will evolve configurations which, unlike the previous construction, have angular momentum. These are the “spinning” initial data, for which

with

and the initial configuration for the $K_{\phi}$ field can be trivially obtained from Eq. (29). Here, $A_{l,m}$ is the amplitude of each $(l,m)$ mode, which can be freely specified, and $\sigma$ is a typical width of the clump. $(r,\theta,\varphi)$ are the standard spherical coordinates around the center of mass of the system. The time dependence $m\varphi+\omega t~{}(m=\pm 1)$ in each term introduces angular momentum on the $z$ axis. Initial data with negative $m$ are corotating with the BHB, while initial data with positive $m$ are antirotating.

To analyze our results, we decompose the evolved field at fixed radial distance into spherical harmonics $Y^{m}_{l}(\theta,\varphi)$ as follows

Note that, since $\phi$ is real, $\phi_{lm}(t)=(-)^{m}\phi_{l,-m}^{\ast}(t)$ where $\ast$ denotes complex conjugation. We will further Fourier-transform $\phi_{lm}(t)$ to check its frequency spectra and compare with the results obtained in Sec. II.1. In order to do so, however, we must note that the frequencies computed in Sec. II.1 were computed in the rest frame of the binary, whereas here we will be extracting the data in the lab frame. To compare the data we then need to change frames once again, which can be done as follows:

The Fourier spectra in the lab frame is computed through [$\phi=\phi(t,r_{\rm ext},\theta,\varphi)$]

Given that $Y^{m}_{l}(\theta,\varphi)=Ne^{im\varphi}P^{m}_{l}(\cos\theta)$, where $P^{m}_{l}$ are the Legendre polynomials and $N$ is the normalization constant, using Eq. (6), we can write

relating the frequencies computed in the lab frame to the ones computed in the comoving frame.

Since in the following section we will be focusing on the real part of the multipolar components let us also write

where $\mathcal{F}_{\Re}\equiv\mathcal{F}[\Re(\phi_{lm})](\omega)$ and Eq. (36) is used in the last step. To connect with our upcoming results we further need to take the real part of Eq. (37),

where we have used

In conclusion, in the lab frame we expect to see a superposition of $\omega\pm m\Omega$ frequencies for each $m$ mode of the corotating frame.

Based on the scales above, we expect that for large enough couplings $M\mu$ the scale ${\cal S}_{i}$ obeys ${\cal S}_{i}<D$, in which case the scalar cloud localizes around each BH but not around the binary. In other words, we expect that BHs sufficiently far apart can support clouds as if in isolation. To test this, we evolve momentarily static Gaussian initial data, Eq. (30), corresponding to configuration nonspin1 in Table 2. For these parameters, the typical size of the cloud around each BH is ${\cal S}_{1}={\cal S}_{2}\sim 8M$, which is smaller than the BH separation. We therefore expect that quasibound states form around each BH and our results confirm this, as shown in Fig. 2. As can be seen, localized structures are apparent with the dominant mode being an $l=1,m=0$ state around each individual BH. Our results show that these configurations remain localized to the binary with negligible variation in its shape and topology for thousands of dynamical timescales, with only a slight variation in amplitude, hence qualifying as true quasibound states.

We will now discuss molecular-like structures—i.e., scalar clouds around BHBs. We focus on BH separations $D=10M,\,60M$, and we fix the mass of the scalar field to $\mu M=0.2$. This corresponds to configurations nonspin2 and nonspin3, respectively. The scales are such that now the size of the quasibound states, if present, would encompass the binary when $D=10M$. We will see that even for $D=60M$ such global states exist.

The evolution of these configurations is shown in Fig. 3. Perhaps the most evident aspect of these simulations is that there is a persistent structure, a cloud or quasibound state of scalar field around the binary for relatively long time. In the context of Fig. 3, the evolution timescale is large enough that the binary performed over ten periods for configuration nonspin2 and two periods for configuration nonspin3. This timescale is orders of magnitude larger than the free fall time, and yet the scalar structure persists throughout the evolution.

The second noteworthy aspect is that the state keeps, roughly, the symmetry of the initial conditions. This is clearly seen in the energy density plots in Fig. 3. However, it is clear also that fine structure arises, clearly seen for larger binary separations, excited by the presence of the binary, an asymmetric perturber. In particular, we observe the excitation of the quadrupole mode by the BHB. This is clearly shown in Fig. 4, where we show the monopole $l=m=0$ and quadrupole $l=m=2$ components of the field at selected “extraction” $r_{\rm ex}/M=40,\,60,\,90$.

Figure 4 shows that this is indeed a quasibound state, decaying exponentially in time, albeit on long timescales. This is exactly as expected from an analysis of massive fields around nonspinning BHs Brito et al. (2015a). The scalar at late times behaves as an exponentially damped sinusoid, as expected for quasibound states. A fit to the late-time behavior (see Fig. 4) shows that the lifetime of such structures is $3\times 10^{3}M$ for $D=10M$, and $10^{4}M$ for $D=60M$. These scales should be compared with the analytical prediction, Eq. (4.13) in Ref. Wong (2020). That expression is applicable to $D=10M$ and yields a timescale $\sim 3.9\times 10^{3}M$, in excellent agreement with our numerical results.

In accordance with our analysis in Sec. II, the scalar field is oscillating with a frequency $\sim\mu$. To quantify the agreement with the nonrelativistic analysis, we computed the Fourier spectrum of the monopole and quadrupole modes at different radii, and averaged the result. We find clear peaks at different frequencies, the dominant ones are shown in Table 4 for $D=60M$ and compared against the nonrelativistic prediction of Sec. II.³³3One needs to pay attention when comparing against the results of Sec. II, since a rotation of axes is involved and, as we mentioned, at finite separations a spherical harmonic mode maps into a superposition of the $(m_{\xi},\,m_{\eta},\,m_{\chi})$ states used in Table 1. The agreement is of order 0.1% or better, providing strong evidence that bound, molecular-like gravitational states do form around BHBs. Notice also that we find two (or more) peaks for $l=m$ modes. For $D=60M$ we can read from the table that their separation in frequency $\Delta^{\omega}_{m}\equiv M(\omega_{m}-\omega_{-m})$ is $\Delta^{\omega}_{m=2}=0.0082$, in good agreement with the expected prediction of Sec. III.3, $\Delta^{\omega}_{m}=2mM\Omega$ ($M\Omega\sim 0.0022$ for this example).

The interpretation of these states as molecular-like is further supported by the spatial profile of the scalar and energy density, shown in Fig. 5. The time evolution of the energy density along the $x$ axis is depicted in the three panels of this figure. The late-time profile around the binary is well described by a density $\sim e^{-2r/(25M)}$. This profile is in accord with the nonrelativistic, single BH expression [Eq. (16)] for the quasibound state. Our results indicate that this is also a good expression, even for the moderately large separations that we studied. Figure 5 shows that the exponential falloff gives rise to a slowly decaying but small tail of energy density for distances $r\gtrsim 100M$. This looks to be a slow, radiative component part of the signal.

In fact, this gravitational system behaves like a rotating molecule in quantum mechanics.⁴⁴4Disclaimer: one should not read too much in this analogy. The shared electron cloud in a molecule is what binds the atoms together. Here, the “atoms” (two BHs) are bound together by gravity and then the “electron cloud” is a solution of the scalar wave equation on that background. The snapshots of the scalar field and energy density shown in Fig. 3 show that the scalar field is not only localized around the BHB, but that it is dragged by it, rotating counterclockwise, along the orbital motion of the binary. The field shows modulations at low-frequencies—in particular, at $2m\Omega$—as can be seen in the figures. Such modulation is the expected for an equal-mass binary. Notice that the signal is almost equally modulated in amplitude at different extraction radii; thus the low-frequency envelope is not the result of beatings caused by overtone excitation. In other words the scalar is indeed gravitationally dragged by the binary. The period of such a pattern is the same as the orbital period of the binary.

The previous evolutions referred to spherically symmetric, momentarily static initial data. We now consider time-asymmetric initial data. Let us focus on dipolar initial data corotating with the binary—i.e., configurations spin1 and spin2 of Table 3. The evolution of the spatial profile of the field and its energy density is shown in Figs. 6 and 7. A dipolar pattern stands out from these plots.

Figures 6 and 7 show that for compact binaries (when the BH separation is smaller than the typical size of the cloud), the energy density of the state acquires a torus-like shape, centered on the binary, and supported by its angular momentum. This is similar to the topology of quasibound states around single BHs. On the other hand, for large BH separations as in Fig. 7, the profile is no longer connected; the torus “breaks up”, and leaves two overdensity clumps of scalar field corotating with the BHB.

The presence of the binary excites other modes with similar symmetries to that in the initial data. In particular, we see a strong octupole $l=m=3$ mode, shown in Fig. 8 for $D=60M$. Notice that the octupolar mode grows from a negligible value to roughly 10% in amplitude of the dipolar component. Notice also the large timescales involved: the amplitude of these components is approximately constant up to timescales of order $\sim 10^{4}M$ or larger. The modulation in the signal is, as we explained before, due to the motion of the binary, and has a frequency $2m\Omega$ as expected for an equal-mass binary.

These results indicate that the evolution drove the system to a quasibound state, a relativistic analogue of the molecular solutions discussed in Sec. II.1 for the nonrelativistic system. Together with the previous results, and as we will insist below, these features indicate that the formation of quasibound states is a robust result for general initial conditions. The analysis of the Fourier-decomposed signal shows a frequency content which is peaked at $M\omega=0.1992,\,0.1948$. This is in good agreement with the nonrelativistic predictions of Sec. II.1 (which yield $0.1994$ and $0.1951$, respectively). One finds (see Table 4) $\Delta^{\omega}_{m=1}=0.0044$, also in very good agreement with the expectations.

The energy density of the configuration is shown in Fig. 9. These results show clearly that the initial conditions are similar to those of a quasibound state, and the system evolution does not take it away significantly from the initial conditions. The asymptotic behavior of the density profile is well described by $(re^{-r/50M})^{2}$. Note that the nonrelativistic prediction [Eq. (16)] for the scalar profile is of the form $\sim e^{-M\mu^{2}r/(\ell+1)}=e^{-r/(50M)}$ for the fundamental mode when $M\mu=0.2$. This is in perfect agreement with our time-evolution results around a BHB.

Similar results hold for the evolution of counterrotating initial data (with respect to the binary, thus data rotating clockwise). We have evolved configuration spin3 of Table 3 and the corresponding profiles of the field and energy density are shown in Fig. 10.

Since the initial profile has angular momentum which is opposite to that of the binary, the scalar field rotates in a direction opposite to the BHB. However, the energy density of the scalar cloud does rotate in the same direction as that of the binary. The asymptotic late-time spatial distribution of the energy density is well described by a radial dependence $(re^{-r/50M})^{2}$, consistent with a nonrelativistic analysis of quasibound states.

As we discuss in Appendix A, selection rules imply that other modes must be excited—in particular, the $l=m=3$ mode. Figure 11 shows precisely this. It is also clear that a very long-lived quasibound state forms (the lifetime is so large, in fact, that we were unable to estimate it).

We investigated other types of initial conditions, such as narrower pulses, with smaller width $\sigma$. This changes the amount of field that is dissipated in the early stages, but a quasibound state always ends up forming, with a phenomenology similar to what we described.

We also studied data with higher initial frequency content: since quasibound states have $\omega\sim\mu$, it is conceivable that high-frequency initial data just dissipates away. Our results show that this does not happen. Instead, again, the system evolves towards a frequency content $\omega\sim\mu$ and a spatial distribution described well by a nonrelativistic quasibound state. It seems that these quasibound states are an attractor in the phase space.

We should also add that quasibound states decay exponentially in space, far away from the binary. At large distances, the dominant behavior is controlled by power-law tails Hod and Piran (1998); Koyama and Tomimatsu (2001, 2002); Witek et al. (2013). Asymptotically, our results are consistent with a late-time decay $\phi_{lm}\sim t^{-p}\sin(\mu t)$ where the exponent $p=l+3/2$ at intermediate times, and $p=5/6$ at very late times. Our results are consistent with analytical predictions Hod and Piran (1998); Koyama and Tomimatsu (2001, 2002); Witek et al. (2013).

Finally, we note that if the quasibound states described above did not exist, the evolution of initial data would lead, on relatively short timescales, to the total depletion of scalar field everywhere. It would quickly be absorbed by the BHs or dispersed to infinity.