Tidal Spin-up of Black Hole Progenitor Stars

Gravitational wave observations indicate the existence of merging black holes (BHs) with high spin ($a\gtrsim0.3$), whose formation pathways are still an open question. A possible way to form those binaries is through the tidal spin-up of a Wolf-Rayet (WR) star by its BH companion. In this work, we investigate this scenario by directly calculating the tidal excitation of oscillation modes in WR star models, determining the tidal spin-up rate, and integrating the coupled spin-orbit evolution for WR-BH binaries. We find that for short-period orbits and massive WR stars, the tidal interaction is mostly contributed by standing gravity modes, in contrast to Zahn's model of travelling waves which is frequently assumed in the literature. The standing modes are less efficiently damped than traveling waves, meaning that prior estimates of tidal spin-up may be overestimated. We show that tidal synchronization is rarely reached in WR-BH binaries, and the resulting BH spins have $a \lesssim 0.4$ for all but the shortest period ($P_{\rm orb} \! \lesssim 0.5 \, {\rm d}$) binaries. Tidal spin-up in lower-mass systems is more efficient, providing an anti-correlation between the mass and spin of the BHs, which could be tested in future gravitational wave data. Nonlinear damping processes are poorly understood but may allow for more efficient tidal spin-up. We also discuss a new class of gravito-thermal modes that appear in our calculations.


INTRODUCTION
The spins of stellar-mass black holes (BHs) are still not fully understood.Most BHs detected from LIGO/Virgo events have low aligned components of their spins (Abbott et al. 2019a;Zaldarriaga et al. 2018;Miller et al. 2020;Roulet et al. 2021;Zevin et al. 2021), which agrees with predictions of efficient angular momentum (AM) transport within the interiors of massive stars.Such processes remove the majority of AM from the stellar core, predicting slowly rotating remnants after core-collapse (Fuller & Ma 2019;Ma & Fuller 2019).These theories are approximately consistent with core-rotation rate measurements of low-mass red giants from asteroseismology (Beck et al. 2012;Mosser et al. 2012;Deheuvels et al. 2014Deheuvels et al. , 2015;;Triana et al. 2017;Gehan et al. 2018), with a few discrepancies (Eggenberger et al. 2019).Yet, among a small fraction of LIGO/Virgo BHs and a ma-jority of high-mass X-ray binaries (Miller & Miller 2015), high BH spins are measured.Therefore it still remains a theoretical challenge to explain the existence of these rapidly rotating objects (see, e.g.discussions in Qin et al. 2019;Belczynski et al. 2021;Fishbach & Kalogera 2022).
A natural scenario to form high-spin BHs is through binary interactions, as nearly all BHs with spin measurements are found via BH mergers or X-ray binaries.One possible progenitor of BH binaries are Wolf-Rayet-BH binaries.Such a system is formed from an ordinary massive binary system, where the primary collapses to a (likely slowly rotating) BH, and then strips off the envelope of the secondary, making it a Wolf-Rayet (WR) star.Tidal interactions during the WR phase could possibly spin up the latter, forming a rapidly spinning BH.Many studies have investigated this scenario and made predictions for the spins of the second-born BHs (Kushnir et al. 2017;Qin et al. 2018;Bavera et al. 2020;Belczynski et al. 2020;Olejak & Belczynski 2021;Fuller & Lu 2022), finding they can be large for sufficiently close binary systems (P orb ≲ 1 day).

Ma & Fuller
However, in most of these studies, the tidal response of the WR star to the BH companion is not calculated directly.Instead, an effective tidal torque calculated from Zahn's theory of dynamical tides (Zahn 1975(Zahn , 1977, see also Goldreich & Nicholson 1989) is often assumed.The basic picture of Zahn's theory is as follows: gravity waves are tidally excited near the convective coreradiative envelope interface inside a star.The waves propagate outwards and damp due to radiative diffusion near the surface of the star.The damping is often so strong that the waves dissipate before reaching the surface and behave as travelling waves rather than standing waves.The energy and AM deposited by the waves can be calculated and translated to an effective tidal torque.While this picture is often assumed in studies of tidally excited waves, it has not been closely examined in binaries involving a WR star.
In this work, we directly solve for oscillation modes of WR stellar models, quantifying their tidal coupling strengths and dissipation rates.We then compute AM transfer rates and model their spin evolution and resulting BH spins, comparing to those from Zahn's theory.The plan of this paper is as follows: in §2 we review the basic formalism of dynamical tides for calculating tidal torques based on stellar evolution models, and we summarize the setups of our models of the WR stars; in §3 and §4 we present our analysis for the tidally excited modes and the stellar spin evolution.We discuss our results in §5 and conclude in §6.

TIDAL TORQUES BY DYNAMICAL TIDES
In classical tidal theory, tides can be decomposed into two components: equilibrium tides and dynamical tides.The former corresponds to the global distortion of the star, while the latter is composed of internal oscillations, which is believed to be a dominant cause of tidal dissipation.From Ma & Fuller (2021), the energy dissipation rate of a tidally forced oscillation mode α excited by the tidal potential of an aligned and circular orbiting secondary is given by , (1) where ω α and γ α are the mode frequency and damping rate, and ω f = m(Ω orb − Ω spin ) is the tidal forcing frequency (measured in the frame co-rotating with the primary), and Ω spin is the star's angular rotation frequency.M 1 and R 1 are the mass and radius of the primary, q = M 2 /M 1 is the mass ratio of the secondary to the primary, a and Ω orb are the semi-major axis and the angular frequency of the orbit.l and m are the mode's angular and azimuthal wave numbers and W lm is an expansion coefficient of the tidal potential.
α is the dimensionless overlap integral describing the spatial coupling between the mode and the tidal potential, which is calculated by the relation Q α = −(2l + 1)δΦ α /(4πω 2 α ) (Fuller 2017), where δΦ α is the surface gravity potential perturbation.The mode angular momentum dissipation rate is related to the energy dissipation by (Fuller 2017) assuming a circular orbit.Hence, by solving for the internal oscillation modes (with ω α , γ α and Q α ) inside the primary, we are able to calculate the energy dissipation and tidal spin-up rate, given a companion mass and orbit.

