Science potential for stellar-mass black holes as neighbours of Sgr A*

We implement the Fisher analysis formalism to estimate the statistical error on various parameters in the waveform. Such analysis assumes the noise of the detector to be Gaussian and stationary, which works better in the limit of high SNR. Let us denote the detector’s noise as $n(t)$; then the detector output can be expressed in terms of noise and signal as

We define the inner product between two quantities $A$ and $B$ in the following manner:

$\tilde{A}(f)$ is the Fourier transform corresponding to $A$, and the asterisk superscript denotes complex conjugate. $P_{n}(f)$ is the one-side detector noise spectral density for which we consider the LISA detector noise given in Ref. Robson et al. (2019) (we ignore the galactic confusion noise as it is small compared to the instrumental for the frequency range we are working in).

The SNR of an event is

while a Fisher matrix element is defined as

In above equation, we denote a derivative with respect to a binary parameter $\theta^{i}$ as $\partial_{i}$. For a Sgr A* EMRI, the orbital frequency can be considered as almost constant during the observation period, so that it is more convenient to use the time-domain waveform which is real-valued. Hence, the Fisher matrix can be simplified in the following way:

where we have considered the contribution from two LISA channels $\{I,II\}$, and approximated the Sgr A* EMRI as a monochromatic source with constant frequency $f$. In practice, the upper limit of above integral is taken to be the observation period which is 10 years in our case. Let us denote the 1-$\sigma$ error bar on Fisher parameter $\theta^{i}$ as $\Delta\theta^{i}=\theta^{i}_{t}-\theta^{i}$ with $\theta^{i}_{t}$ denoting the true value of $\theta^{i}$. If we do not impose any prior information, the root-mean-square of $\Delta\theta^{i}$ can be obtained as the square root of diagonal elements of the inverse Fisher matrix:

The waveform of a precessing Sgr A* EMRI contains nine parameters as follows:

where $\mathcal{A}=(M_{1}M_{2})/(Da)$. Strictly speaking we know the distance $D$ and the sky locations of Sgr A* $(\bar{\theta}_{s},\bar{\phi}_{s})$ from electromagnetic observations, but here we still include $(\bar{\theta}_{s},\bar{\phi}_{s})$ as Fisher variables to address the angular resolution of LISA for such systems. As the measurement error varies for different system configurations, we use the Monte-Carlo method to sample 100 different sets of underlying parameters of Sgr A* EMRI, assuming uniform sky distribution for the angles and uniform distribution for the spin of Sgr A*. The underlying parameters for $(\bar{\theta}_{s},\bar{\phi}_{s}),D,M_{1}$ are set to be their known values, and $M_{2}$ is assumed with the same value in the fiducial model ($M_{2}=20M_{\odot}$). For each set of the underlying parameters we perform the Fisher analysis to compute the measurement uncertainties for various parameters. The histograms summarizing these Fisher analysis results are presented in Sec. II.3.3.

In order to address the detectability of Sgr A* EMRIs, we need both the expected SNR of these systems and the SNR detection threshold. The first quantity is partially discussed in Gourgoulhon et al. (2019) and in Sec. II.3.3, and in this section we compute the threshold SNR, following similar strategy described in Ref. Moore et al. (2019) developed for stellar-mass black hole binaries. Generally speaking, the threshold SNR reflects our tolerance on the false-alarm probability with suitable detection probability. It depends on the number of templates in the template bank against which the observed signal needs to be matched. Let us denote the number of waveforms in the bank by $N_{B}$, and write the template waveforms as $h_{a}(t)=\rho\hat{h}_{a}(t)$ such that $\langle\hat{h}|\hat{h}\rangle=1$, where $a=1,2,\ldots,N_{B}$.

The statistics $\sigma_{a}:=\langle s|h_{a}\rangle$ between the data $s$ and the templates $h_{a}$ follows the following probability distribution:

where $I_{0}$ is the modified Bessel function of the first kind of order $0$. If there is no signal in the data ($s=n$), the probability distribution is simply

We can claim a detection if $\sigma_{a}>\sigma_{\rm thr}$ at least for one template $a$, with $\sigma_{\rm thr}$ denoting a certain threshold determined by a conventional false alarm possibility $P_{F}$

i.e., $\sigma_{\rm thr}=\sqrt{-2\ln P_{F}}$. Following Ref. Moore et al. (2019), we choose a representative value of $P_{F}=10^{-3}/N_{B}$, which is similar to the threshold $P_{F}$ in Abbott et al. (2016). The detection probability as a function of SNR for a given $\sigma_{\rm thr}$ is given by

