Spin orientations of merging black holes formed from the evolution of stellar binaries

Spin orientations of merging black holes formed from the evolution of stellar binaries

Abstract

We study the expected spin misalignments of merging binary black holes formed in isolation by combining state-of-the-art population synthesis models with efficient post-Newtonian evolutions, thus tracking sources from stellar formation to gravitational-wave detection. We present extensive predictions of the properties of sources detectable by both current and future interferometers. We account for the fact that detectors are more sensitive to spinning black-hole binaries with suitable spin orientations and find that this significantly impacts the population of sources detectable by LIGO, while this is not the case for third-generation detectors. We find that three formation pathways, differentiated by the order of core collapse and common-envelope phases, dominate the observed population, and that their relative importance critically depends on the recoils imparted to black holes at birth. Our models suggest that measurements of the “effective spin” parameter will allow for powerful constraints. For instance, we find that the role of spin magnitudes and spin directions in can be largely disentangled, and that the symmetry of the effective spin distribution is a robust indicator of the binary’s formation history. Our predictions for individual spin directions and their precessional morphologies confirm and extend early toy models, while exploring substantially more realistic and broader sets of initial conditions. Our main conclusion is that specific subpopulations of black-hole binaries will exhibit distinctive precessional dynamics: these classes include (but are not limited to) sources where stellar tidal interactions act on sufficiently short timescales, and massive binaries produced in pulsational pair-instability supernovae. Measurements of black-hole spin orientations have enormous potential to constrain specific evolutionary processes in the lives of massive binary stars.

pacs:
1

I Introduction

Gravitational-wave (GW) observations of merging black-hole (BH) binaries have the potential to unveil the fate of massive stars. As they exhaust all the available fuel, stars with initial masses are expected to undergo gravitational collapse. About of them are predicted to form BHs O’Connor and Ott (2011). The detection of stellar-origin BHs in a binary system requires not only the formation of BHs in the first place, but also the occurrence of astrophysical processes that can dissipate enough energy and angular momentum to bring the orbital separation below , where GW damping can drive the binary to merger Peters (1964).

There are two main classes of formation models, depending on whether (i) the two BHs spend their entire lives together as stars, or (ii) they form separately and meet later. In models belonging to class (i), BH binaries are the end product of the life of binaries of massive stars Postnov and Yungelson (2014). Each of the two stars undergoes gravitational collapse and, if the binary is not disrupted, a binary BH is left behind. A common-envelope phase – where the core/remnant of one the two objects sinks into the outer layers of its companion Paczynski (1976) – is typically invoked to dissipate enough angular momentum and produce a merging binary. Models of class (ii) instead require dense stellar environments to facilitate the assembly of multiple BHs and many-body interactions to harden the binary Benacquista and Downing (2013). For a comprehensive review on BH binary formation channels see, e.g., Mandel and Farmer (2018); B. P. Abbott et al. (LIGO and Virgo Scientific Collaboration) (2016a) and references therein.

The most obvious observable to confirm or rule out formation channels is the merger rate, currently constrained to the range B. P. Abbott et al. (LIGO and Virgo Scientific Collaboration) (2017a). These large uncertainties leave ample room for models in both classes to match the observational constraints which, at present, do not allow to confirm or rule out any of the preferred scenarios.

Measurements of the BH masses also tend to be poorly constraining, partly because of a selection bias: more massive systems are visible farther out. This tends to wash out differences in the intrinsic distributions, such that the observable distributions predicted by various models all tend to overlap. In practice, observations could be necessary before strong constraints can be placed using mass measurements Stevenson et al. (2015); Zevin et al. (2017); Kovetz et al. (2017); Barrett et al. (2018) (although sharp features like a mass cutoff will be accessible earlier Fishbach and Holz (2017); Talbot and Thrane (2018); Wysocki et al. (2018a)).

Merger redshifts are also weak observables. They are expected to be set by the star formation history Madau and Dickinson (2014), which is essentially the same in all star-based BH formation models. Notable exceptions include models where older populations of stars are responsible for present-day BHs Kinugawa et al. (2014); Belczynski et al. (2017a), as well as predictions which make use of large-scale cosmological simulations Mapelli et al. (2017); Lamberts et al. (2018).

BH spin magnitudes can be very powerful observables for constraining the physics of individual massive stars, but provide a less effective way to distinguish between stellar-based compact-binary formation channels. Spin magnitudes are expected to be set by stellar collapse dynamics O’Connor and Ott (2011), and should therefore be similar for BHs formed either in galactic fields or dynamically. A possible handle could be provided by dependence on the star’s metallicity, which is expected to impact processes like angular momentum transport and mass loss. Again, these observables might turn out to be particularly useful to constrain specific mechanisms, such as scenarios where previous mergers, rather than stellar collapse, are responsible for forming the merging BHs Gerosa and Berti (2017); Fishbach et al. (2017); Kovetz et al. (2018).

Binary eccentricities may also provide information on some specific models Antonini and Perets (2012); Samsing (2018). However, eccentricities from the most favored scenarios are expected to be too low in the LIGO/Virgo band to provide stringent constraints B. P. Abbott et al. (LIGO and Virgo Scientific Collaboration) (2017b). To this end, LISA observations at low frequencies (when binaries are not yet fully circularized) may turn out to be crucial Nishizawa et al. (2016, 2017); Breivik et al. (2016); Samsing and D’Orazio (2018).

The most promising observables to shine light on BH binary formation are the spin directions. Spins of BHs formed following dynamical encounters are expected to be isotropically distributed. This is because the BH binary evolution is set by the astrophysical environment, whose coupling to the BH spins is known to be negligible Cashen et al. (2017). It is worth pointing out, however, that some angular momentum from the cloud that formed the cluster could be transferred to the stellar spins, thus introducing correlations between their directions Corsaro et al. (2017). Even if present, these correlations are expected to be largely washed out by the many dynamical encounters leading to the formation of the GW sources.