Stellar Models
We built our WR star models with the MESA stellar evolution code (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019)).Instead of using the binary options in MESA, we construct our models as follows: we start with a number of zero-age main-sequence (ZAMS) single star models with a variety of masses, summarized in Table 1.The stars evolve to core hydrogen depletion before the stripping-off process occurs.We simulate this process by artificially removing the outer hydrogen envelope immediately after hydrogen depletion (defined by the time when the central hydrogen fraction drops below 10 −5 ), producing a helium star as the initial setup for the WR star.We then restart the evolution until the end of core helium depletion (when the central helium fraction drops below 10 −5 ), and we output the stellar pulsation parameters to be used later for spin-evolution calculations.Example MESA inlists are available on Zenodo under an open-source Creative Commons Attribution license: doi:10.5281/zenodo.7935443,and the model parameters are summarized in Table 1.
During the helium burning phase, we compute the internal oscillations of the models with the GYRE stellar oscillation code (Townsend & Teitler 2013;Townsend et al. 2018;Goldstein & Townsend 2020).We use the second order Magnus differential scheme to calculate non-adiabatic modes, as it proves to be the most reliable when dealing with low-frequency oscillations.We specify our search to l = m = 2 modes since this is the dominant part of the tidal potential in aligned and circular orbits, with the corresponding W 22 = 3π/10.Example GYRE inputs are available on Zenodo under an open-source Creative Commons Attribution license: doi:10.5281/zenodo.7935443.Once we have the mode solutions, we integrate the spin-orbit evolution with Eq. 1 and 2, summing over all modes.We assume the primary remains rigidly rotating during the evolution, due to the strong AM diffusion inside WR stars (Fuller & Lu 2022).We use our non-rotating mode solutions all along in the integration, as we will see that most systems never get to tidal synchronization, such that the rotational effects can be ignored.
An important process related to the spin-orbit evolution is the large wind mass loss experienced by WR stars (e.g., Sander & Vink 2020), which removes AM from both the spin and the orbit.The mass loss rates of high-mass WR stars are somewhat uncertain (Sander & Vink 2020;Vink 2022), especially at low-metallicity, hence there are few reliable observed/modelled values to compare with.We simulate the mass loss with the "Dutch" wind scheme (Nugis & Lamers 2000) in MESA with η = 0.5 and include its effects in our integration.The mass-loss rate has a strong dependence on the metallicity of the star.However, we find that GYRE was unable to solve the oscillations correctly for some of our massive models at solar metallicity due to MESA's artificial treatment of super-Eddington near-surface layers.We hence used a universal metallicity Z = 0.01 Z ⊙ in all our models so that the stellar structure can be more accurately modeled.Oscillations solved from these models are reasonable approximations since the mode properties are mostly determined by the deep internal structure of the stars which are not strongly dependent on metallicity.
To estimate the evolution and mass loss rates of higher metallicity stars, we increase the wind scaling factor to match the mass-loss rate of an alternative model with the desired metallicity.For instance, to simulate a 10 M ⊙ WR star at solar metallicity, we use a wind scaling factor of η = 4 for our Z = 0.01 Z ⊙ model of the same mass, which produces roughly the same mass loss rate as a Z = Z ⊙ model with η = 0.5.In the following, we will simply reference the models with their desired metallicity, yet the readers should keep in mind that the underlying models actually have Z = 0.01 Z ⊙ and adapted Dutch factors, which are summarized in Table 1.
3. MODE MORPHOLOGY Figure 1 shows some example mode eigenfunctions from a 10 M ⊙ WR star model at solar metallicity during the helium burning phase.We see that there is a distinction between the high-frequency (P ≲ 0.8 d, green line) and low-frequency (P ≳ 0.8 d, blue and red lines) modes.For a typical high-frequency mode, the eigenfunction appears to be a standing low-order gravity wave trapped in the radiative envelope of the star.This is in contrast with the Zahn (1975) model for travelling waves that damp near the stellar surface.When the tidal forcing frequency ω f becomes close to the frequency of one of these modes, a resonance occurs and the energy/angular momentum dissipation becomes dominated by it (cf.Equation 1).The tidal torque contributed by this standing mode is different from Zahn's theory, and we will see in §4 that Zahn's results on spinorbit evolution are significantly altered.
At lower frequency, the modes turn to travelling waves (blue line) as Zahn assumed, since the g-mode dispersion relation indicates an imaginary wave number Im(k r ) ∝ 1/ω 2 α (see Appendix B), i.e., the spatial evanescence becomes larger at lower frequency.For even lower frequency, the eigenfunction becomes a mixed mode which can be separated into two components: a gravity wave inner region of the radiative envelope, and a thermal wave region in the outer envelope.We show in Appendix A and B that the transition occurs around a critical frequency ω crit ≡ (4λN 2 ω T ) 1/3 , where ω T = κ/r 2 is the thermal frequency, and κ = 16σ B T 3 /(3ρ 2 c P κ R ) is the thermal diffusivity.The thermal mode exists where |ω crit | > ω α , while the g-mode exists where |ω crit | < ω α , as seen in the right panel of Figure 1.For higher frequency waves the thermal mode region becomes narrower and disappears.Since the thermal mode components only exist near the very surface of the star, where  the density is very low, we would expect that the mechanical torques are mostly contributed by the travelling g-mode component excited in the deep interior as Zahn's theory assumed.Hence at low tidal forcing frequency the tidal torques should be similar to Zahn's model, as we will see in §4.
Figure 2 shows the periods, damping rates and overlap integrals Q α of all modes we solved for the same 10 M ⊙ model.We see that most modes have a nearly constant period spacing, matching the expectations for g-modes.The damping rates for most modes are at the same order of magnitude, except at high-frequency where the damping is significantly lower.This is due to their low radial wave numbers k r and the damping rate γ α ∝ star k 2 r κ|ξ α | 2 dm.Low-frequency modes become traveling waves whose damping rate is roughly the wave crossing time.
The overlap integral |Q α | typically decreases as the mode frequency decreases, but with significant scatter and with "hills and valleys" as the frequency decreases.Since the on-resonance AM dissipation Jα ∝ γ α |Q α | 2 , we expect to see the same "hills and valleys" features in the tidal synchronization rate, as the tidal forcing frequency is moving across different modes with varying Q α .This is also different from Zahn's theory, which predicts a "smooth" power-law relation for the AM deposition rate as a function of orbital period (Kushnir et al. 2017, or Equation 6 in this paper).
In the frequency ranges where mixed modes appear, we notice that GYRE suffers from numerical convergence problems as it starts looking for extremely highorder modes.We identified some of the mode solutions in that regime as "strange modes", and an example is labelled in Figure 2.These modes often have excess damping rates and unusual winding numbers (mode radial order), and do not obey the usual frequency spacing of g-modes.In addition, their eigenfunctions appear to be artificially truncated as they reach deeper inside the star, unlike other modes with an inner g-mode compo-nent at similar frequencies, which are truncated near the convective core boundary.We are not sure if these modes are physical or caused by numerical artifacts from GYRE, hence we do not include them in the spin-orbit integration.A detailed discussion of these modes is presented in Appendix C.

