Prospects of Finding Detached Black Hole–Star Binaries with TESS

• Ellipsoidal variation: Tidal force due to BH's gravity changes the geometric shape of the star as well as brightness distribution on the stellar surface via gravity darkening (von Zeipel 1924). This causes a light modulation with an amplitude of

  $$\begin{eqnarray}\begin{array}{rcl}{s}_{\mathrm{ev}} & = & {\alpha }_{\mathrm{ev}}\,\displaystyle \frac{{M}_{\bullet }\sin i}{{M}_{\star }}{\left(\displaystyle \frac{{R}_{\star }}{a}\right)}^{3}\sin i\\ & = & 1.89\times {10}^{-2}{\alpha }_{\mathrm{ev}}{\sin }^{2}i{\left(\displaystyle \frac{P}{1\mathrm{day}}\right)}^{-2}\\ & & \times \ {\left(\displaystyle \frac{{\rho }_{\star }}{1{\rm{g}}{\mathrm{cm}}^{-3}}\right)}^{-1}\left(\displaystyle \frac{1}{1+{M}_{\star }/{M}_{\bullet }}\right).\end{array}\end{eqnarray} \tag{ 1 }$$

  Here,

  $$\begin{eqnarray}&&{\alpha }_{\mathrm{ev}}=0.15\,\displaystyle \frac{(15+u)(1+g)}{3-u},\end{eqnarray} \tag{ 2 }$$

  where g is the gravity-darkening coefficient and u is the linear limb-darkening coefficient (Morris & Naftilan 1993). For the circular orbit, the modulation is symmetric with respect to our line of sight, and the signal has the period half the orbital one. In this paper, we neglect the dependence on the eccentricity and set ${\alpha }_{\mathrm{ev}}=1$ just for simplicity, since our purpose is not to give a precise prediction for the expected yield.Figure 1 shows that, at a fixed orbital period, the amplitude of the EV signal saturates with increasing companion masses. This is because ${a}^{3}\sim {M}_{\bullet }$ when ${M}_{\bullet }\gg {M}_{\star }$ and this dependence cancels the companion mass dependence of the tidal force (Equation (1)). Hence it will be difficult to precisely weigh the most massive companions with the EV amplitude alone.
• Doppler beaming: Relativistic aberration of light, time dilation, and Doppler shift caused by the primary's radial acceleration produces the photometric variation with the amplitude of (Loeb & Gaudi 2003)

  $$\begin{eqnarray}\begin{array}{rcl}{s}_{\mathrm{beam}} & = & {\alpha }_{\mathrm{beam}}\,4\,\displaystyle \frac{{K}_{\star }}{c}\\ & = & 2.8\times {10}^{-3}{\alpha }_{\mathrm{beam}}\sin i{\left(\displaystyle \frac{P}{1\mathrm{day}}\right)}^{-1/3}\\ & & \times \ {\left(\displaystyle \frac{{M}_{\bullet }+{M}_{\star }}{{M}_{\odot }}\right)}^{-2/3}\left(\displaystyle \frac{{M}_{\bullet }}{{M}_{\odot }}\right),\end{array}\end{eqnarray} \tag{ 3 }$$

  where ${\alpha }_{\mathrm{beam}}$ is given by integrating the following wavelength-dependent factor

  $$\begin{eqnarray}&&{\alpha }_{\mathrm{beam},\nu }=\displaystyle \frac{1}{4}\left(3-\displaystyle \frac{d\,\mathrm{log}\,{F}_{\nu }}{d\,\mathrm{log}\,\nu }\right)\end{eqnarray} \tag{ 4 }$$

  over the photometric bandpass, and ${F}_{\nu }$ is the spectrum of the star. Again, we set ${\alpha }_{\mathrm{beam}}=1$, which is reasonable for a sunlike star (e.g., Shporer 2017). This signal is anti-phased with the radial velocity variation of the star and has a different orbital-phase dependence from the ellipsoidal variation (see Figure 2). Thus, if both the EV and beaming signals are detected, the consistency of their phases and amplitudes will be useful to confirm the binary origin. An independent measurement of the spectroscopic orbit also helps the interpretation of the light curve because the shape and amplitude of the beaming component are essentially fixed.