Conversely, the spin directions of BH binaries formed in isolation are greatly influenced by the evolutionary paths of their stellar progenitors. The two stars will form a binary BH without prominent interactions with other bodies, thus “carrying memory” of some of the physical mechanisms occurring during their history. Even in the simplest models where stellar spins are initially aligned to the binary’s orbital angular momentum, misalignments are expected to be introduced by recoil velocities imparted to the BHs at birth. These “supernova kicks” tilt the orbital plane, thus introducing some misalignment between the orbital angular momentum and the spin directions Kalogera (2000). Tidal interactions can also influence the spin directions, generically acting towards realigning spins with the orbital angular momentum Hut (1981). After the BH binary is formed, spin directions are further modified by post-Newtonian (PN) spin-orbit and spin-spin couplings during the long inspiral phase before the binary becomes detectable by LIGO and Virgo Apostolatos et al. (1994). PN effects tend to separate different subpopulations, hence greatly improving model distinguishability Gerosa et al. (2013). The effectiveness of BH spin tilts at constraining formation channels was already explored in previous work through astrophysical models Mandel and O’Shaughnessy (2010); Gerosa et al. (2013); Rodriguez et al. (2016); Antonini et al. (2018); Belczynski et al. (2017b), simulated LIGO/Virgo data Gerosa et al. (2014); Vitale et al. (2014); Trifirò et al. (2016); Stevenson et al. (2017); Vitale et al. (2017a); Talbot and Thrane (2017) and actual GW observations O’Shaughnessy et al. (2017); Farr et al. (2017); Wysocki et al. (2018b); Farr et al. (2018); Wysocki et al. (2018a).

This paper presents a comprehensive study of the expected spin direction distributions of BH binaries formed from isolated pairs of stars. Using the StarTrack Belczynski et al. (2008) and precession Gerosa and Kesden (2016) numerical codes, we combine for the first time state-of-the-art evolutions of binary stars to accurate PN spin tracking and coherently model spin evolution from formation to detection (Sec. II). We present forecasts for both approximate one-spin dynamics through the effective spin parameter (which is easier to measure; Sec. III) and genuine two-spin effects (which encode more information; Sec. IV). We then illustrate predictions of our models in terms of the spin morphologies identified in Kesden et al. (2015); Gerosa et al. (2015) (Sec. V). We conclude with prospects for constraining these mechanisms with current and future GW detectors (Sec. VI). Unless otherwise noted, we use geometrical units ().

Our database is publicly available at github.com/dgerosa/spops, where we also provide a convenient python module (called spops) to facilitate its exploration.

Ii Methods: stellar and black-hole evolution

We perform binary star evolutions using the collection of semi-analytic prescriptions implemented in the StarTrack code Belczynski et al. (2002, 2008); Dominik et al. (2012); Fryer et al. (2012); Dominik et al. (2013, 2015); Belczynski et al. (2016a, b). Each evolution results in a BH binary characterized by masses (with ; or alternatively and ), spin magnitudes (with ), directions , and merger-rate weight (see Sec. II.3 below). The direction of each spin is described by a polar angles (relative to the direction of the orbital angular momentum) and by an azimuthal angle in the orbital plane, .

Our suite of models is described below. Each model has a single free parameter (setting the magnitude of the kick velocities) and three flags corresponding to our assumptions on spin magnitudes, tidal interactions and the sensitivity of GW detectors. With the exception of BH kicks and spins, all other assumptions are the same as in model M10 of Belczynski et al. (2016a). We refer the reader to that paper for a more comprehensive description of our population synthesis simulations.

ii.1 Spin magnitudes

As for the BH spin magnitudes, we implement three different models:

  • uniform”: We assume the dimensionless BH Kerr parameters to be uniformly distributed in , independently of the other binary parameters.

    Figure 1: Model collapse for the BH spin magnitude as a function of the BH mass. In this model, heavier (lighter) collapsing stars preferentially form BHs with smaller (larger) spins. Filled circles shows data-points from the simulations reported by Belczynski et al. (2017b) at various metallicities , while empty triangles show our resampled distribution. Dashed (solid) lines illustrate our construction procedure (see text). The hard cutoff at (dotted line) is due to pulsational pair-instaibility supernovae as implemented in StarTrack Belczynski et al. (2016a).
  • collapse”: Simulations of stellar collapse show that stars with large (low) mass tend to form slowly (highly) rotating BHs Belczynski et al. (2017b); Qin et al. (2018). This feature introduces a specific correlation between masses and spins, with potentially critical impact on the predicted GW sources. Here we implement a very simple prescription to qualitatively capture this effect, leaving more robust explorations to future work. We use evolutionary simulations of stars with specific angular momentum transport from Eggenberger et al. (2012); Georgy et al. (2013) as reported in Table 3 of Belczynski et al. (2017b), together with the approximate expression

    (1)

    obtained from Fig. 1 of Belczynski et al. (2017b). Since there are not enough data-points to construct meaningful interpolants, we opt for the following heuristic approach. At low (large) masses, spins appear to be centered about () with a scatter of (). The turnover between the two regimes is at with a scatter of about . Our procedure is illustrated in Fig. 1. We first construct two curves

    (2)

    where (the upper/lower signs refer to the upper/lower limits in Fig. 1). Spins are then generated with random samples uniform in the region in between the two curves. We argue that this model captures some of the key features found in Belczynski et al. (2017b), namely that larger BH masses tend to correlate with smaller spins, while at the same time reflecting the large uncertainties of those results.2

  • max”: In order to maximize the effects of spin precession dynamics and highlight some trends, we also run a set of models where all BHs are maximally spinning ().

