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 stateoftheart population synthesis models with efficient postNewtonian evolutions, thus tracking sources from stellar formation to gravitationalwave 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 blackhole 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 thirdgeneration detectors. We find that three formation pathways, differentiated by the order of core collapse and commonenvelope 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 blackhole 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 pairinstability supernovae. Measurements of blackhole spin orientations have enormous potential to constrain specific evolutionary processes in the lives of massive binary stars.
pacs:
I Introduction
Gravitationalwave (GW) observations of merging blackhole (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 stellarorigin 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 commonenvelope 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 manybody 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 starbased BH formation models. Notable exceptions include models where older populations of stars are responsible for presentday BHs Kinugawa et al. (2014); Belczynski et al. (2017a), as well as predictions which make use of largescale 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 stellarbased compactbinary 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 postNewtonian (PN) spinorbit and spinspin 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 stateoftheart 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 onespin dynamics through the effective spin parameter (which is easier to measure; Sec. III) and genuine twospin 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 blackhole evolution
We perform binary star evolutions using the collection of semianalytic 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 mergerrate 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.

“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 datapoints 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 spinprecession 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, populationsynthesis 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 spinorbit 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 models
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 WolfRayet 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 leadingorder 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 secondborn 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 multitimescale 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 spinprecession 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 metallicitydependent 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 signaltonoise 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 thirdgeneration detector in 40km 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, equalmass, nonspinning BH binary with sourceframe 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 ().
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 
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. Thirdgeneration 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 hangup, a wellknow effect in BH binary dynamics which causes binaries with aligned (antialigned) 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 thirdgeneration detectors: future instruments will detect virtually all stellarmass 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 binarystar 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 fieldbinary formation channel corresponds to the first case, “BH1 CE BH2”: the heavier star collapses first and forms the heavier BH, a commonenvelope 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 massratio reversal is known to have potentially strong impact on the spin distribution Kesden et al. (2010a); Gerosa et al. (2013).
The detection rates and their fraction in each channel are shown in Fig. 3 as a function of the kickvelocity 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 massratio 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 commonenvelope 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 parameterestimation 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 
We first illustrate our predictions for the marginalized distributions of . A summary of our findings is provided in Table 2, where we report detectionweighted medians and 90% confidence intervals of for each of our simulations.
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 thirdgeneration 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 hangup 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 highfrequency range will make us sensitive to lower mass systems which, in the collapse model, have preferentially high spins (cf. Fig. 1). The hangup effect is still present, but turns out to be subdominant.
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 effectivespin 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 realigned 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).
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 spinmagnitude assumption (especially for thirdgeneration 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 effectivespin 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 largescale populationsynthesis simulations.
iii.3 Mass dependence
In Fig. 7 we present predictions for the effective spins of BH binaries with different total sourceframe 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 hangup 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 spinmagnitude 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.
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 commonenvelope 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 commonenvelope 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 commonenvelope 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 stateoftheart astrophysical populations) is that BH spin orientations near merger fall into three wellseparated 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
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.

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

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 largemisalignment 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.

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.

The second supernova (SN2) imparts another tilt to the orbital plane. Before this second kick, the binary already underwent a commonenvelope 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 postSN tilts. 
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 spinspin and spinorbit couplings. These are conditioned to keep constant Racine (2008); Gerosa et al. (2015), such that increases only when decreases, and vice versa.
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 spinspin and spinorbit 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 spinmagnitude 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 .
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 stateoftheart 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.
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.

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

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

The configuration is forbidden; the precession cycle consists of librations about (L).
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 radiationreaction 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