• Self lensing: If the system is eclipsing, the BH passing in front of the stellar disk acts as a lens to gravitationally magnify the background star (e.g., Gould 1995; Sahu & Gilliland 2003; Rahvar et al. 2011). The light curve exhibits periodic pulses with the height (Trimble & Thorne 1969; Maeder 1973; Agol 2003):

  $$\begin{eqnarray}\begin{array}{rcl}{s}_{\mathrm{sl}} & = & 2{\left(\displaystyle \frac{{R}_{{\rm{E}}}}{{R}_{\star }}\right)}^{2}=7.15\times {10}^{-5}{\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)}^{-2}\\ & & \times \ {\left(\displaystyle \frac{P}{1\mathrm{day}}\right)}^{2/3}\left(\displaystyle \frac{{M}_{\bullet }}{{M}_{\odot }}\right){\left(\displaystyle \frac{{M}_{\bullet }+{M}_{\star }}{{M}_{\odot }}\right)}^{1/3},\end{array}\end{eqnarray} \tag{ 5 }$$

  where

  $$\begin{eqnarray}\begin{array}{rcl}{R}_{{\rm{E}}} & = & \sqrt{\displaystyle \frac{4{{GM}}_{\bullet }a}{{c}^{2}}}=4.27\times {10}^{-2}{R}_{\odot }{\left(\displaystyle \frac{P}{1\mathrm{yr}}\right)}^{1/3}\\ & & \times \ {\left(\displaystyle \frac{{M}_{\bullet }}{{M}_{\odot }}\right)}^{1/2}{\left(\displaystyle \frac{{M}_{\bullet }+{M}_{\star }}{{M}_{\odot }}\right)}^{1/6}\end{array}\end{eqnarray} \tag{ 6 }$$

  is the Einstein radius. The duration of the signal is

  $$\begin{eqnarray}\begin{array}{rcl}{\tau }_{\mathrm{sl}} & = & \displaystyle \frac{{R}_{\star }P}{\pi a}\cdot \displaystyle \frac{\pi }{4}=1.8\,\mathrm{hr}\times \displaystyle \frac{\pi }{4}{\left(\displaystyle \frac{P}{1\mathrm{day}}\right)}^{1/3}\\ & & \times \ {\left(\displaystyle \frac{{M}_{\bullet }+{M}_{\star }}{{M}_{\odot }}\right)}^{-1/3}\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)\end{array}\end{eqnarray} \tag{ 7 }$$

  when averaged over the impact parameter, and the signal vanishes if $\cos i\gt {R}_{\star }/a$ for a circular orbit. The pulses have the same period as the orbital period.Since stellar light is also blocked by the eclipsing compact object, the pulses are not observed unless $\sqrt{2}{R}_{{\rm{E}}}$, which increases with the binary separation, is larger than the physical size of the object. This is why all of the known self-lensing binaries containing WDs (Kruse & Agol 2014; Kawahara et al. 2018; Masuda et al. 2019) have orbital periods ranging from months to years, longer than typical eclipsing systems. On the other hand, the physical radius has negligible effects for BHs and NSs, and so they always show pulses regardless of the orbital period. Depending on the detected period the WD case can be excluded by sheer presence of pulses.If the orbital period is measured from multiple pulses, the BH mass is derived from the pulse height and the stellar radius. For the latter, the parallax information from Gaia has already allowed for 10%-level measurements for tens of millions of stars (Andrae et al. 2018). If combined with the spectroscopic effective temperature, the precision better than five percent can be achieved (Berger et al. 2018); this translates into 10%-precision for the BH mass, provided that the pulse height is sufficiently well constrained.

The main purpose of the TESS mission is to find transiting planets around nearby stars. TESS employs four cameras each with a field of view (FOV) of 24° × 24°, which point each sector of the sky for $27.4\,\mathrm{days}$ (two spacecraft orbits). Over its two-year mission, 26 sectors with partial overlaps will be observed to tile almost the entire sky excluding the region near the ecliptic; the region near the ecliptic poles will be most densely covered and continuously observed for one year. TESS will provide images stacked at two minutes cadence for a few 10⁵ targets pre-selected for planet search and asteroseismology, as well as 30 minute cadence full-frame images for all sources.

To enable the selection of optimal targets for the planet search, a catalog of luminous sources on the sky, TESS Input Catalog (TIC, Stassun et al. 2017), has been created. We adopt the magnitude in the TESS bandpass, stellar mass ${M}_{\star }$, stellar radius ${R}_{\star }$, and stellar effective temperature ${T}_{\star }$ in the TIC version 6 (Stassun et al. 2017). We focus on the objects classified as stars, on ecliptic latitude $| l| \gt 6^\circ$ (i.e., outside the gap in the TESS FOV), and with both mass and radius estimated in the catalog. The last selection limits the targets from ∼400 million to ∼20 million low-mass stars (mostly $\lt 2.5\,{M}_{\odot }$) with TESS magnitudes brighter than 15. We do not consider this to be a serious limitation because the quality of the TESS photometry will be limited for fainter stars in any case.