We note that the first two models implement rather conservative assumptions regarding the expected spin-precession dynamics (cf. Gerosa et al. (2013); Stevenson et al. (2017), where only very high spins are considered). Our approach complements that of Belczynski et al. (2017b), where a specific model (their Eq. 3) was assumed for the BH spin magnitude. More work is needed to fully include the impact of the metallicity on the expected BH spin magnitudes, which is here neglected.

ii.2 Spin directions

We assume stellar spins to be initially aligned to the orbital angular momentum of the binary (). This is a conservative assumptions made by most, if not all, population-synthesis models (but see e.g. Albrecht et al. (2012); Postnov and Kuranov (2017)). As the first star collapses and forms a BH, the resulting kick tilts the orbital plane Kalogera (2000); Hurley et al. (2002) and introduces a spin-orbit misalignment (). It is worth pointing out that spin misalignments are induced by asymmetric mass and neutrino emission during core collapse (“natal kicks”), while symmetric mass loss only impacts the binary’s center of mass (“Blaauw kicks” Blaauw (1961)). Supernova kicks are drawn from a Maxwellian distribution with 1D dispersion , independently of the mass of the system. We generate 7 different models3 at and 265 km/s. The largest value km/s corresponds to the observational constraints from pulsar proper-motion measurements Hobbs et al. (2005). We adopt this approach because it constitutes a simple and well-defined one-parameter family of models to illustrate the main trends of the BH spin alignment distributions. More elaborate (and perhaps more physical) prescriptions where, e.g., kicks are suppressed by fallback material Fryer (1999) can be constructed by appropriate mass-dependent mixture of our distributions Wysocki et al. (2018b).

After the first kick, tidal interactions may realign one of the spins. In between the two supernova explosions, the system is formed by a BH and a (perhaps evolved) star Kushnir et al. (2016). Since tidal interactions scale with the cube of the size of the object, tides raised on the star by the BH are much more effective than tides raised on the BH by the star.

In the spirit of introducing only minimal assumptions, we implement three prescriptions Gerosa et al. (2013); Belczynski et al. (2017b) and postpone a more careful treatment of tidal spin alignment to future work Steinle et al. ().

  • alltides”: Tidal interactions align all stellar spins in between the two explosions. This corresponds to setting either , or , , depending of which star explodes first.

  • notides”: None of the spins is realigned.

  • time”: We attempt a physical model for tidal interactions following Kushnir et al. (2016); Zaldarriaga et al. (2018). In particular they estimate the tidal alignment time to be

    (3)

    where and are the masses of the star and the BH, respectively, and is the binary separation. We compute Eq. (3) from the StarTrack data before the second supernova, and compare the result with both the time between the two explosions and the typical lifetime of a Wolf-Rayet star yr Kushnir et al. (2016); Zaldarriaga et al. (2018). The star’s spin is realigned if

    (4)

The azimuthal angles may also evolve in between the two explosions because of relativistic spin precession. We compare the time between the two explosions to the leading-order precession timescale Gerosa et al. (2015)

(5)

and set if or draw randomly if .

It is worth noting that tidal interactions (here considered only regarding the spin directions) are expected to affect the spin magnitude of the second-born BH as well Kushnir et al. (2016); Zaldarriaga et al. (2018); Qin et al. (2018). Tidally locked stars are going to be both aligned and spun up. This behavior is partially captured by the combination of the collapse and time models, as lower mass stars are both assigned high spins and lower . A more systematic study of the effect of tides in BH binaries formation pathways is under development Steinle et al. (). The potential impact of mass transfer on both spin magnitudes and directions is also an avenue of future improvement.

At the second explosion, another supernova kick will further tilt the orbital plane. This finally results in the angles and at BH binary formation. Each system now needs to be evolved from the separation where it forms down the LIGO/Virgo band ( Hz). We use the numerical code precession Gerosa and Kesden (2016), which implements multi-timescale methods to efficiently evolve BH systems over their very long inspirals before merger Gerosa et al. (2015); Kesden et al. (2015). This is a crucial improvement over previous work Gerosa et al. (2013); Stevenson et al. (2017) (but see Rodriguez et al. (2016); Wysocki et al. (2018b), which also use the same code), since it was shown that integrations from separations as large as may be necessary to fully capture spin-precession effects Gerosa et al. (2015).

ii.3 Detectability

For each set of assumptions, our procedure generates a sample of weighted BH binaries Dominik et al. (2013, 2015); Belczynski et al. (2016b). These weights correspond to the merger rate contribution of each evolutionary track, taking into account redshift- and metallicity-dependent star formation history Dominik et al. (2013), as well as the antenna pattern of the GW interferometer Dominik et al. (2015).

Contrary to previous StarTrack studies, we now take into account spin corrections in the calculation of the merger weight. As illustrated below, this is a crucial point to faithfully predict spin distributions which, if neglected, could lead to sizable biases Ng et al. (2018); Baird et al. (2013); Dominik et al. (2015); O’Shaughnessy et al. (2010); Reisswig et al. (2009). We generate GW signals using the IMRPhenomPv2 Hannam et al. (2014) waveform model as implemented in the pyCBC pipeline Usman et al. (2016). We compute signal-to-noise ratios (SNRs) using three different noise curves:

  • LIGO: the expected sensitivity for Advanced LIGO in its design configuration B. P. Abbott et al. (LIGO and Virgo Scientific Collaboration) (2018);

  • Voyager: a planned upgrade designed to maximize the science return within the current LIGO facilities LIGO Scientific Collaboration ();

  • Cosmic Explorer: a proposed third-generation detector in 40-km scale facilities B. P. Abbott et al. (LIGO and Virgo Scientific Collaboration) (2017c).