EVOLUTION OF WR SPINS
We integrate the spin-orbit evolution of WR-BH binaries from the WR star models and oscillations modes we have computed.Throughout the evolution, the orbital AM of the system is lost due to winds from the primary, gravitational radiation and tidal AM transfer: where Jwind,orb = Ṁ1 Ω orb (M 2 a/(M 1 + M 2 )) 2 and Jtide is given by summing over all modes from Equation 2. At short orbits, the orbital decay timescale by gravitational wave radiation is given by Peters (1964) (assuming circular orbits) where q = M 2 /M 1 is the mass ratio.This gives t GW ≈ 206 Myr × (M 1 /10 M ⊙ ) −5/3 (P orb /0.3 d) 8/3 for equal mass binaries (q = 1), much greater than the typical WR lifetime (≲ 1 Myr).Hence, gravitational radiation is not important in our case, but we still include the term JGW = (32/5)(G/a) The primary receives spin AM from the orbit at Jtide , and it loses AM through winds: where Jwind,spin = Ṁ1 Ω spin R 2 1 .The spin of the primary may also change due to the changes of its internal structure and hence moment of inertia.Since the secondary is a BH in our case, its spin is not coupled.
For comparison, we also integrate each evolution based on Zahn's formalism, with an adapted AM transfer rate from Kushnir et al. (2017): where s c = 3/πGρ c |Ω orb − Ω spin |, while r c , ρ c and ρc are the convective core radius, the density at the core boundary, and the average density of the core, respectively.
We construct the integration machinery of the spinorbit evolution as follows: after generating a grid of stellar model snapshots throughout the star's evolution, we solve for oscillation modes for each snapshot with GYRE.We begin our integration at the start of the helium burning phase (defined by the instant when 2% of the core helium burning lifetime has passed).We carefully apply an adaptive time step control to avoid i) sudden crossing of resonance locations; ii) sudden changes of mode frequencies; iii) changes of more than 2% of the total evolution phase lifetime; iv) sudden change of stellar spin by 2%, in one time step.To evaluate physical quantities (e.g.mode frequencies, stellar masses) between two model snapshots, we estimate them by interpolating these snapshots and their corresponding GYRE solutions.In doing so, we track the modes by their radial orders n pg , and only include the mode eigenfunctions existing in both snapshots.We carried out resolution tests with half our selected timesteps and we confirm that the results are nearly identical.
Theories have suggested that the strong magnetic coupling between the stellar core and envelope (e.g.Taylor-Spruit dynamo, Spruit 2002;Fuller et al. 2019) removes the majority of core AM immediately after the main sequence (Ma & Fuller 2019), before the envelope can be stripped off.Hence, we assume initially non-rotating WR stars.We run models with initial orbital periods of 0.3, 0.5, or 0.8 days.We find that longer period orbits exhibit very little tidal spin-up.
In this work we specify our calculations to equal mass binaries (q = 1), since they are the most relevant for binary black holes.For cases with different mass ratios, one would expect from Equation 1 that the tidal dissipation rate (hence the tidal spin-up rate) naïvely scales as q 2 , as we verified in some additional test runs.However, extreme mass ratios could allow for orbital decay and resulting processes such as resonance locking or binary mergers.We hope to generalize our calculations to such systems in future works.
Figure 3 shows the spin and orbital evolution for a few of the systems we studied.Since the WR primary burns helium and loses mass throughout the evolution, the mass and central helium fraction can be seen as time coordinates, as shown.For systems with short initial orbits and high metallicities (P orb,i = 0.3 d, Z = Z ⊙ , left panel), we see that the primaries get significantly spun up, yet they are not tidally synchronized even at the end of evolution.In addition, the final spins of massive models are never higher than what one would expect from Zahn's theory.
When we increase the initial periods (P orb,i = 0.8 d, middle panel of Figure 3), the tidal torques decrease as expected.At these long periods, the tidal torque is dominated by traveling waves and the results agree well with  6) is assumed.Left: Systems with initial orbital periods of 0.3 days and a mass-loss rate equivalent to solar metallicity.Mass loss is very significant for high-mass models and the final spins depart from Zahn's formalism for them.Middle: Systems with initial orbital periods of 0.8 days and a mass-loss rate equivalent to solar metallicity.Mass loss overpowers tidal spin-up and the primaries are not spun up much, consistent with Zahn's results.Right: Systems with initial orbital periods of 0.3 days and a mass-loss rate equivalent to 0.01 solar metallicity.Mass loss is negligible for most systems and the primaries are partially spun up, but not as much as Zahn's formalism predicts.None of these models reach tidal synchronization.Zahn's formalism.However, the tidal torque is unable to compete with mass loss, which almost completely removes the spin AM the star accumulated during the first half of the evolution, leaving a slowly spinning primary.Tidal spin-up is followed by mass-loss induced spindown in the middle of the evolution for solar-metallicity models (Figure 3, left and middle panel) due to an increase in the wind mass loss rate.This occurs when the helium envelope is lost completely, exposing the CO-rich The mode frequencies increase as the star evolves, while the forcing frequency decreases, preventing resonance locking.Ωspin increases rapidly while ω f decreases rapidly at resonance crossings, creating the "step-like" features we see in the spin frequency evolution (Figure 3).
core, and greatly increasing the mass loss rate according to the "Dutch" mass loss prescription (Figure 4).In several thousand years the winds remove the star's spin AM until the mass loss rate decreases somewhat, allowing tidal spin-up to proceed.However, ongoing mass loss and a widened orbit prevent tides from spinning up the star to synchronization.
When we consider short initial periods but move to low-metallicity models (P orb,i = 0.3 d, Z = 0.01 Z ⊙ , right panel of Figure 3), the mass loss becomes negligible and the spin evolution is dominated by tidal effects.The orbits do not change significantly.We see that the primaries get significantly spun up, yet still not reaching tidal synchronization, and the resulting spin is much slower than Zahn's prediction, except for the lowest mass models.This is because the transition period from standing modes to travelling waves increases as the stellar mass increases.Hence the evolution is more likely to depart from Zahn's formalism for more massive primaries (see §5.1).
The evolution of spin frequencies show "step-like" features (most easily seen in Figure 3 middle panel) characterized by sudden increases in spin frequency.This is caused by the resonance-crossing of standing modes with the tidal forcing, as illustrated in Figure 5.When the tidal forcing frequency gets close to one of the mode frequencies, a near-resonance occurs (cf.Equation 1) and the tidal torque drastically increases, leading to high spin-up rate.The occasional crossings of these resonances create the "step-like" features.