We use the fitting formula in Stassun et al. (2017, ver. 20170531) to evaluate the white component of the expected photometric noise as a function of the TESS magnitude presented in Sullivan et al. (2015). We adopt the observing duration $T=27.4\,\mathrm{days}$ (i.e., minimum observing duration) for all of the stars to give a conservative estimate. We always assume 30 minutes cadence observations to deal with all of the stars that will be in the full-frame images.

If the stellar rotation is synchronized with the orbit, the star may show enhanced spot activity with a similar period to that of the self-lensing pulses. Since the duration of self-lensing pulses is much shorter than the orbital period (and hence rotation period of the spin-synchronized star), the rotational modulation can basically be filtered out without affecting the self-lensing signal. However, this removal may not be perfect, and also stellar variabilities may enhance correlated noise on a timescale similar to the duration of self-lensing events. Thus, estimates considering the uncorrelated noise alone (as described in Section 3.1) most likely overestimate the detectability of the self-lensing signal.

To evaluate the detectability in the presence of such effects, we perform injection-and-recovery simulations using the actual light curves of spotted stars from the prime Kepler mission. Because Kepler has a ∼10× larger aperture than TESS, the Kepler light curves include much smaller photon noise and can be considered as a realistic template of astrophysical variability of the stars. Therefore we can simulate the light curves of spotted stars observed with TESS by adding white noise component due to its smaller aperture (as described in Section 3.1) to the actual Kepler light curves.³ Here we assume that the stellar rotation is always synchronized with the orbit of the self-lensing systems, and inject self-lensing signals with the same period as the rotation period of each Kepler star.

We take random segments of $27.4\,\mathrm{days}$ long from the Kepler light curves of rotating stars and detrend them with median filters with the length of $0.25\times {(P/30)}^{1/3}\,\mathrm{days}$, which is typically a few times longer than self-lensing durations and much shorter than rotation period. Then, we run the box least-squares (BLS) algorithm⁴ (Kovács et al. 2002) on these light curves, after randomly injecting Gaussian noise with log-uniform dispersions spanning from 100 to 10,000 ppm (to cover the expected noise range from TESS) and flipping the light curves around their median values (to evaluate any possible source of positive excursion). We then set the threshold value for the peak signal-to-noise of the BLS spectrum, so that no detection is found from 10⁶ trials to achieve a false-positive rate less than 10⁻⁶. Given that ∼10% of sunlike stars exhibit spot modulation detectable with Kepler (McQuillan et al. 2014), this yields the FP rate of ∼10⁻⁷, which corresponds to ${ \mathcal O }(1)$ FP in our sample of 20M stars. Considering the strong dependence of the spot activity on the rotation period, we choose different thresholds for each period bin, and the resulting thresholds turn out to be higher for shorter-period systems. We also count the BLS peaks as detections only if the detected period is close to the measured rotation period of the stars.⁵ Then, we perform the simulations injecting both Gaussian noise and self-lensing pulses of various amplitudes (100–10⁵ ppm; see Figure 1) and phases, and check what fraction of the injected signals are recovered using the thresholds as defined above. We fit the resulting recovery rate as a function of the signal-to-noise of the phase-folded signal $\sqrt{T/P}({s}_{\mathrm{sl}}/{\sigma }_{30\mathrm{minutes}})$ using the gamma cumulative distribution function (Fulton et al. 2017) scaled by the asymptotic value of the recovery rate at high signal-to-noise ratios: the values were ≈100% for P = 0.8–$4.8\,\mathrm{days}$, ≈90% for P = 0.3–$0.8\,\mathrm{days}$, and P = 4.8–$11.9\,\mathrm{days}$, and $\approx 20 \%$ for P = 11.9–$30\,\mathrm{days}$. The resulting function ${C}_{\mathrm{rec}}(\sqrt{T/P}({s}_{\mathrm{sl}}/{\sigma }_{30\mathrm{minutes}}))$ is used to evaluate the expected number of searchable stars.

For a star with mass ${M}_{\star }^{j}$ and radius ${R}_{\star }^{j}$, the effective searchability for a BH/NS companion with period P and mass ${M}_{\bullet }$ on randomly oriented orbits is evaluated as

$$\begin{eqnarray}&&{N}_{\mathrm{eff}}(P,{M}_{\bullet },{M}_{\star }^{j},{R}_{\star }^{j})=\displaystyle \frac{{R}_{\star }}{a}\,{C}_{\mathrm{rec}}\left(\sqrt{n}\,\displaystyle \frac{{s}_{\mathrm{sl}}}{{\sigma }_{30\mathrm{minutes}}}\right),\end{eqnarray} \tag{ 8 }$$