For simplicity we consider single detectors with an SNR threshold of 8, which is a reasonable approximation to mimic realistic data analysis procedures J. Abadie et al. (LIGO and Virgo Scientific Collaboration) (2010). For reference, an optimally located and oriented, equal-mass, non-spinning BH binary with source-frame total mass of will have an SNR smaller than 8 at redshifts for Advanced LIGO, for Voyager and for Cosmic Explorer. Detection rates (in units of ) are then computed as detailed in Belczynski et al. (2016b) (see also Finn and Chernoff (1993); Finn (1996); Dominik et al. (2015)) using the public code gwdet Gerosa ().

Figure 2: Detection rates for LIGO (solid line), Voyager (dashed line) and Cosmic Explorer (dotted line) for models with different kick speed parameters , assuming the time tidal model and the collapse spin model. All other models give qualitatively similar curves; the max spin model yields marginally higher rates for LIGO (cf. Table 1).
   Detector Spins Tides Natal kick
km/s km/s km/s km/s km/s km/s km/s
   LIGO collapse time
   LIGO collapse alltides
   LIGO collapse notides
   LIGO uniform time
   LIGO uniform alltides
   LIGO uniform notides
   LIGO max time
   LIGO max alltides
   LIGO max notides
   Voyager collapse time
   Voyager collapse alltides
   Voyager collapse notides
   Voyager uniform time
   Voyager uniform alltides
   Voyager uniform notides
   Voyager max time
   Voyager max alltides
   Voyager max notides
   3 gen. collapse time
   3 gen. collapse alltides
   3 gen. collapse notides
   3 gen. uniform time
   3 gen. uniform alltides
   3 gen. uniform notides
   3 gen. max time
   3 gen. max alltides
   3 gen. max notides
Table 1: Detection rates in units of yr for all our simulations. Results for LIGO assumes the expected design sensitivity of the instrument Usman et al. (2016). Voyager is a planned instrumental upgrade to be located in the current LIGO facilities LIGO Scientific Collaboration (). Here we use Cosmic Explorer B. P. Abbott et al. (LIGO and Virgo Scientific Collaboration) (2017c) as an illustrative example of what would be possible with third-generation detectors.

The detection rate is a steep function of the natal kick velocity Belczynski et al. (2002). For large values of , more and more stellar progenitor binaries are unbound by natal kicks and fail to form GW sources. This is shown in Fig. 2 for a subset of our models; we find that drops by about a factor between km/s and km/s. Third-generation detectors increase the expected detection rates by a factor of compared to LIGO at design sensitivity B. P. Abbott et al. (LIGO and Virgo Scientific Collaboration) (2017c).

Detection rates for each of our model variations are reported in Table 1. For LIGO, the max spin models predicts higher rates (by about , i.e a factor 1.5) when compared to models with uniform spin distributions. This behavior is due to the orbital hang-up, a well-know effect in BH binary dynamics which causes binaries with aligned (anti-aligned) spins to have a larger (smaller) horizon distance Campanelli et al. (2006); Scheel et al. (2015). For distributions with spins mostly aligned like ours, binaries with large spin magnitudes are therefore easier to detect. This result refines the rough estimate of Dominik et al. (2015), where spins were estimated to increase rates by at most a factor of 3 (cf. also Reisswig et al. (2009)) Interestingly, this rate increase disappears for third-generation detectors: future instruments will detect virtually all stellar-mass BH mergers in the Universe, irrespectively of their spins. For the same reason, models with more misaligned spins (notides) have marginally lower rates than models where all spins are realigned by tidal interactions (alltides).

ii.4 Simplified pathway classifications

StarTrack provides full information on the various processes and the stellar types involved during each phase of the binary-star evolution. For this paper, we found particularly illustrative to simplify the classification of the evolutionary pathways marking the formation of the heavier BH (label “BH1”), the formation of the lighter BH (label “BH2”) and the occurrence of common envelope phases (label “CE”). All of our stellar evolutions can be classified into eight mutually exclusive channels:

1.BH1 CE BH2 5.CE BH1 CE BH2
2.BH2 CE BH1 6.CE BH2 CE BH1
3.CE BH1 BH2 7.BH1 BH2
4.CE BH2 BH1 8.BH2 BH1

The vanilla field-binary formation channel corresponds to the first case, “BH1 CE BH2”: the heavier star collapses first and forms the heavier BH, a common-envelope phase tightens the binary, and finally the companion star forms the lighter BH. The second channel, “BH2 CE BH1”, corresponds to cases where the light BH formed first: such mass-ratio reversal is known to have potentially strong impact on the spin distribution Kesden et al. (2010a); Gerosa et al. (2013).

Figure 3: Detection rates (bottom panel) and normalized rate fractions (top panels) of BH binaries formed via different channels as a function of the natal kicks . In this notation “BH1” and “BH2” stand for the formation of the heavier and lighter BH, respectively, while “CE” stand for the occurrence of a common-envelope phase (cf. Sec. II.4). Results are shown for the time and collapse spin model and weighted by LIGO detection rates. Results for Cosmic Explorer and other spin models are qualitatively similar.

The detection rates and their fraction in each channel are shown in Fig. 3 as a function of the kick-velocity dispersion parameter  Hobbs et al. (2005). For small kicks, most binaries follow the standard picture and evolve through a common envelope phase between the two stellar collapses. At we have . Two thirds of these binaries follow the more standard pathway where the large BH is formed first, while the rest undergo mass-ratio reversal.