Comparison to Zahn's Formalism
To understand why and how the tidal evolution differs from Zahn's formalism, in the upper panel of Figure 6 we show the tidal torques calculated for a 10 M ⊙ WR star model for different tidal forcing periods P tide = 2π/Ω tide with these two approaches.For short tidal periods, the tidal torque has sharp peaks at standing g-mode frequencies, in contrast to the power-law dependence of Zahn's prediction.This is caused by low damping rates of standing modes at short mode periods (cf. Figure 2 and Figure 6 lower right panel).
As the mode and orbital frequencies evolve, resonance crossings occur, hence the accumulated tidal spin-up must be evaluated by integrating Jtide over time.Zahn's theory for travelling waves generally overestimates the tidal spin-up in this case, as the resonance peaks are narrow, and during the majority of evolution, the tidal torque is much less than what Zahn's formalism estimates.This helps to explain the departure of final spins shown in Figure 3 for short-period systems (left and right panel) from Zahn's predictions.
When the systems are in long-period orbits, or already at a stage where Ω spin becomes comparable to Ω orb (close to tidal synchronization), the tidal forcing periods become long and the tidal torques are mostly contributed by modes at low frequencies.These modes, in contrast to their high-frequency partners, have large damping and are essentially travelling waves (cf. Figure 2 and Figure 6 lower right panel).Therefore, the tidal torque is no longer dominated by resonance with an individual mode, but instead has contributions from many strongly damped modes, effectively forming a "continuum" (Figure 6, red dot).This continuum formed by traveling waves is exactly what Zahn's formalism assumes, so the torques should be similar to Zahn's formalism, which is confirmed in Figure 6.This explains the consistency between our results and Zahn's for longperiod initial orbits (Figure 3 middle panel).
We note, however, that the transition period between standing waves and travelling waves depends on the stellar mass: for higher mass models, the frequency range for standing waves extends to longer periods, as shown for the 38 M ⊙ model in the lower left panel of Figure 6.Hence, the tidal evolution of massive WR stars departs more strongly from Zahn's formalism, as we see in the left and right panels of Figure 3.
This distinction is caused by the different structures of low and high-mass WR stars.For higher mass stars, a larger fraction of the total internal pressure is contributed by radiation pressure since they are hotter and more luminous.Radiation pressure, however, contributes smaller buoyancy forces because the Brunt-Väisälä frequency N is zero for a star supported purely by radiation pressure.Indeed, the Brunt-Väisälä frequencies within our high-mass models are smaller than those within our low-mass models (Figure 7).This increases the g-mode period spacing (proportional to N −1 ) of high-mass stars and decreases the radial wave number at the same frequency (as k r ∝ N ).Hence, the mode damping rate γ α ∝ star k 2 r κ|ξ α | 2 dm is also decreased.A secondary effect is that the convective cores are larger in more massive stars, making the radiative envelopes and g-mode cavities narrower.These combined effects make the resonance peaks narrower and more widely spaced for higher mass models, and further from the travelling wave limit of Zahn's formalism.