where ${R}_{\star }$/a corresponds to the eclipse probability for a circular orbit, and n is the number of pulses in the data. The number n can be either ${n}_{0}\equiv [T/P]$ or ${n}_{0}+1=[T/P]+1$ depending on the orbital phase, where $[x]$ is the greatest integer that does not exceed x. The corresponding probabilities are ${p}_{{n}_{0}}\,=[T/P]+1-T/P$ and ${p}_{{n}_{0}+1}=T/P-[T/P]$, respectively. We count the number of searchable stars for both n₀ and ${n}_{0}+1$, and average them with the weights ${p}_{{n}_{0}}$ and ${p}_{{n}_{0}+1}$. We set ${p}_{{n}_{0}=1}=0$ to count only the pulses that are observed at least twice; although, such longer-period systems turned out to be difficult to search in any case.

The interference with stellar activity is a more serious issue for the phase-curve signal, because the star-spot modulation includes signals with both periods $P/2$ and P, if the rotation is synchronized with the orbit. This means that a simple periodicity search will introduce a large number of FPs from spot modulations, because there exist more spotted stars than the BH/NS systems. Here, we consider a search focusing on EV/beaming signals that are likely too large to be produced by a spotted star with a given rotation period and mass, and examine if the search could provide reasonable number of candidates without being completely swamped by false positives.

To set the threshold, we compute Fourier transform of the Kepler light curves of all of the spotted stars in McQuillan et al. (2014) and measure the amplitudes of the peaks corresponding to the rotation period s_P and its first harmonic ${s}_{P/2}$. We then divide the sample into the period bins as we used for orbital periods, as well as temperature bins ($\lt 4000$, 4000–5000, 5000–6000, $\gt 6000\,{\rm{K}}$), and take the maximum values of s_P and ${s}_{P/2}$ in each bin. We set the threshold for the EV and beaming signals, ${s}_{\mathrm{ev},\mathrm{th}}$ and ${s}_{\mathrm{beam},\mathrm{th}}$, so that the Fourier amplitudes from the EV and beaming signals are larger than ${s}_{P/2}$ and s_P, respectively. This allows us to pick up phase-curve signals that are likely too strong to be caused by spot modulations, taking into account the dependence of their amplitudes on the stellar rotation period and effective temperature. These thresholds are much larger than the noise level and statistical false positives are negligible. Thus, assuming that sunlike stars observed by TESS have similar properties as the Kepler sample, the FP rate due to spotted stars would roughly be less than one in 10⁵. The rate corresponds to some 100 FPs in our 20M sample stars.⁶ Although this number of possible FPs is not small, if a similar (or larger) number of BH/NS candidates can be detected with this choice of thresholds, the search would still provide meaningful information to select promising candidates, especially considering that further vetting (e.g., presence or absence of any RV variations) would readily be possible for many stars using the archival spectroscopic data (e.g., Thompson et al. 2018; Gu et al. 2019). We will see that this can be the case in Section 4.2.

Following Equations (1) and (3), the above detection criteria, ${s}_{\mathrm{ev}}\gt {s}_{\mathrm{ev},\mathrm{th}}$ and ${s}_{\mathrm{beam}}\gt {s}_{\mathrm{beam},\mathrm{th}}$ translate into

$$\begin{eqnarray}&&\begin{array}{l}1-{\cos }^{2}i\gt \displaystyle \frac{{s}_{\mathrm{ev},\mathrm{th}}}{{s}_{\mathrm{ev},\cos i=0}}\equiv {\beta }_{\mathrm{ev}},\\ \sqrt{1-{\cos }^{2}i}\gt \displaystyle \frac{{s}_{\mathrm{beam},\mathrm{th}}}{{s}_{\mathrm{beam},\cos i=0}}\equiv {\beta }_{\mathrm{beam}},\end{array}\end{eqnarray} \tag{ 9 }$$

where the subscript $\cos i=0$ refers to the amplitude for a system with $\cos i=0$. Thus, averaging over random orbital inclination (i.e., uniform $\cos i$), the effective searchability of each star is given by $H(1-{\beta }_{\mathrm{ev}})\sqrt{1-{\beta }_{\mathrm{ev}}}$ for the EV signal and $H(1-{\beta }_{\mathrm{beam}})\sqrt{1-{\beta }_{\mathrm{beam}}^{2}}$ for the beaming signal, where H is the Heaviside step function ($H(x)=1$ for $x\gt 0$ and 0 otherwise). Considering that either of the signal, if above the threshold, would be sufficient to flag the system to be anomalous, we take the larger of the two values to evaluate the searchability of each star