If kicks are larger, the majority of binaries are found in the “CE BH1 BH2” channel. In this regime, systems are typically unbound by the first explosion (causing a drop in the rates), unless a common-envelope phase takes place before the explosion. Common envelope shrinks the orbital separation by orders of magnitude, thus dramatically increasing the chance of the binary surviving the first natal kick. In particular, in the extreme case km/s we obtain .

Channels with zero or two common envelope phases are always subdominant, and represent at most of the population.

Iii Results: effective spin

The spin parameter which is currently best measured in GWs B. P. Abbott et al. (LIGO and Virgo Scientific Collaboration) (2016b) is the effective spin Damour (2001); Racine (2008); Kesden et al. (2010a)

(6)

(this is equivalent to in the notation of Gerosa and Kesden (2016); Gerosa et al. (2015); Kesden et al. (2015); Gerosa et al. (2017); Zhao et al. (2017)). The effective spin is a constant of motion at 2PN order Racine (2008); Kesden et al. (2010a), and is therefore an excellent parameter to parameterize the dynamics because it depends very weakly on the frequency/time at which it is measured by parameter-estimation algorithms.

As evident from the definition (6), measurements of are inevitably plagued by a degeneracy between the spin magnitudes and their directions: small (large) values of could be realized by either small (large) spin magnitudes or large (small) misalignment angles. For strategies to maximize the astrophysics that can be inferred from measurements of alone, see e.g. Farr et al. (2017, 2018); Ng et al. (2018); Farr et al. (2018).

iii.1 Marginalized distributions

   Detector Spins Tides Natal kick
km/s km/s km/s km/s km/s km/s km/s
   LIGO collapse time
   LIGO collapse alltides
   LIGO collapse notides
   LIGO uniform time
   LIGO uniform alltides
   LIGO uniform notides
   LIGO max time
   LIGO max alltides
   LIGO max notides
   Voyager collapse time
   Voyager collapse alltides
   Voyager collapse notides
   Voyager uniform time
   Voyager uniform alltides
   Voyager uniform notides
   Voyager max time
   Voyager max alltides
   Voyager max notides
   3 gen. collapse time
   3 gen. collapse alltides
   3 gen. collapse notides
   3 gen. uniform time
   3 gen. uniform alltides
   3 gen. uniform notides
   3 gen. max time
   3 gen. max alltides
   3 gen. max notides
Table 2: Medians of the marginalized effective-spin distributions in each of our model variations. Errors refer to the 5 and 95 confidence levels, respectively. All percentiles reported in this table have been weighted with the corresponding detection rates.

We first illustrate our predictions for the marginalized distributions of . A summary of our findings is provided in Table 2, where we report detection-weighted medians and 90% confidence intervals of for each of our simulations.

Figure 4: Marginalized effective-spin distribution predicted by our three spin magnitude modes (collapse,max,uniform) as detectable by LIGO and Cosmic Explorer. Results are shown for km/s and the time model of tidal interactions.

Among our model variations, the spin magnitudes have the largest impact on . Fig. 4 shows the distributions of predicted by our three different spin models, uniform, collapse and max, assuming km/s and the time model of tidal interactions. The max model predicts a sharp peak at , while the uniform model presents a broader peak at . The collapse case acts like a rough mixture of the two, with light BHs presenting preferentially large spins , while the spins of heavier BHs can span a wider range.

An interesting feature of these distributions is their dependence on the detector sensitivity. For the max model, switching to a third-generation detector decreases the typical effective spin: the median in goes from to (assuming the time, km/s distribution as shown in Fig. 4). The same is true for the uniform model, where the median decreases from for LIGO to for Cosmic Explorer. The orbital hang-up effect (cf. Sec. II.3) causes a selection bias on GW measurements: binaries with negative (positive) have a shorter (longer) waveform and therefore are harder (easier) to detect Ng et al. (2018). A detector with better sensitivity reduces this selection bias, thus pushing the median of the detectable events to lower values.

The collapse model behaves in the opposite way: better instruments will detect larger ’s (median increasing from to ). This is because instrumental improvements in the high-frequency range will make us sensitive to lower mass systems which, in the collapse model, have preferentially high spins (cf. Fig. 1). The hang-up effect is still present, but turns out to be subdominant.

Figure 5: Median values (dashed lines) and 90% confidence intervals (solid lines and shaded areas) for as a function of natal kicks and tidal interactions (time, alltides, notides). Results are shown for uniform spin magnitudes and weighted with LIGO detection rates.

Fig. 5 illustrates the effect of tidal interactions on the predicted values of . In our models tides only affect the spin orientations and, as expected, produce larger values. Within the context of these models, tides are less important than natal spins to predict the effective-spin distributions. Notably, tides mainly affect the small- tail of the population. As illustrated at length below, negative values of are hard to explain in the alltides model (where all stellar spins are re-aligned in between the two explosions), while they are relatively easy to accommodate with both the notides and the time models.

iii.2 On the sign and symmetry of

From Eq. (6), it is obvious that only largely misaligned spins can produce negative values of . It has been suggested that a single confirmed measurement of a system with could rule out isolated BH formation for that event in favor of dynamical interactions Rodriguez et al. (2016). One of the events observed so far (GW151226 B. P. Abbott et al. (LIGO and Virgo Scientific Collaboration) (2016c)) has at very high confidence. Some of the other events present more posterior weight at negative values, but cannot be ruled out. We stress that these significance assessments have to be taken with care as they depend on the Bayesian prior used in the analysis Vitale et al. (2017b).

Figure 6: Fraction of binaries with negative effective spin as a function of natal kicks (x-axis), tidal interaction (colors) and spin-magnitude model (linestyles). Top (bottom) panels shows results for LIGO (Cosmic Explorer).