and the threshold SNR $\rho_{\rm thr}$ is obtained for a given detection probability of $P_{D}(\rho_{\rm thr})$ which is taken as $0.9$ here.

We now implement Fisher analyses to estimate the number of templates in the bank $N_{B}$. In our case, the first 3 models parameters (the amplitude $\mathcal{A}$ and the source location $\bar{\theta}_{s},\bar{\phi}_{s}$) do not affect $N_{B}$, because the amplitude which only shows up as a prefactor of each template, therefore can be searched over each template at negligible cost and the source location is well known. Let us define the following Fisher matrix for the remaining 6 model parameters $(\bar{\theta}_{sp},\bar{\phi}_{sp},\alpha_{p,0},\lambda_{p},\Omega_{p},{\color[rgb]{0,0,0}f})$:

Let us also define a prior range on a Fisher parameter $\theta^{i}$ by $\delta\theta^{i}$, and a modified Fisher matrix $\tilde{\Gamma}$,

Then $N_{B}$ can be approximated as the following integral over the Fisher variables:

Eq. (27) can be evaluated using the Monte-Carlo integration method where the input Fisher parameter $\theta^{i}$ is distributed over the prior range $\delta\theta^{i}$ Cornish and Porter (2005):

Here $\mathrm{Det}\,\tilde{\Gamma}_{\beta}$ is the determinant of the of the matrix $\tilde{\Gamma}$ for sample $\beta$, $N$ is the number of samples, and $V$ is the volume of the search space.

We first examine the threshold SNR for a detection of GWs from the Sgr A* EMRI. As explained in Sec. II.3.2, for detection purpose, we only need to consider the waveform containing seven free parameters. In addition, as the cost of searching over amplitude is negligible, we consider only the last six parameters in Eq. (20). The volume of search space is then $V=\delta\bar{\theta}_{sp}\,\delta\bar{\phi}_{sp}\,\delta\alpha_{p,0}\,\delta\lambda_{p}\,\delta\Omega_{p}\color[rgb]{0,0,0}{\delta f}$. To evaluate Eq. (28), we apply Monte-Carlo sampling by distributing the direction of spin and angular momentum uniformly over the 2-sphere, $\alpha_{p,0}$ between 0 to $2\pi$, and dimensionless spin parameter between 0.01 to 1. The lower end of spin is set to be nonzero as the Fisher matrix tends to me more singular for small spins, which affects the accuracy in computing the determinant of the Fisher matrix. We distribute the frequency $f$ uniformly between $1.7\times 10^{-5}$ Hz to $2\times 10^{-5}$ Hz, corresponding to a range of orbital separation $a=85M_{1}$ to $a=95M_{1}$. As a consistency check, We implement two independent sets of 100 samples to perform the integration in Eq. (28) and find a consistent threshold SNR $\rho_{\rm thr}\color[rgb]{0,0,0}{\approx 9.8}$.

In order to estimate the statistical distribution of measurement errors, we perform a separate Monte-Carlo study with the waveform containing all 9 parameters in Eq. (20). We distribute $\bar{\theta}_{sp}$, $\bar{\phi}_{sp}$, $\alpha_{p,0}$, $\color[rgb]{0,0,0}{f}$, and $\lambda_{p}$ in a similar manner as in the case of threshold SNR computation. However, we distribute the magnitude of dimensionless spin parameter from 0 to 1 uniformly as the calculation of determinant is not needed here. For each set of Monte-Carlo samples we compute the Fisher matrix, and collect all the data in histograms.

For the $M_{2}$ and range of $a$ assumed in this fiducial model, we find 80% of the samples produce SNR greater than the threshold SNR $\rho_{thr}$, while the median SNR is about 12 (see Fig. 4). Note that the SNR linearly scales as $M_{2}$ if it is different from $20M_{\odot}$, and its dependence on $a$ is discussed in Gourgoulhon et al. (2019). We also present a histogram distribution of the relative errors in measuring $\Omega_{p}$ in Fig. 5, which shows that we can constrain $\Omega_{p}$ (and hence the magnitude of spin) with a median relative error of 2%. The direction of spin can also be measured to percent level accuracy, as shown in Fig. 6 and 7. This level of accuracy is more than an order of magnitude better than the spin measurement accuracy of radio interferometry, which has to account for all the modelling uncertainties of the accretion flow. In addition, Fig. 2 and Fig. 3 show that the LISA resolution uncertainties for the sky angles are on the level of $\Delta\bar{\theta}_{s}\sim 0.033,\Delta\bar{\phi_{s}}\sim 0.022$ radian, or $\Delta\Omega/4\pi\sim 6\times 10^{-5}$. This means that the chance of having a double white dwarf source having similar signal strength (less than parsec scale distance) and sky locations is negligible. The expected measurement error of the opening angle $\lambda_{p}$ and initial precession angle $\alpha_{p}$ are shown in Fig.  8, and 9. On the other hand, the frequency can be measured very precisely with an error bar of the order $10^{-10}$ Hz.