$$\begin{eqnarray}&&\begin{array}{l}{N}_{\mathrm{eff}}(P,{M}_{\bullet },{M}_{\star }^{j},{R}_{\star }^{j},{T}_{\star }^{j})=\max \left[H(1-{\beta }_{\mathrm{ev}})\sqrt{1-{\beta }_{\mathrm{ev}}},\right.\\ \quad \left.H(1-{\beta }_{\mathrm{beam}})\sqrt{1-{\beta }_{\mathrm{beam}}^{2}}\right].\end{array}\end{eqnarray} \tag{ 10 }$$

Since we focus on detached systems, we exclude the stars filling their Roche lobes. We compute the effective Roche-lobe radius for the star by (Eggleton 1983)

$$\begin{eqnarray}&&{R}_{{\rm{L}}}(q,a)=\displaystyle \frac{0.49{q}^{2/3}}{0.6{q}^{2/3}+\mathrm{ln}(1+{q}^{1/3})}\,a\equiv {\tilde{R}}_{{\rm{L}}}(q)a,\end{eqnarray} \tag{ 11 }$$

where $q={M}_{\star }/{M}_{\bullet }$. We require that ${R}_{\star }\lt {R}_{{\rm{L}}}$, or $a\gt {R}_{\star }/{\tilde{R}}_{{\rm{L}}}$. We also exclude the systems where the GW emission causes rapid orbital decay compared to the main-sequence lifetime of the star. We compute the decay time as ${t}_{\mathrm{GW}}=3.3\,\times {10}^{8}\,\mathrm{Gyr}\times \tfrac{{(a/\mathrm{au})}^{4}}{{M}_{\bullet }{M}_{\star }({M}_{\bullet }+{M}_{\star })/{M}_{\odot }^{3}}$ (Peters 1964) and the main-sequence lifetime as ${t}_{\mathrm{MS}}=10\,\mathrm{Gyr}{({M}_{\star }/{M}_{\odot })}^{-2.5}$, and omit the systems with ${t}_{\mathrm{GW}}\lt {t}_{\mathrm{MS}}/2$. Taking all of these into account, the number of searchable stars for a set of $(P,{M}_{\bullet })$ is computed by

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}{N}_{{\rm{s}}{\rm{e}}{\rm{a}}{\rm{r}}{\rm{c}}{\rm{h}}{\rm{a}}{\rm{b}}{\rm{l}}{\rm{e}}}(P,{M}_{\bullet })=\displaystyle \sum _{j}{N}_{{\rm{e}}{\rm{f}}{\rm{f}}}(P,{M}_{\bullet },{M}_{\star }^{j},{R}_{\star }^{j},{T}_{\star }^{j})\\ \,\cdot \,H({a}^{j}-{R}_{\star }^{j}/{\tilde{R}}_{{\rm{L}}}^{j})\cdot H({t}_{{\rm{G}}{\rm{W}}}^{j}-{t}_{{\rm{M}}{\rm{S}}}^{j/2}),\end{array}\end{array}\end{eqnarray} \tag{ 12 }$$

where the index j runs over all of the stars and ${N}_{\mathrm{eff}}$ is computed for each of the self-lensing signal (Equation (8)) and phase-curve signal (Equation (10)). The results are shown in Figure 3 for the self-lensing signal (left panel) and for the phase-curve signal (right panel).

4.1.1. Estimate Based on Field Binaries

Here, we construct the population of BH–star binaries assuming that their properties follow those of field binaries (see Mashian & Loeb 2017). We pick up "BHs" from a power-law mass function, ${I}_{\mathrm{BH}}({M}_{\bullet })\propto {M}_{\bullet }^{-2.3}$, with a minimum mass of $5\,{M}_{\odot }$. The mass function is normalized such that all of the stars with $\geqslant 20{M}_{\odot }$ end up as BHs. This gives the occurrence of BHs as

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{\mathrm{BH}}}{{{dM}}_{\bullet }}=H({M}_{\bullet }-5\,{M}_{\odot }){I}_{\mathrm{BH}}({M}_{\bullet }).\end{eqnarray} \tag{ 14 }$$

Then, we assign companion stars to these BHs, based on the occurrence of massive star binaries in the field (e.g., Sana et al. 2012):

$$\begin{eqnarray}&&{f}_{\mathrm{companion}}(q,P){dq}\,{dP}=C\,\displaystyle \frac{{q}^{0}}{P}\,{dq}\,{dP},\end{eqnarray} \tag{ 15 }$$