Fig. 6 shows the predicted rate fraction of BH binaries with in each of our models. As expected, misalignments are larger for larger kicks and, consequently, increases as increases. The typical fraction of binaries with negative effective spins detectable by LIGO ranges from to (with the exception of km/s, where by construction).

Our results shows that isolated pairs of stars can explain single events with , in disagreement with the main claim made by Rodriguez et al. (2016) (but see their Fig. 3). Obviously, since we are assuming that stars are initially aligned with the orbital angular momentum, BH spins cannot be misaligned if kicks are not present. The fiducial model of Rodriguez et al. (2016) heavily suppresses kicks for BHs compared to neutron stars, thus effectively preventing misalignments. A more conservative statement is the following: single detections with would point towards dynamical interaction, if stellar spins are initially aligned and BH kicks are heavily suppressed. Even moderate kicks of km/s allow for .

Together with kicks, tidal interactions are important to determine the sign of . Higher (lower) fractions of negative effective spins are predicted for the notides (alltides) model, while the time model lies somewhere in between. Notably, is largely independent of the spin-magnitude assumption (especially for third-generation detectors).

This suggests that, at least in the context of well specified astrophysical models like ours, measurements alone can partially break the degeneracy between spin magnitudes and spin directions encoded in Eq. (6). Natal spins mainly determine the broad shape of the distribution (Fig. 4), while alignment processes have a clean impact on the low- tail (Fig. 6).

Although negative values of are possible, our distributions are far from being symmetric (cf. e.g. Table 2 where all medians are ). On the contrary, dynamical formation channels predicts spins isotropically distributed (although see Corsaro et al. (2017)), which corresponds to a marginalized effective-spin distribution symmetric about . Our models suggest that the symmetry of the is a robust indicator to distinguish isolated binary formation from dynamical interactions. We therefore confirm the ideas put forward by Farr et al. (2017, 2018) with large-scale population-synthesis simulations.

Figure 7: Effective spins and detection rates in bins of total mass . We show results for our three spin models (left: collapse; middle: max; right: uniform) assuming km/s, the time model for tides, and the LIGO sensitivity curve. Thick solid (dashed) lines show median (90% confidence interval) of , as reported on the right y-axis. Light histograms show the cumulative detection rates in each mass bin, as reported on the right y-axis.
Figure 8: Medians of as a function of the total mass for all of our model variations. The top panel shows results weighted by LIGO detection rates, while the bottom panel assumes a third-generation detector (Cosmic Explorer). Colors differentiate our three spin-magnitude models (collapse: blue; max: orange; uniform: green). In each series, the various lines are obtained by varying over tidal interactions (time, alltides, notides) and kick magnitudes ( km/s).

iii.3 Mass dependence

In Fig. 7 we present predictions for the effective spins of BH binaries with different total source-frame masses. Results for the collapse model directly reflect the injected relationship between BH masses and spins. In the absence of this correlation, a simpler trend emerges, namely that kicks more easily misalign light systems. This is especially evident in the max case because all BH spin magnitudes are equal.

Fig. 7 also shows the expected detection rates as a function of . The three panels are constructed with the very same StarTrack evolution ( km/s), which predicts a mass spectrum peaking at about . This is the strongest feature visible in all three distributions shown in Fig. 7.

Differences in between the three panels are a direct consequence of the spinning waveform model used to compute the horizon distance. This effect was mostly neglected in previous StarTrack studies, with the exception of Dominik et al. (2015). For partially aligned systems like ours, the orbital hang-up effect facilitates the detection of BH binaries with large spin magnitudes. The max rates are therefore higher than those predicted by the uniform model. In model collapse, where heavy BHs spin slower compared to less massive ones, detection rates at large (low) are suppressed (enhanced). As expected, the behavior changes around , which is roughly twice the value of the turnover of Fig. 1.

Fig. 8 shows medians of as a function of for all kick, spin and tide variations. The main takeaway here is that distributions are qualitatively very similar for all models of tidal interactions and natal kicks, and only depend on the spin-magnitude variation. This finding further stresses one of the points made above: measurements alone can provide powerful constraints on BH natal spins, even in the presence of misalignment processes.

Figure 9: Detection rates in bins of effective spin divided into formation channels. Here “BH1” and “BH2” stand for the formation of the heavier and lighter BH, respectively, while “CE” stand for the occurrence of a common-envelope phase (cf. Sec. II.4). Results are shown for the notides, uniform spin model as detectable by LIGO. Other models have qualitatively similar results. Natal kicks are varied in the three panels as indicated in the legend.

iii.4 Constraints on formation channels

Finally, we present our predictions for in the eight different formation channels introduced in Sec. II.4. Fig. 9 shows results for some of our uniform, notides models. It is illustrative to look at this variation in particular because, as described above (cf. Figs. 4 and 5), it maximizes the fraction of binaries with far from unity. We stress, however, that the trends described here are illustrative of all of our distributions.

As already shown in Fig. 3, the fraction of “standard” binaries with a common-envelope evolution in between the two supernovae decreases with the kick velocity dispersion parameter . For km/s kicks unbind most binaries at the first supernova, unless the binary separation was already tight because of an earlier common-envelope phase. Binaries in those channels (see in particular “CE BH1 BH2” in Fig. 3) are largely unaffected by kicks. Their orbital angular velocity is so large that they not only remain bound, but also roughly aligned.

Largely misaligned binaries all belong to the more standard “BH1 CE BH2” and “BH2 CE BH1” channels, independently of . This result illustrates a clean prediction of our models: binaries with small are formed following very specific pathways, namely those which present a common-envelope phase between the formation of the two BHs (and not earlier).