Resonance Locking
When the tidal dissipation is dominated by resonant modes, an important process called resonance locking may occur (Witte & Savonije 1999).However, we argue that this is unlikely to occur for WR-BH binaries.Resonance locking can happen when the mode's frequency evolves at the same rate as the tidal forcing frequency: We see from our example evolution tracks (Figure 3) that in most cases Ωorb < 0 and Ωspin > 0, which means ωf < 0, in contrast to the fact that ωα > 0 due to stellar At long forcing periods (e.g., the red dot), the tidal torque is no longer dominated by individual modes, but arises from a multitude of highly-damped travelling modes.This is exactly Zahn's assumption, and hence the torques are similar at long periods.Lower Left: The same as the upper panel but for a more massive 38 M⊙ Wolf-Rayet model with 90% central helium abundance.We see that the mode-period spacing becomes larger, and the standing wave region extends to longer tidal forcing period.Lower Right: Eigenfunctions of the most resonant mode at the green and red dots in the upper panel.We see clearly that one is a standing mode while the other is a travelling mode.
evolution (Figure 5).Hence the above relation never holds and resonance locking can never happen.Instead, the system rapidly passes through resonances, creating the step-like features in Figure 3.We note that when mass loss dominates the spin evolution, we could occasionally have Ωspin < 0 (Figure 3, left and middle panels) and resonance locking may happen during this phase.This may prevent a star from rapidly spinning down, but it cannot cause tidal spin-up.However we find that during these phases the Brunt-Väsäilä frequency of the star increases rapidly, such that ωα exceeds ωf even at resonance, in contrast to the resonance locking criterion.Hence, resonance locking does not appear to occur in any of our modeled systems.

Implications for BH Spins
A rapidly rotating WR star can probably collapse to a fast-spinning BH, forming a high-spin binary BH system.If angular momentum is conserved during the corecollapse process, the dimensionless spin parameter of the resulting BH is In Figure 8, we show the resulting black-hole spins of our WR star models, assuming that they preserve their masses and angular momenta after helium burning and during the core-collapse process.We also show the predicted BH spins with Zahn's formalism.We see that lower-mass systems can form faster-spinning BHs, as their tidal spin-up is more efficient.It is only in ultra-short orbits that these systems form fast-spinning BHs.
For solar metallicity systems, tidal spin-up cannot overcome AM loss from winds, resulting in low spins for systems starting at long (P orb,i ≳ 0.5 d) orbital periods.Low-metallicity (1% Z ⊙ ) systems with 0.5 d ≲ P orb,i ≲ 1 d produce larger BH spins with 0.1 ≲ a ≲ 0.8, compared to a ∼ 0.01 predicted by single star evolution models (Fuller & Ma 2019).For high-mass systems, the spins are much smaller than Zahn's predictions.
If the companion black hole is assumed non-spinning, our predicted BH spins will be roughly compatible with some LIGO measurements with moderate spins (0.1 ≲ χ eff ≲ 0.5) (Abbott et al. 2019a(Abbott et al. ,b, 2021)), but would have a tough time matching any events with large χ eff .The relationship between orbital period, mass, and spin is different than what Zahn's theory predicts.Whereas we typically find higher spins for M BH ≲ 10 M ⊙ , Zahn's theory predicts smaller spins for lower mass BHs.There may be an anti-correlation between mass and spin (Safarzadeh et al. 2020) which would support our new models.A mass-spin correlation from future LIGO-VIRGO data will help distinguish between these models.None of our high-mass models predict spins comparable to some high-spin measurements (a ≳ 0.9) from X-ray binaries (Miller & Miller 2015), and the uncertainty of such measurements is still under debate (Belczynski et al. 2021;Fishbach & Kalogera 2022).However, those measurements are for the first-born BH, while our models only apply to the second-born BH.
Previous works on tidal interactions between WR-BH binaries have predicted black hole spins similar to our "Zahn's results" in Figure 8, in which Zahn's formalism is assumed (Qin et al. 2018;Bavera et al. 2020;Belczynski et al. 2020;Olejak & Belczynski 2021;Fuller & Lu 2022).These results likely overestimate the black hole spins when standing waves are present, which applies primarily to massive BHs (M BH ≳ 10 M ⊙ ).Detmers et al. 2008 also investigated tidal spin-up of WR stars, but used different prescriptions for tidal dissipation, winds, and orbital AM losses.Unlike our results and those listed above, they found that tidal spin-up coupled with mass loss frequently caused the orbits to decay and instigate mass transfer.This outcome is more likely with small companion masses (e.g., neutron stars) whose orbits must decay more in order to tidally spin-up the WR star.