To study the nonrelativistic description of the superradiant cloud, we will work in spherical coordinate systems $(r,\theta,\phi)$ and Cartesian coordinate $(x,y,z)$ with Sgr A* located at the origin. We consider a free complex scalar field, in which case the cloud is axisymmetric around the direction of spin of Sgr A* Zhang and Yang (2020); Baumann et al. (2019). In the nonrelativistic limit, the axion field can be expressed as

$(n,l,m)$ label the stationary eigenmodes of the field and $\omega$ is the eigenfrequency. The radial function $R_{nl}$ takes the form of that of hydrogen atom, which we can conveniently write as

Here we introduce the parameter $\tilde{x}=\alpha^{2}r/M_{1}$ with $\alpha=M_{1}\mu$ and $\mu$ is the mass of the scalar field. $L_{n-\ell-1}^{2\ell+1}\left[\frac{2\tilde{x}}{n}\right]$ is the generalized Laguerre polynomial of degree $(n-l-1)$ ¹¹1We implement the convention where $L_{n}^{a}(x)$ satisfies the differential equation $xy^{\prime\prime}+(a+1-x)y^{\prime}+ny=0$.. We choose the constant $\mathcal{B}$ in Eq. (37) such that the total mass of the cloud $M_{c}$ is normalized to $\alpha M_{1}$, namely,

Such an approximation is reasonable when $\alpha\ll 1$ and the initial BH spin is close to the maximal spin Baumann et al. (2019). Using Eq. (39) we then find $\mathcal{B}=M_{1}$.

In the next step we derive the gravitational potential generated by the bosonic cloud. Since the potential satisfies the Poisson equation, we decompose it in terms of spherical harmonics in the following manner Poisson and Will (2014); Ferreira et al. (2017):

with

We further define the harmonic components of the density $\rho$ as

For simplicity we restrict ourselves to the fastest growing mode, which corresponds to $n=2$, $l=1$, and $m=1$. For this particular mode, using Eq. (37) and Eq. (39) to Eq. (42), non-zero contributions to cloud density $\rho$ come from

Finally, using Eq. (43) and Eq. (41) to Eq. (40), we obtain the gravitational potential of the cloud as

with

and

It is convenient to express the components of the torque applied on Sgr A* EMRI orbit using the Cartesian coordinates $(x,y,z)$. Implementing the Newtonian potential in Eq. (44), the torque on object $M_{2}$, is

Such a torque will cause the orbit of $M_{2}$ to precess around the direction of spin of Sgr A*. However, since the precession frequency is expected to be much smaller than the orbital frequency, we are interested in the torque averaged over one period. In this case the point mass can be replaced by a constant density ring of the same radius to receive the gravitational torque. Let us choose the z-axis to be along the direction of spin and the x-axis to be along the ascending node of the binary. We also assume at $t=0$ object $M_{2}$ is located along x-axis. The position of $M_{2}$ can be expressed as

Here $\Omega$ is the orbital angular frequency and $\lambda_{p}$ is the angle between z-axis and the orbital angular momentum. Applying above equation to Eq. (III.2.2) and averaging over one orbital period we obtain

Next, we can express the precession frequency as a small change of precession angle $\Delta\Theta$ over small time $\Delta t$ as

Using Eq. (49) to above equation leads to a precession frequency of

which depends on $\alpha$ through $P_{2}$. We need to choose a suitable range for $\alpha$, which we discuss next.

To allow a detectable effect for Sgr A* EMRI, the growth period of the superradiant cloud needs to be much smaller than the Hubble time. In the limit of $\alpha\ll 1$, the growth period of the cloud is given by the Detweiler’s approximation Baumann et al. (2019); Detweiler (1980) :

If we set $\tau_{g}$ to be $\mathcal{O}(10^{8})$ years, above equation suggests $\alpha\simeq 0.04$, which we will consider as a lower limit on $\alpha$. On the other hand, the cloud may deplete due to GW radiation if the field is real-valued. Notice that in this case the calculation for precession rate still applies, as the cloud rotates faster than the precession rate, and a rotation-averaged density profile for real-valued field is the same as that for complex-valued field. Taking into account the GW radiation, the lifetime of the cloud can be approximated by the following expression which is valid for $\alpha<0.1$ Baumann et al. (2019):