This observation can be rephrased as follow: if kicks are large, binaries in the “BH1 CE BH2” and “BH2 CE BH1” channels are either unbound or, if they survive, they are largely misaligned. This behavior can be explained with some simple kinematics. For a circular orbit, the spin misalignment angle imparted by a kick is (Kalogera, 2000; O’Shaughnessy et al., 2017)

(7)

where is the kick velocity, is the orbital velocity and is the direction of the orbital angular momentum before the explosion. Since is drawn from a Maxwellian distribution, the component and are Gaussian. In the limit of large , this implies that is uniformly distributed (Marsaglia, 1972). As the kick increases, more binaries are unbound while the distribution of misalignments flattens.

Iv Results: spin directions

We now explore predictions of our models for the individual directions of the two spins. As already outlined in Sec. II, the mutual orientations of the two spins and the orbital angular momentum can be described by three variables: and are the angles between the two spins and the orbital angular momentum, and is the angle between the projections of the two spins onto the orbital plane (see e.g. Fig. 1 of Gerosa et al. (2015) for a schematic representation). The angles and are polar angles, defined in the range , while is an azimuthal angle defined in the range . However, precession cycles are symmetric in PN dynamics Gerosa et al. (2015), so that we can consider without loss of generality.

The punchline of this Section (which generalizes the toy model of Gerosa et al. (2013) to state-of-the-art astrophysical populations) is that BH spin orientations near merger fall into three well-separated subpopulations. These classes of BH binaries carry the imprint of specific physical processes driving the evolution of their stellar progenitors.

iv.1 Evolution of the spin tilts

Figure 10: Evolution of the spin orientations along the lives of BH-binary progenitors detectable by LIGO. The top (bottom) subpanel in each plot shows the tilt () of the object forming the more (less) massive BH. All binaries are aligned before the first supernova (SN1), which imparts a first tilt to both spins. Tidal interactions can realign one of the spins in between the two explosions. The second kick (SN2) sets the spin misalignment angles at BH binary formation. These orientations then evolve under the influence of relativistic spin-spin and spin orbit couplings until they become detectable in GWs (roughly at Hz). At each stage, the median of the distribution is marked with a red line; the blue boxes (bars) include 50% (90%) of the detection rate. Thin gray lines show individual evolutionary tracks for the 100 binaries with the highest detection rates in each sample.

First, we illustrate the evolution of the spin angles during the various steps of binary stellar evolution. There are five key stages where spin directions can change. These are listed below and illustrated in Fig. 10, where we track changes of the two tilt angles along each stage by separating progenitors that form the more () and less () massive BHs.

  1. Our initial assumption is that primordial misalignments are negligible, i.e. at the beginning of each evolution.

  2. The first stage where spin tilts can change is the supernova that forms the first BH (SN1). This is typically, but not always, the collapse event where the more massive BH is formed (c.f Sec. II.4). It turns out that the kick imparted at the first explosion is the dominant effect setting the spin directions in the entire evolution, and all other stages play a subdominant role (cf. O’Shaughnessy et al. (2017), where this consideration was used to estimate from GW151226 data). On average, larger kicks introduce larger misalignments. However, this trend is mitigated by the fact that larger kicks also unbind binaries. Only the harder binaries in the sample survive strong kicks, and those same binaries are harder to tilt. At this stage, the median in the angles (both members receive the same tilt) is , and this number changes only weakly among our kick and spin models. The large-misalignment tail of the tilt distributions, on the other hand, depends strongly on : tilts require km/s. This is consistent with the results already presented in Fig. 6, where indeed curves steepen at about km/s.

  3. After the first explosion, the system is formed by a BH and a star. At this stage, tidal interactions can realign the stellar spin. All stars are realigned in the alltides case where, consequently, one of the two spin misalignment angles drops to zero. In the majority of the cases, tides enforce between the two explosions, because the first explosion typically forms the most massive BH. However, if the less massive BH is formed first, tidal alignment enforces . Both spin misalignments are unchanged in the notides models, where tidal realignment is assumed to be completely inefficient.

  4. The second supernova (SN2) imparts another tilt to the orbital plane. Before this second kick, the binary already underwent a common-envelope phase which greatly tightened the separation.4 Because of their larger orbital velocities, binaries are much harder to tilt at this stage compared to the first explosion. This second kick is virtually irrelevant for km/s. In the case of larger kicks where the “CE BH1 BH2” channel dominates, on average the second explosion increases the misalignment angle (this is trivially true for the spin that was previously realigned by tides). After the second kick, the tilt angles and are in general not equal to each other. Even in the notides cases where before the SN, spin precession might affect the azimuthal angles ( if , c.f. Sec. II) and consequently the post-SN tilts.

  5. Finally, PN evolutions Gerosa and Kesden (2016) are used to propagate binaries from BH formation (after SN2) to detection, here assumed to happen when the GW emission frequency drops below Hz. PN evolutions could last Gyrs, where binaries undergo many precession cycles. The tilt angles are modified by relativistic spin-spin and spin-orbit couplings. These are conditioned to keep constant Racine (2008); Gerosa et al. (2015), such that increases only when decreases, and vice versa.

    Figure 11: Marginalized distributions of the spin angles , and at Hz for binaries detectable by LIGO. We assume natal kicks of km/s and a variety of assumptions on tides and spin magnitudes at collapse. Distributions of and are peaked at 0, with widths . Distributions of at detection carry clear imprints of the underlying assumptions on tidal interactions and realignment.

We now explore in detail the spin orientations at the last stage, when binaries become detectable in GWs.

iv.2 Spin angles at detection

Figure 11 shows marginalized distributions of , and at Hz for a subset of our models. The tilt angle distributions peak at and , which is the initial assumption for all stars in our models. Misalignments are introduced by natal kicks. The typical widths of the and distributions varies from for km/s to for km/s. These curves are largely independent of the chosen spin model.