Nonlinear dissipation
Throughout the paper, we have assumed that the tidally forced modes are linear.However, this is not true when a mode is close to resonance, especially for massive models with larger on-resonance mode ampli- .The Brunt-Väisälä frequency (solid lines) and the ratio of radiation pressure to total pressure (dashed lines) of two WR models of 10 (blue) and 38 M⊙ (red), the same models in Figure 6.The radiation pressure fraction is higher for the more massive model, making its Brunt-Väisälä frequency lower.This causes its resonance peaks to be narrower and more separated, as seen in Figure 6.
tudes.To examine how nonlinearity may affect our conclusions, we estimate the nonlinear damping rate for a mode α in Appendix D. For weakly nonlinear modes, an approximate nonlinear damping rate may be where the numerator is the mode nonlinearity (i.e., the peak value of dξ r /dr within the star, which is much less than unity for a weakly nonlinear mode) and the denominator is the wave crossing time of the envelope.We rerun our evolution models with this nonlinear damping rate, and we find that nearly all models achieve more tidal spin-up compared to linear damping only.
In Figure 9, we show the predicted BH spins for our systems with nonlinear damping, compared to the predictions from Zahn's theory.We see that for very shortperiod orbits (P ∼ 0.3 d), the strong tidal forcing triggers substantial nonlinear damping, spinning up the BHs to nearly the same rotation rates as predicted by Zahn's theory, where maximum damping occurs.Nonlinear effects are the most significant for low-metallicity and high-mass (M ≳ 30 M ⊙ ) models, which have lower order g-modes dominating their tidal processes and less linear damping (see discussion in §5.1).
However, our nonlinear damping model is crude, so these predictions are not very reliable.A more detailed study of the nonlinear interactions has to be carried out to establish firm conclusions for the final BH spins.

Caveats
Throughout this paper, we have assumed very efficient angular momentum transport within the WR star, such that it remains rigidly rotating.This is justified by the asteroseismically callibrated models of magnetic angular momentum transport (Fuller & Ma 2019) that predict nearly rigid rotation during the helium-burning phase (Fuller & Lu 2022).However, if angular momentum transport is inefficient, gravity waves damping near the stellar surface will preferentially spin up those layers until they are synchronized.This will create a critical layer at which subsequent waves are absorbed (Goldreich & Nicholson 1989), synchronizing the star from the outside inwards.Recent works have investigated the formation of critical layers and subsequent absoroption of incoming waves (Su et al. 2020;Guo et al. 2023), though they do not include magnetic torques that may allow angular momentum transport to prevent critical layer formation.If a critical layer can form, it will absorb outgoing waves, such that Zahn's model applies once again.
We have ignored the influence of the Coriolis force in our calculations.This will become significant once the star's have been partially spun up and the tidal forcing frequency becomes smaller than the rotation frequency.However, the prograde ℓ = m = 2 modes that dominate the tidal interaction have eigenfunctions that are only slightly changed by Coriolis forces (see, e.g., Fuller 2014), so we don't expect any of our conclusions to be greatly affected.
We have adopted the "Dutch" wind models with an artificial scaling factor to simulate the mass-loss rates for Wolf-Rayet stars at different metallicities.However, the wind physics for these stripped stars are highly uncertain (Vink 2022), and different wind models could result different rates for the removal of spin angular momentum from the primary, introducing uncertainties in the final black hole spins.Nevertheless, we don't expect these uncertainties to exceed the differences between our solar-metallicity models and the 0.01 Z ⊙ models, as they represent extreme cases of large and negligible massloss, respectively.Hence, the final black hole spins with the "correct" wind physics should lie between the data points representing models with the same initial mass and periods but different metallicities in our Figure 8.The conclusion that these black holes are not spun up to maximal rotation appears robust.
Our stellar models were run at low metallicity in order to reliably calculate the near-surface structure and mode eigenfunctions.Higher metallicity stars will have somewhat different structure and mode eigenfunctions near the surface, particularly around the iron group opacity peak.If this significantly affects mode damping rates, then the tidal synchronization efficiency would be similarly altered.We set up additional test models with 0.02, 0.03 and 0.2 Z ⊙ and find that the oscillation Porb,i = , , days mass-loss equivalent to: : Z , : 0.01Z filled: ours, empty: Zahn's formalism 0.3 0.5 0.8 Figure 8.The dimensionless spin of resulting black holes for our Wolf-Rayet star models if their masses and angular moat the end of helium burning is preserved.Zahn's predictions for the same system are also shown, and are assumed 1 (maximum rotating) if the progenitors have J > GM 2 /c.For short initial orbits, the models typically predict higher spins than individual stellar evolution models, where the spins could be as low as 10 −2 .However, the spins are usually lower than Zahn's predictions, especially for high-mass systems.
mode parameters (frequencies, damping rates and overlap integrals) show no significant differences, nor specific trends towards higher metallicities, hence we expect our treatment to be appropriate.However, these models all have weak winds, while the strong winds in solar-metallicity models may also alter the eigenmode properties.Ro & Matzner (2016); Ro (2019) present detailed models of the transition from the hydrostatic star the hydrodynamic wind in the near-surface layers.Future work should investigate how those types of stellar models affect mode eigenfunctions and damping rates.
Finally, our calculations are performed by summing up the contribution of individual tidally excited oscillation modes.If non-resonant modes outside our computed frequency range contribute to the tidal dissipation, or if our eigenmode calculations miss highly non-adiabatic thermal modes that contribute to the dissipation, then the tidal dissipation rate could possibly increase.It would be interesting to compare to calculations performed by directly computing the forced tidal response, as outlined in (Sun et al. 2023).