where $q\leqslant 1$ is the binary mass ratio, P is the orbital period, and C is the normalization constant that fixes the binary fraction (see below). Thus, the occurrence of BHs with stellar companions is

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{dN}}_{\mathrm{BH} \mbox{-} \mathrm{star}}}{{{dM}}_{\bullet }\,{{dM}}_{\star }\,{dP}}=\displaystyle \frac{1}{{M}_{\bullet }}\displaystyle \frac{{{dN}}_{\mathrm{BH} \mbox{-} \mathrm{star}}}{{{dM}}_{\bullet }\,{dq}\,{dP}}\\ \quad =\ \displaystyle \frac{1}{{M}_{\bullet }}\,\displaystyle \frac{{{dN}}_{\mathrm{BH}}}{{{dM}}_{\bullet }}\,{f}_{\mathrm{companion}}\left(\displaystyle \frac{{M}_{\star }}{{M}_{\bullet }},P\right)=\displaystyle \frac{C}{{M}_{\bullet }}\,\displaystyle \frac{{{dN}}_{\mathrm{BH}}}{{{dM}}_{\bullet }}\,\displaystyle \frac{1}{P}.\end{array}\end{eqnarray} \tag{ 16 }$$

Here, we assume circular orbits for simplicity. This yields

$$\begin{eqnarray}&&\begin{array}{l}{p}_{\mathrm{field}}(P,{M}_{\bullet }| {M}_{\star })=\displaystyle \frac{{{dN}}_{\mathrm{BH} \mbox{-} \mathrm{star}}/{{dM}}_{\bullet }\,{{dM}}_{\star }\,{dP}}{{{dN}}_{\mathrm{star}}/{{dM}}_{\star }}\\ \quad =\ \displaystyle \frac{C}{{M}_{\bullet }}\,\displaystyle \frac{H({M}_{\bullet }-5\,{M}_{\odot }){I}_{\mathrm{BH}}({M}_{\bullet })}{I({M}_{\star })}\,\displaystyle \frac{1}{P},\end{array}\end{eqnarray} \tag{ 17 }$$

where ${{dN}}_{\mathrm{star}}/{{dM}}_{\star }=I({M}_{\star })$ is the initial stellar mass function from Kroupa (2001). We choose C so that the binary fraction integrated over $P=0.1\,\mathrm{days}$ to $P={10}^{3.5}\,\mathrm{days}$ is 0.5, following Sana et al. (2012). Strictly speaking, it is not clear whether the normalization remains the same for binaries containing BHs. The choice is for an optimistic scenario in which formation of a BH does not lead to any loss of binary companions.

4.1.2. Estimate Based on a Simple Model of the Common-envelope Evolution

The above estimate does not take into account any interactions in the binary. If they go through the common-envelope (CE) phase, the low-mass companion may merge with the BH progenitor. Even if it survives, the binary orbit should dramatically shrink during the process. Here, we examine the outcome of this CE evolution. Unlike in Yamaguchi et al. (2018), we do not consider the case where the mass transfer occurs stably, since our focus is on systems with initially large mass ratios.

We first follow the same procedure as in Section 4.1.1 to sample the population of binaries consisting of a BH progenitor with mass ${M}_{\bullet {\rm{i}}}$ and its stellar companion with mass ${M}_{\star }$: $p({P}_{{\rm{i}}},{M}_{\bullet {\rm{i}}}| {M}_{\star })$. Then, we use the fitting formulae provided by Hurley et al. (2000) to compute the core mass ${M}_{\bullet {\rm{i}},\mathrm{core}}$ of the BH progenitor as well as its maximum radius ${R}_{\max }$ (typically ∼1000–$3000{R}_{\odot }$) achieved during its evolution. We then assume that all of the sampled binaries with initial semimajor axes ${a}_{{\rm{i}}}={\left[{P}_{{\rm{i}}}^{2}G({M}_{\bullet {\rm{i}}}+{M}_{\star })/4{\pi }^{2}\right]}^{2/3}$ smaller than ${R}_{\max }$ goes through the CE phase. The companion survives if the orbital energy is sufficiently large to completely strip the envelope of the BH progenitor; we assume only the stripped core of the BH progenitor is left when this happens. The semimajor axis after this process, ${a}_{{\rm{f}}}$, is computed by (Webbink 1984):

$$\begin{eqnarray}&&\alpha \left(\displaystyle \frac{{{GM}}_{\bullet {\rm{i}},\mathrm{core}}{M}_{\star }}{2{a}_{{\rm{f}}}}-\displaystyle \frac{{{GM}}_{\bullet {\rm{i}}}{M}_{\star }}{2{a}_{{\rm{i}}}}\right)=\displaystyle \frac{G({M}_{\bullet {\rm{i}}}-{M}_{\bullet {\rm{i}},\mathrm{core}}){M}_{\bullet {\rm{i}}}}{\lambda {R}_{\mathrm{Roche},{\rm{i}}}},\end{eqnarray} \tag{ 18 }$$