If we choose an $\alpha$ of 0.1, which is the limit of validity of above expression, the cloud survives about 1 billion years for the Sgr A* EMRI. Hence it is reasonable to assume 0.1 as the upper bound on $\alpha$. A range of $\alpha\in[0.04,0.1]$ corresponds to the mass of the scalar field $\mu\in[13,32.6]\times 10^{-19}$ eV.

In Fig. 10 we present $\Omega_{p}$ as a function of $\alpha$ for $\lambda_{p}=\pi/4$, and also compare it with spin induced precession with same opening angle and a spin of $S=0.5M_{1}^{2}$. We find that the dark matter induced precession frequency is overall smaller than that of spin-induced precession, although they are almost comparable near $\alpha=0.1$. Nevertheless, since in both cases the precession happens around the direction of spin, the modification to GW phase is similar and in principle it is not possible to distinguish the two effects.

In addition to the orbital precession caused by the Newtonian potential of the axion cloud, a dynamical friction force is also applied on object $M_{2}$, which can be expressed as Hui et al. (2017); Chandrasekhar (1943)

In case of Sgr A* EMRI where the orbit is circular, the velocity $v$ of $M_{2}$ relative to the wave function can be approximated as the orbital velocity, and $C_{\Lambda}$ has a simplified form Zhang and Yang (2020)

where $k=\mu v$ with $\mu$ representing the mass of the scalar field. We again restrict ourselves to the $n=2$, $l=1$, and $m=1$ mode and use Eq. (55) and $\rho=\Psi^{*}\Psi$ to Eq. (54) to obtain the rate of change of energy loss due to dynamical friction as $\dot{E}_{DF}=Fv$. Similar to III.1, Energy loss due to dynamical friction and GW emission in this case can be taken as small and constant, which enters the GW phase as $\delta\varphi=\dot{f}\pi t^{2}$. The waveform we consider in this case is

where $\varphi_{p}(t)$ is the polarization phase and $\varphi_{D}(t)$ is the Doppler phase as in section II.1, and

Quantities with subscript “0” are evaluated for a circular binary with $a=90M$. The second term on the right side is the contribution to frequency evolution from dynamical friction: $\dot{f}_{DF}$.

In Fig. 11 we present $\dot{f}_{DF}$ as a function of $\alpha$ for a Sgr A* EMRI in the equatorial plane $\theta=\pi/2$, as well as the 1-$\sigma$ upper bound on $\dot{f}_{DF}$ obtained from Fisher analyses. The Fisher analysis is performed with the parameters $\ln\mathcal{A}$, $\bar{\theta}_{s}$, $\bar{\phi}_{s}$, $\bar{\theta}_{L}$, $\bar{\phi}_{L}$, $\color[rgb]{0,0,0}{f_{0}}$, and $\dot{f}_{DF}$. Here ($\bar{\theta}_{L},\bar{\phi}_{L}$) is the direction of orbital angular momentum of Sgr A* EMRI, for which we consider ($\pi/4,3\pi/4$) as a fiducial value. The error of $\dot{f}_{DF}$ is not sensitive to the fiducial value of $\dot{f}_{DF}$. We find the 1-$\sigma$ upper bound on $\dot{f}_{DF}$ is approximately $\color[rgb]{0,0,0}{3.47\times 10^{-18}}\,\mathrm{s}^{-2}$, which roughly corresponds to $\color[rgb]{0,0,0}{\alpha=0.095}$, which means an effect to the waveform phase coming from $\color[rgb]{0,0,0}{\alpha\in[0.095,0.1]}$ or $\color[rgb]{0,0,0}{\mu\in[3.1,3.3]\times 10^{-18}}$ eV can exceed the statistical uncertainty. The chance that axion mass lies within such narrow range is small. Fig. 11 also shows that the dynamical friction due to cold DM spike around Sgr A* can produce $\dot{f}_{DF}$ which is comparable to the Axion DM case only when $\alpha\leq 0.054$. Note that this is an overestimate, as we assumed static cold dark ($\xi=1$) while in reality $\xi$ is expected to be smaller Coogan et al. (2021); Kavanagh et al. (2020). On the other hand, $\dot{f}_{DF}$ in Axion DM case can be greater than the GW radiation induced $\dot{f}$ if $\alpha>0.087$.