CONCLUSION
In this work, we investigate the dynamical tidal spinup of Wolf-Rayet stars from black hole companions.We build Wolf-Rayet star models with different metallicities, and then calculate their oscillation mode frequen- Figure 9.The dimensionless spin of resulting black holes for our Wolf-Rayet star models with the nonlinear damping rates we estimate in Appendix D. Zahn's predictions are also shown.The black hole spins are significantly increased when including nonlinear damping, especially for high-mass black holes and short-period systems.Hence, nonlinear effects could be an important factor in these systems.
cies, damping rates, and eigenfunctions.We use these to integrate the coupled spin-orbit evolution of the binary based on the tidal excitation of these oscillation modes.We also make predictions for the resulting BH spins upon core-collapse of the Wolf-Rayet star.
We study the properties of the oscillation modes and find that at shorter orbital period, the tidal forcing is mostly contributed by standing g-modes, in contrast to the usual assumption of travelling waves proposed by Zahn (1975).The standing g-mode spectra contributes a resonance structure, and during most of the spinorbit evolution phase, the tidal response lies between resonances and the interaction strength is weaker than Zahn's prediction.The tidal forcing transits to Zahn's travelling wave limit at longer periods, in which Zahn's estimate is more accurate.However, the specific transition frequency depends on the stellar masses, and the structure for more massive stars (supported significantly by radiation pressure) tend to have lower transition frequencies, allowing systems in longer-period orbits to evolve differently compared to Zahn's prediction.
We find that it is difficult to tidally synchronize Wolf-Rayet stars during helium-burning.For solar-metallicity Wolf-Rayet stars, strong winds tend to remove the majority of angular momentum deposited by tides, leaving slowly spinning stars and black holes.At low metallicity, the stellar wind is weaker and the stars are significantly tidally spun up, yet still less than Zahn's prediction, especially for massive stars and short-period orbits.
Tidal interactions can significantly spin up the resulting BHs compared to single-star models.Yet the predicted black hole spins a are still ≲ 0.4 for all but our shortest period (P orb ≲ 0.5 d) models.These predictions are consistent with some low/moderate-spin measurements from LIGO/Virgo black hole merger events, but cannot explain high-spin X-ray binaries events since only the second-born black hole has large spin in these models.
We have discussed a new class of gravito-thermal modes that appear in our calculations, yet we do not reach a firm conclusion whether these modes are physical or caused by numerical artifacts.Future work should investigate the origins of these modes, and any effect they could have on tidal spin-up.
There are also caveats to our work.We have assumed rigid rotation of the Wolf-Rayet star during our spin-orbit evolution calculations, as expected if there are strong internal AM transport processes caused by magnetic torques.However, weak AM transport could enable surface critical layer formation, allowing Zahn's model to apply.We did not realistically calculate the near-surface structure of our solar-metallicity models, which could alter our estimate of the mode damping rates.We also point out that nonlinear damping effects could be significant for our most massive models, which can produce more tidal dissipation than our predictions from linear theory.This should be studied and improved in future work.
It is straightforward to verify that Equation A12 reduces to the energy equation in Ma & Fuller (2019) when η = 0, i.e. radiation is neglected.
We now consider the dynamics of the fluid.The perturbed momentum equation reads: where we again used the WKB approximation ∇ r → ik r .The continuity equation with the incompressible approximation1 gives ∇ When the angular dependence of perturbation variables are expanded in spherical harmonics, we have ∇ 2 ⊥ → −λ/r 2 where λ = l(l + 1).Combining the above equations with the energy equation, we arrive at the dispersion relation where H ≡ P/(ρg) is the pressure scale height.With the WKB approximation, k ≈ k r and k r H ≫ 1, the first term in the bracket of the right hand side can be neglected, and the dispersion relation becomes B. GRAVITY AND THERMAL MIXED MODES When gas pressure is non-negligible, we always have 1 − η ∼ 1 and the second term in the bracket of the left hand side of Equation A16 can usually be neglected under WKB approximation k r H ≫ 1.The dispersion relation hence reduces to the quadratic equation where ω T = κ/r 2 is the thermal frequency, and is the critical frequency between different types of modes.The solution to Equation B17 is ) which has two important limits.
1. High-frequency region (|ω crit | ≪ |σ|): the solution reduces to The "+" sign solution further reduces to the (radial) thermal diffusion solution k 2 r κ ≃ σ.The "−" sign solution reduces to k 2 r = −λN 2 /r 2 σ 2 , which is the g-mode dispersion relation.With σ = i(ω + iγ), under the weakly damped limit γ ≪ ω, we have  10.kr for the thermal mode on the complex plane.The four solutions correspond to rapidly increasing, slowly increasing, rapidly decreasing and slowly decreasing thermal modes, respectively.
Since ω crit explicitly depends on the local stellar properties, modes of a given frequency can behave as either gravity or thermal waves in different parts of the star, as we see in Figure 1.Modes can therefore behave as "mixed modes", with gravity mode character in the core of the star where thermal diffusion is unimportant, and thermal mode character near the surface of the star where thermal diffusion is very important.Such modes have rarely been examined in asteroseismology because their high damping rates mean that they will not be visible as stellar pulsation modes.However, these damping rates also mean they could be very important for energy dissipation via tidal excitation.