where ${R}_{\mathrm{Roche},{\rm{i}}}={R}_{{\rm{L}}}({M}_{\bullet {\rm{i}}}/{M}_{\star },{a}_{{\rm{i}}})$ with ${R}_{{\rm{L}}}$ defined in Equation (11), and we adopt $\alpha \lambda =1$ (Belczynski et al. 2002). Finally, we assume that the remnant core of the BH progenitor loses mass via stellar wind and SN explosion before it becomes a BH. For the wind mass loss during the Wolf–Rayet star phase, we follow Vink (2017). At the core collapse, the remnant BH mass and SN ejecta mass are determined by the formula in Belczynski et al. (2002). Both mass loss processes change the semimajor axis.

The above procedures define a transformation between $({P}_{{\rm{i}}},{M}_{\bullet {\rm{i}}})$ and $(P,{M}_{\bullet })$. The final distribution after the evolution is then computed by

$$\begin{eqnarray}&&{p}_{\mathrm{CE}}(P,{M}_{\bullet }| {M}_{\star })=\left|\displaystyle \frac{\partial ({P}_{{\rm{i}}},{M}_{\bullet {\rm{i}}})}{\partial (P,{M}_{\bullet })}\right|\,{p}_{\mathrm{field}}({P}_{{\rm{i}}},{M}_{\bullet {\rm{i}}}| {M}_{\star }).\end{eqnarray} \tag{ 19 }$$

Figure 4 shows the BH occurrence rates for the field-binary model and the CE model with αλ = 1 described in Section 4.1. The values are averages of $p(P,{M}_{\bullet }| {M}_{\star }){\rm{\Delta }}{M}_{\bullet }{\rm{\Delta }}P$ over ${M}_{\star }$ of the stars searchable with self-lensing in each bin; the differences are mostly within a factor of a few when the stars searchable with phase-curve signals are used instead. Figures 5 and 6 show the estimated numbers of detectable BHs based on the occurrences in Figure 4 and the number of searchable stars in Figure 3. Both field-binary and CE models predict some 10 BH companions detectable with self-lensing, and some 100 with phase-curve variations. This is the consequence that the BH occurrence in the relevant region of the parameter space is about 10⁻⁴ (Figure 4), while the effective numbers of stars searchable with self-lensing and phase-curve signals are ∼10⁵ and ∼10⁶, respectively. Curiously, the estimated occurrence appears to be compatible with the discovery of a BH/NS system by Thompson et al. (2018) among ≳10⁵ stars from Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2017); although, the period we focus on is much shorter. While the resulting occurrences are similar between the field-binary and CE models, the detectable binaries from the CE models were initially in much wider orbits (P = 10²–10⁴ days) and surrendered their orbital energies to survive the CE evolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These calculations suggest that the number of searchable stars with TESS may indeed be sufficient to actually detect BH/NS companions. This also indicates that even an upper limit on the occurrence rate from the null detection, if quantified, will be able to provide a meaningful observational constraint on the poorly understood process of the CE evolution. We also note that the searchability increases both in terms of the number and period range, if individual stars are observed for a longer duration in the overlapping regions of the observing sectors or during a potential extended mission (Bouma et al. 2017).

In these population models, the typical mass of the stars with detectable BH companions is 1–$2\,{M}_{\star }$, mainly by construction of the sample. The brightest star with a BH companion is predicted to have a V-band magnitude of ∼10 for the self-lensing population, and will be even brighter for the phase-curve sample. These magnitudes correspond to the distance of a few tens to a few hundreds of parsecs. If actually found, these systems will therefore be among the nearest known BHs.

For a sanity check, we applied the same procedure to the Kepler sample, using stellar properties and the combined differential photometric precision (Koch et al. 2010) of the stars in the Kepler input catalog (Mathur et al. 2017). Adopting the same detection thresholds, we found that no detection is expected with either self-lensing or phase-curve signals, consistently with the null detection of such systems so far. The conclusion for the self-lensing BH systems is insensitive to the adopted threshold because the null detection is expected simply due to their low occurrence rate, as pointed out by Agol (2003). The interpretation of the null detection of phase-curve signals is more subtle. Our population models predict that there would be some 10 such systems on close-in orbits that would have been detectable in the absence of stellar variabilities: thus, either they have been swamped in (or confused with) the stellar variabilities and await confirmation with Doppler observations, or our model overestimates the number of BH systems on close-in orbits. In the latter case, the detection with TESS still remains to be plausible even if the actual rate is 10 times smaller than our estimate; if the overestimation factor is 100 or larger, the detection with the TESS data may not be feasible.