On the contrary, the behavior of the angle strongly depends on the spin variation. If all binaries are subject to tidal realignment (alltides), the distribution of at 20 Hz is strongly peaked at , while a less prominent peak at is present if tides are suppressed (notides). If only some of the binaries are realigned (time), three distinct subpopulations are present, which pile up at and . Assumptions on the spin magnitude also play a visible role in Fig. 11. PN spin-spin and spin-orbit couplings are weaker for lower spins and the peaks in are consequently less pronounced. It is worth noting, however, that the impact of tides can be clearly seen even for the uniform spin-magnitude model.

The reason for this peculiar behavior of the precessional phase lies in the PN evolution that binaries undergo between formation and detection. Spin precession naturally separates populations that formed with different tilt angles , into different distributions for . This is a well known PN effect, first discovered by Schnittman Schnittman (2004) and later explored in detail by Kesden et al. (2010a, b); Berti et al. (2012); Gerosa et al. (2013, 2015) (see also Afle et al. (2018); Zhao et al. (2017); Gerosa et al. (2017); Correia (2016); Gupta and Gopakumar (2014) for later investigations). Tidal interactions affect whether the tilt angle of the second formed BH is set to zero between the two explosions, thus strongly impacting its tilt angle at formation. As the binary inspirals towards merger, PN spin evolution tend to mix up the and distributions and separate the subpopulations in the variable .

Figure 12: Spin angles for three sub-channels describing whether tidal realignment was efficient or not (“tides on” vs. “tides off”) and whichever BH formed first (“BH1 BH2” vs “BH2 BH1”). The left panel shows a two-dimensional histogram of the tilt angles and at BH formation, where the three subpopulations are clearly separated. Those same binaries are evolved to detection ( Hz) and the variable at that point is shown in the right panel. PN evolution naturally clusters binaries from different sub-channels to “orthogonal” regions in the parameter space. This figure was generated assuming the LIGO sensitivity, km/s, the time tidal alignment model and the max spin model.

To better illustrate and quantify this behavior, let us divide our binaries in three subpopulations. First, we select binaries that were not realigned by tidal interactions between the two explosions (“tides off”). This will correspond to (0%) of the binaries in the notides (alltides) models, and to some other fractions in the time cases. For the remaining binaries that do undergo tidal realignment (“tides on”), we track whether BH1 forms before/after BH2 (thus grouping the eight channels of Sec. II.4 into two). This results in three mutually exclusive subchannels:

1.BH1 BH2 tides on
2.BH2 BH1 tides on
3.tides off

Fig. 12 shows the resulting distribution of the spin angles. Crucially, we pair distributions of and at BH formation to the distribution of at detection. Tides and the order of BH formation strongly separate the tilt angle distributions: if tidal realignment is prevented, while () if tides are present and the realigned star ends up forming the primary (secondary) BH. Because of the long PN evolution before merger Kesden et al. (2015); Gerosa et al. (2015); Gerosa and Kesden (2016), these three populations are found with preferential values of as they enter the LIGO band: for “BH2 BH1 tides on”, for “tides off” and for “BH1 BH2 tides on”.

These findings confirm the toy model developed by some of the present authors Gerosa et al. (2013), which only considered a few fiducial sources (compare e.g. their Fig. 1 to Fig. 12 in this paper). Despite a substantial extension to a much larger population, including state-of-the-art initial conditions provided by StarTrack, and more complex models for tidal alignment, this simple approach provides an accurate description of several key features of the population of spinning binaries.

Motivated by the identification of these three, coarsely identified classes, in the next Section we employ another tool developed in Kesden et al. (2015); Gerosa et al. (2015) to characterize precessing sources: their spin morphology.

V Results: spin morphologies

The spin morphology is a better tool to quantify spin precession in merging BH binaries. It was first introduced by Kesden et al. (2015); Gerosa et al. (2015) (see also Gerosa (2018) for a concise introduction). Here we briefly review the main concepts behind spin morphology, and then explore the implications for our populations of detectable BH binaries.

v.1 A slowly evolving feature

The qualitative shape of the precession cones of the two spins and the orbital angular momentum can be classified into three mutually exclusive classes based on the evolution of the precessional phase . Fig. 13 shows the evolution of during single precession cycles at Hz for some indicative binaries from our distributions. For simplicity, we define a precession cycle to start and end at configurations where the three vectors , and are coplanar. These corresponds to either or . There are therefore three discrete possibilities.

  1. Both configurations and are allowed, and the angle circulates in the full range during each precession cycle (C).

  2. The configuration is forbidden; the precession cycle consists of librations about (L).

  3. The configuration is forbidden; the precession cycle consists of librations about (L).

Figure 13: Evolution of the angle during (half of) a precession cycle at Hz for a sample of detectable BH binaries. We select some binaries among those with higher detection rates from our model with km/s, the time tidal alignment model and the max spin magnitude model (as in Fig. 12). Each line is shaded according to the LIGO detection rate of the corresponding source.

This classification elucidates the results already presented in Figs. 11 and 12. Binaries in the two librating morphologies L and L spend more time close to the coplanar configurations, and are thus more likely to be found with either or , respectively. Binaries in the circulating morphology C behave in the opposite way. For these binaries, the “azimuthal velocity” is larger at and lower at . Sources naturally spend more time where is lower, and are thus more likely to be found with .

The most notable feature about the spin morphology is its slow variation. In BH binary systems, spins vary on both the short precession timescale Apostolatos et al. (1994) and the longer radiation-reaction timescale Peters (1964). The individual spin angles , and all vary on (the same is true for other quantities typically used to parametrize spin precession, like