C. STRANGE MODES
When solving for high-order, low-frequency mixed modes, GYRE occasionally returns solutions which we identified as "strange modes".An example is given in Figure 11.The strange modes are usually distinct in the following aspects: 1) the modes have higher damping rates, often one order of magnitude larger than the normal modes.This causes stronger spatial evanescence as the waves propagate inwards, as seen from Equation B22. 2) The modes have unusual winding numbers n pg (defined in Takata (2006), and treated as mode radial orders in GYRE), departing from the normal n pg -period relation of normal modes (Figure 11,right panel).3) The modes do not obey the uniform period spacing shared by normal g-modes.4) The strange mode eigenfunctions seem to be artificially truncated in the radiative envelope, once they reach a minimum amplitude.We confirm that there are no special physical conditions inside the star where they are truncated.Resolution tests also show that the strange mode solutions do not converge even at very high spatial/frequency resolution.In addition to the normal modes (black lines, including standing g-modes, travelling g-modes and mixed modes, discussed in §3), a strange mode solution (red line) exists.This mode has much higher damping rate because it is localized near the stellar surface, as indicated by its rapidly decreasing amplitude towards the core.Right: The periods and winding numbers (npg) for mode solutions, showing an outlying strange mode solution.
We guess that the strange modes are effectively gravito-thermal mixed modes that are trapped in the near surface region where the waves behave as thermal waves (|ω crit | ≫ |σ|).Because they are trapped in the surface layers, their damping rates are much larger than normal modes, similar to the acoustic strange modes found at high frequencies (Glatzel & Kaltschmidt 2002) .Because their eigenfunctions evanesce so rapidly towards the core, their amplitudes apparently drop below the numerical precision of GYRE near the core, causing the artificial behavior of the eigenfunction at small radii seen in Figure 11.This also causes the value of n pg computed by GYRE to be incorrect, and explains why their exact frequencies/eigenfunctions do not converge at high spatial resolution.
In our calculations, the strange modes only exist in the low-frequency range of mode spectra.Hence, if actually physically present, these modes will only be relevant at the late stage of spin-orbit evolution, when the star has already been significantly spun up.Hence, we believe our main results to be robust against the uncertainties surrounding strange modes, but these modes should be studied in more detail in future work.

D. NONLINEAR DAMPING OF MODES
Our tidal calculations are based entirely on linear theory, yet under certain circumstances nonlinear effects could be important.The dominant nonlinear term in the fluid momentum equation is ξ • ∇ξ ∼ ξ(dξ r /dr), hence dξ r /dr serves as an approximate measure of linearity, which only holds when dξ r /dr ≪ 1.For our most massive models, the modes become nonlinear very close to resonance, so nonlinear effects will be important during resonance crossings.While developing a complete nonlinear theory is beyond the scope of this work, here we propose an ad-hoc estimate of the nonlinear damping rates of modes.

Figure 1 .Figure 2 .
Figure 1.Left: Example mode eigenfunctions for a 10 M⊙ Wolf-Rayet star model at solar metallicity during helium burning.For high-frequency modes, we see a standing g-mode (green line) excited near the convective core (red region) boundary, in contrast with Zahn's assumption of travelling waves.As the frequency decreases, the modes become travelling gravity waves (blue line), damping near the surface (Zahn's formalism).When the frequency continues to decrease, the modes become mixed modes with a travelling g-mode component and a thermal mode component.The red star marks the transition point, calculated by the local maximum of the eigenfunction.Right: A detailed look at the transition points between g-modes and thermal modes (stars).The lines show the frequencies of all modes and the colored lines correspond to the example modes in the left panel.The transition points agree well with the theoretically derived ones where ωα = ωcrit (the boundary between two propagation regions, cf.Appendix B).At higher frequency the thermal mode region (shaded) gets narrower and disappears.

Figure 3 .
Figure 3.The spin and orbital evolution for our Wolf-Rayet-BH binaries.All systems have equal mass companions initially.The solid and dashed lines show the spin and orbital frequencies, respectively.Line colors indicate the evolving central helium mass fraction.The stars mark the end of evolution (core helium depletion), and the red "ZF" symbols show the spins if Zahn's formalism (Equation6) is assumed.Left: Systems with initial orbital periods of 0.3 days and a mass-loss rate equivalent to solar metallicity.Mass loss is very significant for high-mass models and the final spins depart from Zahn's formalism for them.Middle: Systems with initial orbital periods of 0.8 days and a mass-loss rate equivalent to solar metallicity.Mass loss overpowers tidal spin-up and the primaries are not spun up much, consistent with Zahn's results.Right: Systems with initial orbital periods of 0.3 days and a mass-loss rate equivalent to 0.01 solar metallicity.Mass loss is negligible for most systems and the primaries are partially spun up, but not as much as Zahn's formalism predicts.None of these models reach tidal synchronization.

Figure 4 .
Figure 4.The mass and spin evolution of a 38 M⊙ Wolf-Rayet star at solar metallicity, with an initial orbit of 0.3 days.The shaded regions show the dominant composition as a function of mass coordinate (right axis).At an age of ∼ 3.3 Myr, mass loss exposes the carbon-rich core, greatly enhancing the mass loss and spin down rates.

Figure 5 .
Figure5.The evolution of the mode frequencies (black lines) and tidal forcing frequency (red line) of a WR-BH binary with a 26 M⊙, 1.0 Z⊙ WR star, an equal mass companion and an initial orbit of 0.8 days.The mode frequencies increase as the star evolves, while the forcing frequency decreases, preventing resonance locking.Ωspin increases rapidly while ω f decreases rapidly at resonance crossings, creating the "step-like" features we see in the spin frequency evolution (Figure3).

Figure 11 .
Figure11.Left: All mode eigenfunctions for mode periods between 0.1 to 2 days solved by GYRE for a 10 M⊙ Wolf-Rayet star model at solar metallicity during helium burning.In addition to the normal modes (black lines, including standing g-modes, travelling g-modes and mixed modes, discussed in §3), a strange mode solution (red line) exists.This mode has much higher damping rate because it is localized near the stellar surface, as indicated by its rapidly decreasing amplitude towards the core.Right: The periods and winding numbers (npg) for mode solutions, showing an outlying strange mode solution.

Table 1 .
Parameters of our Wolf-Rayet star models.We fixed the metallicities of all models to Z = 0.01 Z⊙ and adapted their Dutch wind scaling factors to match the massloss rates for the desired metallicities.See discussions in the main text.