Our sample consists of 18 X-ray binaries with dynamically confirmed BH companions in Table 1 of Remillard & McClintock (2006), with two systems in the Large Magellanic Cloud being excluded. Distances to the systems are adopted from Table 4.1 of McClintock & Remillard (2006) except for GS 1354−64 (Casares et al. 2004)⁷ and XTE J1650−500 (Orosz et al. 2004).

We adopt $P=0.8\,\mathrm{days}$ and ${M}_{\bullet }=7\,{M}_{\odot }$ as the representative values of the above BHXB sample, and assume a sunlike star with ${M}_{\star }=1\,{M}_{\odot }$ and ${R}_{\star }=1\,{R}_{\odot }$ for a typical companion. For these parameters, the amplitudes of the self-lensing and phase-curve signals are $\approx 9\times {10}^{-4}$ and $\approx 6\times {10}^{-3}$, respectively, where the detectability of the latter is limited by the beaming signal (see Figure 1). Assuming that these signals are detectable when the signal-to-noise ratio of the phase-folded signal is larger than 10 for the white noise model in Section 3.1, these amplitudes are roughly translated into the limiting magnitudes of $V\approx 11$ and $V\approx 15$, respectively, or the maximum searchable distances of ${d}_{\max }\approx 0.25\,\mathrm{kpc}$ and ${d}_{\max }\approx 1.3\,\mathrm{kpc}$ for the assumed sunlike companion. These are shown as vertical bands in Figure 7. Note that ${d}_{\max }$ for self-lensing is rather sensitive to the property of the stellar companion because the signal scales as ${R}_{\star }^{-2}$ (Equation (5)). It varies from $\approx 0.1\,\mathrm{kpc}$ to $\approx 0.4\,\mathrm{kpc}$ for A to K dwarf companions. The dependence is much smaller for the beaming signal (see Equation (3)).

Zoom In Zoom Out Reset image size

Given the cumulative distance distribution of BHXBs, N(d), the value of $N({d}_{\max })$ roughly corresponds to the number of such systems whose self-lensing (if present) and phase-curve signals would be detectable with TESS, assuming no other source of optical light variations. We estimate N(d) as $N(d)\,={{fN}}_{\mathrm{confirmed}}(d)$, where ${N}_{\mathrm{confirmed}}(d)$ is the cumulative distance distribution of the above confirmed BHXB sample and f is a correction factor for the incompleteness of the X-ray binary search. Arur & Maccarone (2018) evaluated the detectability of X-ray and optical signals as well as X-ray outbursts for the above BHXB sample, and concluded that the completeness of the detection is $\approx 1/30$ on average. We show N(d) with f = 30 motivated by this result as a thin histogram in Figure 7, along with the filled histogram for ${N}_{\mathrm{confirmed}}(d)$. For $d\lesssim 1\,\mathrm{kpc}$ at which no BHXB has been detected, we extrapolate the distribution assuming $N(d)\propto {d}^{2}$, i.e., they are dominated by the disk population, and their space density is roughly constant. Although the validity of these assumptions is uncertain, this extrapolation only matters the estimate regarding the self-lensing population.

Figure 7 shows that the phase-curve signal will be detectable among several tens of BHXB systems for f = 30. While we have focused only on BHXBs, the number of potentially accessible systems doubles if we also consider X-ray binaries with confirmed NSs or candidate BHs (McClintock & Remillard 2006; Casares et al. 2017). Moreover, there are a few hundreds of sources whose detailed properties are still unclear (Liu et al. 2006, 2007). Thus, we expect that TESS photometry will be useful for characterizing at least some X-ray binaries with phase-curve variations. For the self-lensing signal, on the other hand, $N({d}_{\max })$ is ${ \mathcal O }(1)$ or smaller. This number needs to be down-weighted by the eclipse probability, which is typically ${ \mathcal O }(0.1)$, to estimate the number of actual detections. Thus, the detection of self-lensing signals in X-ray binary like systems is not promising for f = 30.

If we set f = 1000 (dotted line in Figure 7), we expect a few hundred systems showing detectable phase-curve signals and a few with self-lensing signals (after taking into account the eclipse probability). These numbers are roughly consistent with the estimates in the corresponding (P, ${M}_{\bullet }$) range from our population models in Section 4 (Figures 5 or 6). This suggests that our population model is roughly consistent with the observed BHXB population if there exist $\approx 30$ times more detached BH–star systems than X-ray systems. The value appears to be compatible with the result of a population synthesis study (Tutukov & Yungelson 2002) in the order of magnitude sense.