Hot subdwarf B sdB stars are stripped heliumburning stars that are often found in close binaries where they experience strong tidal interactions. The dissipation of tidally excited gravity waves alters their rotational evolution throughout the sdB lifetime. While many sdB binaries have wellmeasured rotational and orbital frequencies there have been few theoretical efforts to accurately calculate the tidal torque produced by gravity waves. In this work we directly calculate the tidal excitation of internal gravity waves in realistic sdB stellar models and integrate the coupled spinorbit evolution of sdB binaries. We find that for canonical sdB M SUBsdBSUB  0.47 M SUBSUB binaries the transitional orbital period below which they could reach tidal synchronization in the sdB lifetime is 0.2 day with weak dependence on the companion masses. For lowmass sdBs M SUBsdBSUB  0.37 M SUBSUB formed from more massive progenitor stars the transitional orbital period becomes 0.15 day. These values are very similar to the tidal synchronization boundary 0.2 day evident from observations. We discuss the dependence of tidal torques on stellar radii and we make predictions for the rapidly rotating white dwarfs formed from synchronized sdB binaries.
['B subdwarf stars', 'Stellar oscillations', 'Stellar rotation', 'Tidal interaction', '129', '1617', '1629', '1699', 'Astrophysics - Solar and Stellar Astrophysics', 'Astrophysics - High Energy Astrophysical Phenomena']
Tidal Spinup of Subdwarf B Stars
{'Tidal Spin-up of Subdwarf B Stars': '<!-- image -->  \n1 TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA  \n2  \nDepartment of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA', 'ABSTRACT': 'Hot subdwarf B (sdB) stars are stripped helium-burning stars that are often found in close binaries, where they experience strong tidal interactions. The dissipation of tidally excited gravity waves alters their rotational evolution throughout the sdB lifetime. While many sdB binaries have well-measured rotational and orbital frequencies, there have been few theoretical efforts to accurately calculate the tidal torque produced by gravity waves. In this work, we directly calculate the tidal excitation of internal gravity waves in realistic sdB stellar models and integrate the coupled spin-orbit evolution of sdB binaries. We find that for canonical sdB ( M sdB = 0 . 47 M ⊙ ) binaries, the transitional orbital period below which they could reach tidal synchronization in the sdB lifetime is ∼ 0 . 2 days, with weak dependence on the companion masses. For low-mass sdBs ( M sdB = 0 . 37 M ⊙ ) formed from more massive progenitor stars, the transitional orbital period becomes ∼ 0 . 15 days. These values are very similar to the tidal synchronization boundary ( ∼ 0 . 2 days) evident from observations. We discuss the dependence of tidal torques on stellar radii, and we make predictions for the rapidly rotating white dwarfs formed from synchronized sdB binaries.  \nKeywords: B subdwarf stars(129), Stellar oscillations (1617), Stellar rotation(1629), Tidal interaction (1699)', '1. INTRODUCTION': "Hot subdwarf B (sdB) stars, first observed by Humason & Zwicky (1947), are compact and faint stars with surface temperatures between 20,000 and 40,000 K and masses below 0 . 5 M ⊙ (Heber 2009, 2016; Zhang et al. 2009). These stars have helium-burning cores and thin envelopes (Gotberg et al. 2018), and they are thought to be stripped cores of helium-burning red giants, whose envelopes were previously lost due to some binary interactions (Han et al. 2002, 2003; Pelisoli et al. 2020). About half of the observed sdB systems are found in binaries (Maxted et al. 2001; Napiwotzki et al. 2004; Geier et al. 2011; Copperwheat et al. 2011), with many of them in close ( P orb ≲ 10 d) orbits. This suggests that a prior common envelope phase might be responsible for the ejection of their envelopes, as well as the inspirals of their companions to their current orbital configuration (Kruckow et al. 2021).  \nCorresponding author: Linhao Ma  \nlinhaoma@princeton.edu  \nFor binaries with short orbital periods, tidal interactions can shape both the migration of their orbits and the spin evolution of individual stars. Historically, sdB binaries are sometimes assumed to have reached tidal synchronization, such that their orbital parameters can be derived from measurements of sdB rotation rates (see, e.g., Kudritzki & Simon 1978; Geier et al. 2010), even if the companions (typically white dwarfs or M dwarfs) are too faint to be seen. However, this assumption has been seriously challenged by observations from the past decade, especially those with high-precision measurements with TESS and Kepler /K2, where both spin and orbital frequencies are available (see the summary of observation results in Figure 3). These studies have found that sdB binaries with orbital periods as short as ∼ 7 hours are not always synchronized (Silvotti et al. 2022). Nevertheless, these emerging new data provide an excellent opportunity to test the theoretical modeling of tidal interactions in sdB binaries.  \nFor stars with convective cores and radiative envelopes like sdBs, tidally excited gravity waves in their envelopes are thought to be the most efficient form of tidal interaction (Zahn 1977). These waves are excited by the tidal  \npotential from the companion, and when they propagate through the stellar interior, the fluid damps via radiative diffusion, exerting effective tidal torques that transfer the angular momentum from the orbit to the stellar spin (Zahn 1975). This classical theory of dynamical tides was originally proposed for massive stars, and several works have calculated the tidal evolution of sdB binaries with this model or its adaptations (Geier et al. 2010; Pablo et al. 2012; Preece et al. 2018).  \nHowever, Ma & Fuller (2023) recently pointed out that an assumption in Zahn's model may not be true for stripped helium-burning stars, like Wolf-Rayet stars in the case of massive stars and sdBs in the case of low-mass stars. While Zahn (1975) assumed that the waves are all efficiently damped when they propagate to the stellar surface, in these stripped stars they may be otherwise reflected and form standing waves. This is particularly true for high-frequency gravity waves excited by short-period orbits, with less efficient radiative damping in stellar envelopes. By direct calculations of tidally excited oscillations with radiative damping, Ma & Fuller (2023) showed that real tidal torques should have more complicated frequency dependence than the simple power-law relation derived from Zahn's model. This new approach brings concerns to the existing predictions of sdB rotation rates based on Zahn's tidal calculation.  \nIn this paper, we build sdB models and calculate their tidal evolution with the method in Ma & Fuller (2023). We carry out direct calculations of stellar oscillations and their tidal response, and we find that standing waves indeed exist in these sdBs. The tidal torques are hence different from Zahn's, and our results for sdB rotation rates are consistent with the observed trends of tidal synchronization for these systems. The manuscript is organized as follows: in Section 2 we describe the physics of sdB spin-up from tidally excited g-mode oscillations; in Section 3.1 we describe our model setup; and in 3.2 and 3.3 we describe how we calculate the oscillation modes and the binary evolution. We show the results for tidal torque calculations in 4.1 and the spin-orbit evolution of sdB binaries in 4.2. We discuss the limitations of our models and the various related physical processes in Section 5. We finally conclude in Section 6.", '2. TIDAL PHYSICS': "For subdwarf B binaries, the tidal dissipation inside the sdB star is thought to be responsible for its tidal evolution. In this picture, the tides are excited by the tidal gravity potential from the companion star, which is usually an M-dwarf (dM) or a white dwarf (WD). When the binary orbit is faster than the stellar spin,  \nthe tidal dissipation transfers energy and angular momentum from the orbit to the star.  \nFor rotating sdB stars with convective cores and radiative envelopes, three possible tidal dissipation mechanisms could be at work: namely, 1) the turbulent viscous dissipation of equilibrium tidal bulges in the stellar core (Eggleton et al. 1998; Eggleton 2006); 2) the turbulent viscous damping of tidally excited inertial waves restored by the Coriolis force (Ogilvie 2013; Wu et al. 2024); and, 3) the radiative damping of tidally excited gravity waves in the envelope (Zahn 1975, 1977). Studies have found that the first scenario is usually inefficient for close-in subdwarf binaries, as the orbital periods might be shorter than the convective turnover time in the core, such that convective viscous dissipation is suppressed (Preece et al. 2019). The second scenario is likely unimportant for the same reason, and the small convective core sizes of sdB stars further suppress inertial wave dissipation. We hence focus on the last scenario to be the dominant process for tidal evolution.  \nWe sketch the physics picture of radiative dissipation of tidally excited gravity waves in Figure 1. Gravity waves, propagating in the radiative envelope of the star, can be tidally excited by the orbit of the companion. In the outer envelope with large thermal diffusion, these waves damp partially or fully by radiative diffusion, releasing their energy and angular momentum, and hence exert a tidal torque on the star. The orbital angular momentum is thus transferred to the stellar spin. In previous studies, radiative dissipation is often assumed to be efficient so that these gravity waves damp completely in the radiative envelope (Zahn 1975), while in principle they might reflect back and instead form standing waves, i.e. oscillation modes (Ma & Fuller 2023). Hence, a realistic estimate of tidal torques requires calculation of individual stellar oscillation modes.  \nFor an aligned and circular orbit, the tidal torque for a tidally excited oscillation mode α is given by (Ma & Fuller 2023):  \nτ α = -mω α γ α q 2 M 1 R 2 1 | W lm Q α | 2 ω 2 f ( ω α -ω f ) 2 + γ 2 α ( R 1 a ) 2( l +1) , (1)  \nwhere ω α and γ α are the mode frequency and growth rate (with γ α < 0 for damped modes, and the corresponding τ α > 0), and ω f = m (Ω orb -Ω spin ) is the tidal forcing frequency (measured in the frame corotating with the sdB), and Ω spin is the sdB's angular rotation frequency. M 1 and R 1 are the mass and radius of the sdB, q = M 2 /M 1 is the mass ratio of the companion to the sdB, 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  \nFigure 1. A sketch of the physics of sdB tidal spin-up. Gravity waves, propagating in the radiative envelope, can be tidally excited by the gravity from the orbiting companion. In the hydrogen-rich outer envelope, the waves damp (either partially or fully) and deposit their angular momentum into the star, transferring angular momentum from the orbit to the stellar spin. The color scale shows the hydrogen fraction in the radiative envelope.  \n<!-- image -->  \nand W lm is an expansion coefficient of the tidal potential. Q α ≡ ( G/R l +3 1 ) ⟨ ξ α |∇ ( r l Y lm ) ⟩ /ω 2 α is the dimensionless overlap integral describing the spatial coupling between the mode and the tidal potential, where ξ α is the displacement vector of mode α , normalized by ⟨ ξ α | ξ α ⟩ ≡ ∫ star ξ ∗ α · ξ α ρdV = M 1 R 2 1 . In practice, Q α is calculated by the relation Q α = -(2 l +1) δ Φ α / (4 πω 2 α R 2 1 ) (Fuller 2017), where δ Φ α is the surface gravity potential perturbation. For weak damping ( γ α < ω α ), the excitation of individual oscillation modes is independent from each other, such that the total tidal torque can be expressed as  \nτ tide = ∑ α τ α . (2)  \nHence, by solving for the internal oscillation modes (with ω α , γ α and Q α ) inside the sdB, we are able to calculate the torque and the angular momentum transfer rate, given a companion mass and orbit. We stress that there are no free parameters in estimating the strength of tidal torques with this method.", '3. METHODS': 'We calculate the tidal evolution of sdB binaries with the similar method developed for Wolf-Rayet-black hole binaries in Ma & Fuller (2023). We first build realistic single evolving sdB models throughout their heliumburning lifetime (Section 3.1). We then solve for stellar oscillations based on these models to estimate the tidal torques (Section 3.2). Finally, we numerically integrate the coupled spin-orbit evolution of sdB binaries with interpolation between the previously calculated sdB models and tidal torques, with different choices of initial  \nbinary parameters, i.e. the initial orbital periods and companion masses (Section 3.3).', '3.1. Stellar Models': 'We build single sdB models with the MESA stellar evolution code (r12778; Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023). We build two sdB models to represent the two types of sdBs formed from progenitors of different masses, summarized in Sections 3.1.1 and 3.1.2. We turn on element diffusion for 1 H, 4 He, 12 C, 14 N, and 16 O in the MESA models, to account for correct treatments of gravitational settling and radiative acceleration for these atoms. The MESA inlists are available on Zenodo under an open-source Creative Commons Attribution license: doi:10.5281/zenodo.13388526.', '3.1.1. 0 . 47 M ⊙ Canonical SdB Model': "This model represents the most abundant sdBs ('canonical' sdBs) that are formed from low-mass ( M ≲ 2 M ⊙ ) main-sequence progenitor stars. When these stars evolve off main-sequence, they start hydrogen shell burning which deposits helium into their helium core, until they reach the tip of the red giant branch (TRGB) when the helium core exceeds 0 . 46 M ⊙ . At this moment, an off-center helium flash is triggered and the helium burning propagates to the center of the core, while the star loses most of its envelope through binary processes (e.g., a common envelope event), leaving a core helium-burning sdB star with a little of its envelope ( ∼ 0 . 01 M ⊙ ) retained. SdB stars formed this way are insensitive to the masses of their progenitors, and have universal masses of ≈ 0 . 47 M ⊙ . We establish  \nthis sdB model by evolving a 1 . 2 M ⊙ star from zero-age main-sequence (ZAMS) to the TRGB, and then apply an artificial stellar wind to remove its envelope, until the off-center helium burning propagates to the stellar center. This happens at the moment when the envelope mass reaches the desired mass (see details in 3.1.3), due to the specific wind scaling factor we chose. The model then becomes a zero-age canonical sdB model and we evolve it until core helium depletion.", '3.1.2. 0 . 37 M ⊙ Low-mass SdB Model': 'For stars of ∼ 2 -3 M ⊙ , they can also ignite core helium burning without forming a fully degenerate core, at helium core masses less than 0 . 46 M ⊙ . Hence, they usually form sdBs of lower masses compared to canonical sdBs. To simulate this scenario, we evolve a 2 . 7 M ⊙ star from ZAMS to the TRGB, and then remove its envelope by a similar artificial wind until its envelope mass reaches the desired mass (see details in 3.1.3). The model then triggers central helium burning as a zero-age sdB star. The sdB model we build this way has a mass of 0 . 37 M ⊙ .', '3.1.3. Envelope Mass Setup': 'SdBs are known to retain a small amount of hydrogen envelope above their helium cores. The amount of hydrogen can be constrained from their spectroscopic properties, and is found to be between 0 . 001 - 0 . 005 M ⊙ (see, e.g., Figure 10 of Kupfer et al. 2015). We hence adjust the artificial winds such that the stellar models start core helium burning (zero-age sdB) when they have 10 -3 M ⊙ hydrogen left. After that moment, we turn off stellar winds as the envelope-stripping phase is considered completed. We note that, real sdB stars can possibly retain more hydrogen than 10 -3 M ⊙ . However, we found many unstable stellar oscillations in sdB models with more massive envelopes, and we are not able to calculate the tidal dissipation of these modes with our current method. We discuss the influence of envelope masses and these unstable modes in more detail in Section 5.1.', '3.1.4. Convective Core Boundary Setup': "The excitation of gravity waves is sensitive to the size of the convective core, which in turn can be sensitive to how its boundary is treated in stellar evolution models. Unlike the standard convective-overshooting paradigm that has been established for main-sequence stars (see, e.g., the MIST project; Choi et al. 2016), overshooting parameters for stars with a helium-burning convective core, like sdBs, are poorly constrained. Nevertheless, asteroseismic measurements of core heliumburning stars suggest the existence of bigger cores com-  \npared to theoretical modeling (Constantino et al. 2015; Bossini et al. 2017; Noll et al. 2024). We hence turn off overshooting for our stellar models in the core heliumburning phase, instead applying the 'predictive mixing' scheme for convection (Paxton et al. 2018). This allows for a steady growth of the convective core during core helium-burning, more consistent with asteroseismic observations than other choices for convective mixing. Furthermore, the predictive mixing scheme helps prevent 'breathing pulses' at late stages of the core heliumburning phase, which may split the convective core to create small radiative zones, in which very high-order gravity waves can be trapped. Breathing pulses have been argued to be numerical artifacts (Bauer & Kupfer 2021), and we aim to avoid the associated difficulties in computing gravity modes.", '3.1.5. Rotation Setup': 'When stars evolve off main-sequence, their core contracts and spins up, while their envelope expands and spins down. The shear created between the core and the envelope could trigger hydrodynamical and magnetohydrodynamical instabilities, which transfer some of the core angular momentum to the envelope, forming slowly rotating stellar cores (Fuller et al. 2019; Tripathi et al. 2024). Asteroseismic measurements of red clump stars have shown that their core rotation periods are typically ∼ 100 days (Mosser et al. 2012). Therefore, if these stellar cores form sdBs, they should also be slowly rotating.  \nWe applied the modified Taylor-Spruit torque as described in Fuller et al. (2019) in our stellar models, and we found that the stellar models at the end of the envelope-stripping phase rotate slowly, with rotation rates insensitive to the initial rotation at ZAMS. The slow rotation rates are consistent with the slow sdB rotation rates measured in wide binaries, where tidal effects are not important (see, e.g., Figure 3). We hence set the sdB models to be nonrotating at the start of their helium burning phase. Since we only compute our spinorbit evolution by post-processing of the stellar models, without actually updating their rotation rates in MESA (see details in 3.3), the single sdB models remain nonrotating throughout their lifetime.', '3.2. Calculation of Oscillation Modes': 'We calculate the internal oscillation modes for the individual snapshots of our sdB models with the GYRE stellar oscillation code (Townsend & Teitler 2013; Townsend et al. 2018; Goldstein & Townsend 2020). We solve for non-adiabatic oscillations which account for radiative damping in the oscillation equations. We use the second order Magnus differential scheme, as it proves to be the most reliable when dealing with low-frequency  \noscillations. We specify our search to l = m = 2 modes as 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.13388526. When solving for modes, we find that the Brunt-Vaisala frequency (Vaisala 1925) profiles solved from MESA are sometimes not consistent with the density and pressure profiles from the same model, which may lead to inaccurate mode solutions. We hence slightly adjust the stellar profiles with the process described in Appendix A for our stellar models. We checked that the change of stellar structure due to this process is negligible.  \nIn principle, we need to sum over all modes to get the total tidal torque through Equation 2. This is not practically possible as there are an infinite number of modes which could be excited at all frequencies. Nevertheless, we note from Equation 1 that for a given tidal forcing frequency ω f , typically only the few nearly resonant modes with ω α close to ω f contribute significant torques. Torques from other non-resonant modes are usually negligible due to the ( ω α -ω f ) 2 term in the denominator of Equation 1. We hence restricted our mode solutions to a finite period range, namely from 0 . 005 days to 0 . 5 days, and we hence found a finite number of modes. We can then calculate the total tidal torques as long as the forcing frequency ω f is between 2 π/ (0 . 5 d) = 12 . 57 d -1 and 2 π/ (0 . 005 d) = 1257 d -1 . The ω f calculated from our spin-orbit evolution usually lies well within this range, except for some systems that reach tidal synchronization, whose ω f should approach zero. We hence stop the evolution when ω f reaches the minimum mode frequency 12 . 57 d -1 . We checked that our binary models reaching this condition are at least at 80% synchronization, so we define all systems with Ω spin ≥ 0 . 8 Ω orb as tidally synchronized.', '3.3. Spin-Orbit Evolution': "We integrate the spin-orbit evolution of the sdB binaries from the stellar models and tidal torques we computed. For simplicity, we ignore the tidal dissipation in the companion star (see discussion in Section 5.7). Assuming circular orbits, the orbital angular momentum of the system is lost due to gravitational wave (GW) radiation and tides, while the sdB receives spin from the tidal torque:  \n˙ J orb = -τ GW -τ tide , (3)  \n˙ J spin = τ tide , (4)  \nwhere τ GW = (32 / 5)( G/a ) 7 / 2 c -5 M 2 1 M 2 2 √ M 1 + M 2 is the effective torque by GW radiation (Peters 1964). The  \nGW orbital decay timescale is then given by  \nT GW = 5 c 5 (1 + q ) 1 / 3 P 8 / 3 orb / (64(4 π 2 ) 4 / 3 G 5 / 3 M 5 / 3 1 q ) ≈ 180 Myr ( P orb / 1 h) 8 / 3 (5)  \nfor an equal mass sdB binary ( q = 1), with a 0 . 46 M ⊙ canonical sdB star, comparable to the sdB lifetime of ∼ 150 Myr for very short-period systems. This means GWorbital decay needs to be included in the spin-orbit evolution.  \nAs we expect efficient angular momentum transport during the core-helium burning phase (Fuller et al. 2019; Fuller & Lu 2022; see discussions in 3.1.5), we assume rigid rotation for the sdB star, with a uniform rotational frequency Ω spin . We discuss the case of differential rotation in Section 5.3. The coupled spin-orbit evolution can then be integrated by:  \n˙ Ω spin = ˙ J spin I spin -Ω spin ˙ I spin I spin , (6)  \n˙ Ω orb = ˙ J orb I orb -Ω orb ˙ I orb I orb = -3 ˙ J orb I orb , (7)  \nwhere I spin is the moment of inertia of the sdB star, I orb = µa 2 is the moment of inertia of the orbit and we made use of Kepler's Third Law to simplify Equation 7. This means the spin of the sdB star may also change due to the changes of its internal structure and hence moment of inertia.  \nWe make use of the integration machinery constructed in Ma & Fuller (2023), with the same interpolation method. As sdB lifetime is typically longer than the Wolf-Rayet stars in Ma & Fuller (2023), we choose the integration timestep to be 0 . 1 times the values derived from the timestep control method described in Ma & Fuller (2023). We integrate the evolution from 1 year after the start of the sdB helium-burning phase, and we stop when the model depletes its core helium (defined by the time when the core helium fraction drops below 1%) or when the system reaches ω f ≲ 12 . 57 d -1 (see Section 3.2). The initial rotational period of sdBs is set to be 60 days, to match the observed values from single sdB stars (Silvotti et al. 2022). While single sdBs may not represent a fair sample of binary sdBs at birth, this assumed initial rotational frequency is very low and never important for systems that reach synchronization. As there are only limited theoretical constraints on the initial binary parameters of sdBs, we vary the companion masses uniformly between 0 . 1 M ⊙ and 0 . 8 M ⊙ , and choose the initial orbital periods to be between 1 to 18 hours, aiming to cover the observed parameter space of close-in sdB binaries (Schaffenroth et al. 2022). We checked our results are robust against different choices of timestep resolution.", '4. RESULTS': 'In this section, we show the results of our tidal torque calculation and spin-orbit evolution.', '4.1. Tidal Torque': "In Figure 2, we show our calculated tidal torque magnitude with its dependence on the period of tidal forcing ( P f ≡ 2 π/ω f ). The calculation is based on the 0 . 47 M ⊙ sdB model with a companion of 0 . 4 M ⊙ , when its central helium fraction drops to 60%, and we assume the sdB is nonrotating such that ω f = m Ω orb . The torque generally has a complicated dependence on the forcing period. By plotting the torque contributions from each individual oscillation mode J α , we see that this dependence is caused by summing over the resonance peaks of many modes with different frequencies. When the tidal forcing frequency gets close to one of the mode frequencies, the ( ω α -ω f ) 2 term in Equation 1 vanishes, and the total torque becomes dominated by the strong resonance peak of that mode. Therefore, the frequencies/periods of these peaks are the frequencies/periods of individual oscillation modes inside the star. These peaks have a nearly uniform period spacing, a feature expected for gravity (g) modes. In the right panel of Figure 2, we show some example eigenfunctions of these modes, and we can see that they are indeed g modes that propagate inside the radiative envelope of the star.  \nPrevious studies involving the tidal dissipation of gravity waves in sdBs usually use Zahn's model for dynamical tides (Geier et al. 2010; Pablo et al. 2012; Preece et al. 2018), which assumes these waves are efficiently damped as they reach the stellar surface ('travelingwave limit'). We see from the right panel of Figure 2 that this is clearly not always the case for individual resolved stellar oscillations. The blue line shows the eigenfunction of an example oscillation at short ( P α ≲ 0 . 14 days) period. We see that instead of efficiently damping near the stellar surface, the wave reflects back at the stellar surface and forms a standing wave with nodes. This means Zahn's picture may overestimate the mode damping rate, hence the tidal dissipation.  \nFor comparison, we show the tidal torques calculated based on Zahn's formalism with a modified formula given by Kushnir et al. (2017):  \nτ Zahn = β 2 GM 2 2 r c ( r c a ) 6 s 8 / 3 c ρ c ¯ ρ c ( 1 -ρ c ¯ ρ c ) 2 , (8)  \nwhere 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. β 2 is a dimensionless coefficient solved from stel-  \ns, and different main-sequence and WolfRayet stellar models have β 2 ≈ 1 (Kushnir et al. 2017). Its dependence on the tidal forcing period is mostly the power-law term in s 8 / 3 c a -6 , as seen in Figure 2. This is clearly different from the resonance peak dependence we find from realistic mode calculations. We see that at short periods, if the binary orbit has an off-resonance tidal forcing period (i.e., not close to any stellar oscillation modes), the real tidal torque can be orders of magnitude lower than Zahn's prediction. On the other hand, if the orbit is on resonance, the torque might be significantly larger than Zahn's prediction.  \nHowever, when the tidal forcing is on resonance with one of the oscillation modes, the mode amplitude becomes so large that it can trigger nonlinear wave dissipation. In this scenario, the oscillation mode excites a number of nearby daughter and granddaughter modes, and the overall damping rate by this sea of coupled modes could be much larger than the radiative damping of individual modes (Barker & Ogilvie 2011; Weinberg et al. 2012). The level of nonlinearity can be estimated by the second-order nonlinear term in the momentum equation ξ · ∇ ξ ∼ ( dξ r /dr ) ξ : when dξ r /dr > 1, nonlinear effects become very strong.  \nWith the above criterion, we showed where tidally excited modes are strongly nonlinear in Figure 2 with dashed grey lines. We see that nonlinear effects mostly affect on-resonance torques at short ( P tide < 0 . 07 days) periods. For these torques, the resonance peaks will be smoothed out when additional nonlinear damping is present, and their actual magnitude may not be as large as seen in Figure 2.  \nAs the star and the orbit evolves, both the oscillation mode period (hence the location of the resonance peaks) and the tidal forcing period change over time, so the system can quickly pass through resonances (as long as resonance locking does not happen, see discussions in Section 5.2). Since the resonances are narrow, the system spends more time out of resonance than in resonance (i.e., with torques much weaker than Zahn's prediction), the accumulated angular momentum received by the sdB star can still be less than the predictions from Zahn's theory.  \nAt longer periods, gravity waves have larger wave numbers, and are expected to damp more efficiently with radiative damping. With larger γ α , the γ 2 α term becomes more important in the denominator of Equation 1, smoothing out the resonance peaks. We see in Figure 2 that this is exactly the case for the tidal torques at long period ( P f > 0 . 14 d), when the individual modes damp so much that the resonance structure gets smoothed out. In addition, the shape of the example eigenfunction (red  \n<!-- image -->  \nFigure 2. Left : The tidal torque calculated for a 0 . 47 M ⊙ sdB model with a companion of 0 . 4 M ⊙ , when the central helium fraction is 60%, assuming the sdB is nonrotating. The thick lines show the total torque calculated by summing over contributions from individual tidally excited g modes (thin lines). We also show the torque calculated from Zahn's formalism for comparison. At short periods, the total tidal torque is dominated by resonance peaks from individual weakly-damped g modes, different from the power-law dependence on forcing period of Zahn's formalism. However, on-resonance torques at short periods (thick dashed grey lines) may trigger nonlinear dissipation, and the real torques may not be as large as seen in the plot. At longer periods, g modes are more efficiently damped and the torque agrees better with Zahn's model. Right: the mode eigenfunctions for an example weakly-damped standing mode (blue line) and an example strongly-damped mode (red line) in the left panel. The weakly-damped mode has nodes in its eigenfunction, while the strongly-damped mode efficiently damps near the surface.  \n<!-- image -->  \nline) shown in the right panel of Figure 2 becomes closer to a traveling wave that efficiently damps near the stellar surface, as Zahn's formalism assumes. The tidal torque's dependence on tidal forcing period also gets closer to Zahn's power-law dependence as expected. This further shows that Zahn's traveling wave picture is a limit case of realistic tidal torques at long periods.", '4.2. Tidal Synchronization': 'With the tidal torques calculated, we are able to integrate the coupled spin-orbit evolution of our models. In Figure 3, we show the calculation for our 0 . 47 M ⊙ canonical sdB model and 0 . 37 M ⊙ low-mass sdB model with a fixed companion mass of 0 . 4 M ⊙ , with different initial orbital periods. The rotational and orbital periods are evaluated at the end of the spin-orbit evolution (defined in Section) 3.3, and some ultra-short-orbit synchronized systems that reaches P orb < 0 . 02 d due to gravitational wave orbital decay are not shown. We see that for the 0 . 47 M ⊙ sdB model, all systems with orbital periods less than ∼ 0 . 2 days reach tidal synchronization, while for the 0 . 37 M ⊙ sdB model, the synchronization period becomes ∼ 0 . 15 days.  \nTo compare with observations, we also plot the measured rotational and orbital periods for short-period sdB binaries in Figure 3. The different colors show  \nthe sdB rotation rates derived from spectral line measurements (HS 0705+6700, Drechsel et al. 2001; CD-30 11223, Vennes et al. 2012; SDSS J162256.66+473051.1, Schaffenroth et al. 2014; PTF1 J0823+0819, Kupfer et al. 2017; PTF1 J011339.09+225739.1, Wolz et al. 2018; ZTF J2130+4420, Kupfer et al. 2020a; ZTF J2055+4651, Kupfer et al. 2020b; SDSS J082053.53+000843.4, Schaffenroth et al. 2021; HW Vir, Esmer et al. 2021; and EPIC 216747137, Silvotti et al. 2021), asteroseismic p-mode frequency splitting (NY Vir, Charpinet et al. 2008; Feige 48, Van Grootel et al. 2008; V1405 Ori, Reed et al. 2020; and HD 265435, Pelisoli et al. 2021), or g-mode frequency splitting (KIC 11179657 and KIC 2991403, Pablo et al. 2012; FBS 1903+432, Telting et al. 2014; KIC 7664467, Baran et al. 2016; EQ Psc and PHL 457, Baran et al. 2019; KIC 2438324, Sanjayan et al. 2022; TYC1 4544-2658-1, Silvotti et al. 2022; and PG 0101+039, Ma et al. 2023), respectively.  \nWe see that all the observed systems with P orb ≲ 0 . 2 d are close to tidal synchronization, while all but two sys-  \nFigure 3. The observed trend of tidal synchronization (defined by 0 . 8 ≤ Ω spin / Ω orb ≤ 1, shaded grey region) for short-period sdB binaries (crosses) versus the calculation results from binary spin-orbit evolution (dots). The red, purple and green crosses indicate rotational measurements from spectral line broadening, p-mode frequency splitting and g-mode frequency splitting. The red and blue dots indicate the modeled sdB binaries with the 0 . 47 M ⊙ and the 0 . 37 M ⊙ sdB primary, respectively. All modeled binaries have a 0 . 4 M ⊙ companion and initial orbital periods ranging from 1 to 18 hours. For 0 . 47 M ⊙ sdB binaries, we see that all systems with orbital periods less than ∼ 0 . 2 days reach tidal synchronization, while for 0 . 37 M ⊙ sdB binaries, this synchronization period becomes ∼ 0 . 15 days due to weaker torques on these smaller sdBs. These results match the observed trends of sdB tidal synchronization. In addition, most systems with 0 . 3 d ≲ P orb ≲ 0 . 6 d are spun-up to rotational periods of a few days, which also agree with the observed period range of partially synchronized binaries.  \n<!-- image -->  \nms 1 above this period are not synchronized. This matches strikingly well with the theoretical prediction from our sdB models. In addition, the models with  \n0 . 3 d ≲ P orb ≲ 0 . 6 d are tidally spun-up to a rotational period of a few days, also consistent with the observed partially-synchronized systems in that period range. Hence, our theoretical calculations agree with the observation data.  \nNote that in Figure 3 the spin and orbital frequencies are shown at the end of the spin-orbit evolution, while measurements for realistic systems usually occur when the sdB still undergoes core-helium burning. Therefore, it is illustrative to show the tidal synchronization  \nFigure 4. The ratio between the tidal synchronization timescale ¯ T sync and the sdB lifetime T EHB interpolated between different choices of companion masses and initial orbital periods. The blue regions show the parameter space where ¯ T sync < T EHB . The brackets on the y axis indicate the typical companion masses for sdB+WD or sdB+dM systems (Schaffenroth et al. 2022). Left : The results for the 0 . 47 M ⊙ sdB model, with T EHB = 164Myrs. We see that for systems in orbits less than ∼ 0 . 15 -0 . 22 d, the synchronization timescale is less than the stellar lifetime, meaning these systems are likely to be observed as tidally synchronized. The results have weak dependence on companion masses. Right : The results for the 0 . 37 M ⊙ sdB model. The critical orbital period below which systems become synchronized now becomes ∼ 0 . 10 -0 . 17 d.  \n<!-- image -->  \ntimescales calculated for our models, and to compare them with the sdB lifetime T EHB . If the synchronization time is shorter, then we expect those systems are likely to reach tidal synchronization to be observed. As we consider systems with 80% synchronization as synchronized (see discussions in Section 3.2), we define the tidal synchronization timescale throughout the whole evolution as:  \n¯ T sync ≡        t Ω spin / Ω orb =0 . 8 , if synchronized 0 . 8 T EHB ( Ω spin Ω orb ) -1 final , if not (9)  \nwhere t is the stellar age since the start of sdB core helium-burning. We then run a grid of spin-orbit evolution for both of our 0 . 47 M ⊙ and 0 . 37 M ⊙ sdB models, with initial orbital periods of (1 , 2 , 3 , 4 , 5 , 6 , 7 , 8) hours and companion masses of (0 . 1 , 0 . 2 , 0 . 3 , 0 . 4 , 0 . 5 , 0 . 6 , 0 . 7 , 0 . 8) solar-masses, and calculated their tidal synchronization timescale defined by Equation 9.  \nWe show the interpolated results for the ratio ¯ T sync /T EHB based on this grid in the whole parameter space in Figure 4. On the y -axis we also label the typical companion masses for sdB+dM and sdB+WD binaries (Schaffenroth et al. 2022). We see that this ratio ranges from 10 -1 . 5 to 10 1 . 5 , with weak dependence on the masses of the companion star. This is an expected result from our tidal torque formula (Equation 1): the torque  \ndepends on the mass-ratio (hence the secondary mass) as τ α ∝ q 2 a -6 , where a = ( G ( M 1 + M 2 ) / Ω 2 orb ) 1 / 3 ∝ (1 + q ) 1 / 3 , hence τ α ∝ q 2 / (1 + q ) 2 . For typical sdB binaries, q ranges from 0.3 to 1.5, and the corresponding torque scaling is maximally different by only a factor of ∼ 7. In contrast, ¯ T sync /T EHB varies by a factor of ∼ 10 3 over the period range shown in Figure 4, due to its strong dependence on semi-major axis.  \nThe blue colored regions in Figure 4 show the parameter space where ¯ T sync /T EHB < 1, or where the binaries are expected to be synchronized. We see that for 0 . 47 M ⊙ sdB binaries, systems with initial orbital periods less than ∼ 0 . 15 -0 . 22 days have synchronization timescales shorter than the 164 Myrs sdB lifetime, while for 0 . 37 M ⊙ sdB this period becomes ∼ 0 . 10 -0 . 17 days for its 444 Myrs lifetime, depending on the companion masses. As binaries at these orbital periods produce weak GW emission, their orbital periods are nearly constant, and these results confirm the critical orbital periods for synchronization shown in Figure 3.  \nWe note from Figure 4 that, even though systems below the synchronization periods can reach tidal synchronization in the sdB lifetime, in most of the parameter space their synchronization timescales are not less than the corresponding sdB lifetime by one order of magnitude. This is especially true for sdB+dM binaries with M companion ≲ 0 . 3 M ⊙ , and we can see that ¯ T sync < 0 . 1 T EHB is only achieved for those binaries in P orb ≲ 0 . 05 d ≈ 1 hour orbits. This is consistent with  \nthe findings that sdB binaries with small companions (dMs or brown dwarfs) can be slightly sub-synchronized even at orbital periods less than ∼ 2 . 5 hours (e.g., SDSS J162256.66+473051.1, a 64% synchronized system with P orb = 1 . 67 h, Schaffenroth et al. 2014; and SDSS J082053.53+000843.4, a 65% synchronized system with P orb = 2 . 3 h, Schaffenroth et al. 2021).  \nSchaffenroth et al. (2021) further point out that synchronized binaries locate further away from the zero-age extreme horizontal branch (ZAEHB) on the log g -T eff diagram compared to these sub-synchronized systems, suggesting that those synchronized binaries might be older. Our findings that the tidal synchronization timescales at small orbital periods are shorter than the sdB lifetime (but not by orders of magnitude), agrees with this explanation.  \nHistorically, the orbital inclinations and companion masses are hard to acquire for non-eclipsing sdB binaries. Some works hence assume tidal synchronization for short-period binaries, and derive the orbital parameters from the orbital periods by setting P orb = P rot (e.g., Geier et al. 2010). However, if we can assume tidal synchronization for systems with ¯ T sync < 0 . 1 T EHB , we see that this method should only apply to binaries with P orb ≲ 1 h. This is much shorter than the synchronization period ( P orb , sync = 1 . 2 d) assumed in Geier et al. (2010), meaning that in their work the companion masses/inclinations might be over/underestimated.  \nAdditionally, binaries with P orb ≲ 1 h can undergo significant orbital decay due to gravitational-wave radiation, and it is questionable whether these systems can ever reach 100% tidal synchronization, as tides at subsynchronization may not be strong enough for Ω spin to fully catch up with Ω orb (as in the case of WD binaries, see, e.g., Scherbak & Fuller 2024). Mass-transfer may also happen for these binaries, making their evolution more complicated (Bauer & Kupfer 2021).', '5.1. Tidal Torque Scaling with Stellar Radii': "We saw in Section 4.2 that the binary orbital period required to reach tidal synchronization is shorter for the 0 . 37 M ⊙ low-mass sdB, compared to the 0 . 47 M ⊙ canonical sdB. This means the tidal torque must be weaker for low-mass sdBs. To explain the reason, we consider an equal-mass binary ( q = 1), and rewrite the mode torque of equation 1 as  \nτ α = -f ( ω f ) γ α J α C α S ( R 1 ) , (10)  \nwhere  \nis a dimensionless function describing the resonance dependence of the torque on ω f , whose scaling should be of similar order for different stellar models. The quantity  \nJ α ≡ mω α ⟨ ξ α | ξ α ⟩ (12)  \nis the angular momentum of the oscillation mode α with azimuthal wave number m , which roughly scales as J α ∝ M 1 R 2 1 for stars of similar structure at the same mode frequency. For dissipating modes, the rate for mode α to dissipate its angular momentum (i.e., exerting a torque) is γ α J α . The quantity  \nC α ≡ | W lm Q α | 2 ⟨ ξ α | ξ α ⟩ /M 1 R 2 1 (13)  \ndescribes the dimensionless coupling of the tidal potential and the oscillation mode, and S ( R 1 ) ≡ ( R 1 /a ) 2( l +1) is the scaling of the tidal forcing strength.  \nWith this notation, we can see that the quantity γ α J α C α represents the rate at which mode α deposits its angular momentum into the star (i.e., the tidal torque), per tidal forcing strength. We plot this quantity for the 0 . 47 M ⊙ and 0 . 37 M ⊙ sdB models in Figure 5, and we see that they have similar orders of magnitude without clear dependence on the different stellar models. This is expected as these sdB stars have very similar internal structures.  \nTherefore, the main scaling of the physical tidal torque comes from the forcing strength S ( R 1 ). As shown in right panel of Figure 5, the 0 . 37 M ⊙ sdB is more compact than the 0 . 47 M ⊙ sdB due to its lower mass, and its radius R 1 is smaller by a factor of ∼ 2. For l = 2 modes, this produces a difference in the tidal torque by ( R 0 . 37 /R 0 . 47 ) 6 ∼ 1 / 64. Since the 0 . 37 M ⊙ sdB has a smaller moment of inertia by a factor of a few, its synchronization time scale is about ten times longer. We can confirm this result from Figure 4: the tidal synchronization timescale at P ∼ 0 . 2 d for 0 . 47 M ⊙ sdB binaries is T EHB , 0 . 47 = 164 Myrs, roughly 10 times shorter than the timescale of 0 . 37 M ⊙ sdB binaries at the same period (a few times T EHB , 0 . 37 = 444 Myrs).  \nThe above analysis can also provide insight into the tidal torques for sdBs with different hydrogen envelope masses. Detailed sdB modeling has shown that hydrogen envelope masses range from 0 . 001 to 0 . 005 M ⊙ (Kupfer et al. 2015). This small amount of hydrogen never affects the core structure of the helium-burning sdBs, but it can greatly change the stellar radius.  \nIn the right panel of Figure 5 we show the stellar radii for some 0 . 47 M ⊙ sdB models that retained more hydrogen than 10 -3 M ⊙ . Compared to the original 10 -3 M ⊙ hydrogen model, we can see that even a slight increase of hydrogen could increase the sdB radius by a factor  \nf ( ω f ) ≡ ω 2 f ( ω α -ω f ) 2 + γ 2 α (11)  \n<!-- image -->  \nFigure 5. Left : The magnitude of γ α J α C α (defined in the main text) for different sdB models. This quantity determines the tidal torque of each mode per unit tidal forcing strength, and we see that they are of similar orders of magnitude without clear dependence on the different stellar models. Hence, the physical torque should be roughly proportional to the tidal forcing strength which scales as R 6 . Note that sdB models with hydrogen masses greater than 10 -3 M ⊙ have some unstable modes, and their tidal excitation can not be treated with our current method. Right : Stellar radii as a function of age for different sdB models. More massive systems with more hydrogen left in the envelope have larger radii, and the tidal torques are expected to be stronger based on the R 6 scaling.  \n<!-- image -->  \nof ∼ 1 . 5 -2. We further plot the γ α J α C α calculated for these models in the left panel of Figure 5, and we see that despite some scatter, they are similar to the 10 -3 M ⊙ hydrogen model. We hence expect the τ ∝ R 6 scaling roughly holds for these models, and the tidal torque for these larger sdBs could be larger by a factor of ∼ 10 -50, which will increase the synchronization transitional period.  \nWe note, however, there is a reason that we did not actually compute the tidal torques for these more extensive sdB models. We see in Figure 5 that for sdB models with M H > 10 -3 M ⊙ , there exists a period range where the mode growth rate γ α is positive, i.e. where the modes are unstable . This is caused by the so-called κ -mechanism in these stars, where the partial ionization of iron creates an opacity bump, generating self-excited oscillations (Charpinet et al. 1996, 1997). Our torque in equation 1 only holds for damped oscillation modes, and it is unclear how these self-excited oscillations would interact with tidal forcing (see, e.g., Fuller 2021). These unstable modes have periods of ∼ 0 . 05 -0 . 1 days, so they could be very important for the tidal evolution of sdBs in ∼ 0 . 1 -0 . 2 day orbits. Future works should investigate how these modes will behave under tidal excitation.  \nThe above analysis also explains the discrepancies between the tidal synchronization period we calculated and those estimated by Preece et al. (2018), who applied  \nZahn's traveling wave limit. The synchronization periods we found for the 0 . 47 M ⊙ canonical sdB model are longer than the periods from their work, which means the tidal torque in our cases is stronger. This might be because Preece et al. (2018) used an sdB model with only 10 -4 M ⊙ hydrogen left, whose radius is smaller than our model. Hence, even though the gravity waves are more damped with Zahn's traveling wave limit in their models, the tidal torque can still be weaker due to its strong dependence on the stellar radius.", '5.2. Resonance Locking': 'In binary systems, if the tidal torque consists of many resonance peaks from individual modes, a process called resonance locking may occur (Witte & Savonije 1999, 2001). In this scenario, the forcing frequency of the binary enters a resonance with one of the oscillation modes α , and stays as  \nω f ≡ m (Ω orb -Ω spin ) ≃ ω α (14)  \nthroughout the binary lifetime. As ω α evolves on its own timescale which is independent of the binary separation, this scenario may result in very different binary evolution history compared to other tidal theories.  \nTo see whether resonance locking can happen for sdB binaries, we write the evolution of forcing frequency as  \n˙ ω f = m ( 3( τ tide + τ GW ) I orb -τ tide I spin ) , (15)  \nwhere we substitute Equations 3, 4, 6 and 7, and neglect the ˙ I spin term as the stellar structure and moment of inertia only changes slowly during the evolution. To maintain a resonance lock, we must have ˙ ω f = ˙ ω α . For helium-burning subdwarfs, their g-mode frequency increases over time, hence the necessary (but not sufficient) condition for resonance locking to occur is ˙ ω f > 0, or  \n1 + τ GW τ tide > I orb 3 I spin . (16)  \nFor binaries of order-of-unity mass ratios, I orb = µa 2 ∼ M 1 a 2 ≫ M 1 R 2 1 > I spin . The above relation hence never holds for realistic sdB binaries, unless τ GW ≫ τ tide . This can happen either for already close-to-synchronization binaries or for wide binaries, where τ tide becomes very small in both cases. In the former case, the low-frequency oscillating g-modes that contribute most to the tidal torque should be very efficiently damped (see Section 4.1), which prevents resonance peaks from forming. In the latter case, tidal evolution is not important, because it would occur on a gravitational wave inspiral time which is very long for wide binaries. We hence do not expect resonance locking to happen for sdB binaries, which is confirmed with our numerical spin-orbit evolution calculations.', '5.3. Differential Rotation': 'As we expect very efficient angular momentum transport inside the sdBs, we assume they are rigidly rotating in our spin-orbit evolution calculations. Observationally, asteroseismology can measure the internal rotation of stars via frequency splittings of g-mode and p-mode oscillations (Aerts et al. 2010). Since g-modes mainly probe the deeper region of the star, while p-modes probe the outer layers, a difference between the rotational rates derived from their frequency splittings may suggest the level of differential rotation between the stellar core and the outer layers.  \nWith Kepler /K2, there have now been a handful of pulsating sdBs with both p-mode and g-mode frequency splitting measured. Kern et al. (2018) reports that for the sdB+WD system KIC 11558725, the rotational rate derived from p-mode splitting is P p = 40 . 2 ± 0 . 3 days, while the rate from g-mode splitting is found to be P g = 45 . 1 ± 7 . 8 days, showing that KIC 11559725 is roughly rigidly-rotating. Similar results are found for sdBs EPIC 220422705 (Ma et al. 2022; where P p ∼ 29 days and P g ∼ 32 days), and PG 0101+039 (Ma et al. 2023; where P p = 8 . 60 ± 0 . 16 days and P g = 8 . 81 ± 0 . 06 days). We note that the slightly faster rotation rates measured from pmode splitting for these systems are consistent with our tidal spin-up picture with dissipating gravity waves: as  \ngravity waves mostly dissipate in the outer envelopes of the star (Figure 1), they exert local tidal force mostly in these regions, and angular momentum is subsequently transported inward.  \nHowever, the above systems all have slow rotation rates compared to the typical rotational periods ( P rot ≲ 0 . 3 d) of tidal synchronization calculated from our models. This means the tidal torque are weak for these systems. We pay particular attention to two systems: 1) V1405 Ori (also EPIC 246683636 or KUV 04421+1416), an sdB+dM binary with an orbital period of 0.498 days (Reed et al. 2010). Reed et al. (2020) obtained the p-mode and g-mode splitting from K2 observations, and determined a p-mode derived rotation rate of 0 . 555 ± 0 . 029 days and (marginally) a g-mode derived rotation rate of 4 . 2 ± 0 . 4 days. If the P g measurements are reliable, this appears to be a differentiallyrotating sdB whose envelope is almost synchronized while the interior is not; 2) TIC 441725813 (also TYC 4427-1021-1, 2MASS J17045838+7304433, or Gaia DR3 1655107708129775744), a system with recent asteroseismic measurements from TESS. Su et al. (2024) derived a P p = 17 . 9 ± 0 . 7 days and P g = 85 . 3 ± 3 . 6 days from the data, and their initial analysis showed that it is possibly in a binary orbit of 0 . 28 days, slightly longer than the synchronization orbital period we found in theories. If the orbit is confirmed, this will be a sub-synchronized sdB star with substantial differential rotation.  \nThe above systems suggest that for strong tidal torques acting in the sdB outer envelope, the assumption of rigid rotation might break down. As the stellar core, which carries most of the stellar mass and moment of inertia, is weakly coupled to the envelope in this case, our calculated tidal synchronization orbital period might be shorter than that needed to synchronize the envelope.  \nWhile the current sample is limited, we expect further observations with TESS can provide us more sdB pulsators with differential rotation measured (Baran et al. 2023, 2024; Uzundag et al. 2024). It is beyond the scope of this work to develop methods for the spin-orbit evolution of differentially rotating sdB models, but we comment that it might be a crucial factor to understand sdB tidal spin-up.', '5.4. Implications for Rotation Periods of CO WDs': 'Rotating sdBs can potentially form rapidly rotating carbon-oxygen (CO) white dwarfs after their nuclear burning stops. If the spin angular momentum is conserved after the core-helium-burning phase, the rotation periods of CO WDs are then given by:  \nP rot , WD = I WD I sdB , end P rot , end , (17)  \nwhere I sdB , end and P rot , end are the moment of inertia and the rotation period of the sdB at the end of its helium-burning phase (defined by the time when the central helium fraction drops below 1%), and I WD is the moment of inertia of the CO WD. Evolving our 0 . 47 M ⊙ sdB model until it forms a CO WD, we find that the stellar moment of inertia decreases from I sdB , end = 3 × 10 51 g cm 2 to I WD = 1 . 8 × 10 50 g cm 2 in 50 Myrs, due to the shrinking of the star after nuclear burning stops.  \nSince tidal synchronization occurs at orbital periods less than ≈ 0 . 2 d, rapidly rotating CO WDs formed from tidally synchronized sdBs should have P rot , WD ≲ (1 . 8 × 10 50 / 3 × 10 51 ) × 0 . 2 d ≈ 17 min. At orbital periods P orb ≲ 1 h, the gravitational-wave decay timescale becomes less than the sdB lifetime, so the sdB does not form a WD before mass transfer with its companion. This gives a lower-limit of P rot , WD ≳ (1 . 8 × 10 50 / 3 × 10 51 ) × 1 h ≈ 4 min if mass transfer has not occurred. We hence crudely estimate that rapidly rotating CO WDs formed from synchronized sdB stars can have rotation periods between 4 to 17 minutes. This corresponds to rotation rates roughly a hundred times larger than ordinary WDs.  \nIn the above analysis, we ignored any tidal torques after the core-helium-burning phase. Since the star can rotate much faster than the orbit as it contracts, tidal dissipation may spin it back down, producing longer rotation periods than those listed above. Future works should investigate this scenario to have more realistic estimates of the rotation rates of CO WDs originating from sdB binaries. Nonetheless, future observations of rapidly rotating CO WDs in close binaries may indicate that they were tidally spun up during an sdB phase of evolution.', '5.5. Limitations with the Mode Decomposition Method': "Throughout this work, we calculate the total tidal torque by expanding it into a summation of tidal torques from individual oscillation modes (Equation 2). Townsend & Sun (2023) recently pointed out a potential issue with this 'mode decomposition' method. They found that the magnitude of off-resonance torques could be different from the magnitude calculated from direct solving of fluid equations under the same tidal potential. Dewberry & Wu (2024) introduced a correction to the method for viscous damping addressing this issue, yet it cannot directly apply to the radiative dissipation of tidally excited g-modes.  \nNevertheless, when tidal torques are dominated by resonant modes, the correction to the tidal torque magnitude is likely only significant when the tidal forcing fre-  \nquency is off-resonance (Townsend & Sun 2023). Therefore, our results should still be valid as long as most tidal-spin up is caused by on-resonance torques. Future works should improve our method to investigate whether off-resonance corrections change our results.  \nWe did not include nonlinear dissipation of modes in our calculations (see discussions in Section 4.1). There are no theoretical works on estimating this dissipation on core helium burning stars (except for some toy models, see, e.g., Ma & Fuller 2023). Nevertheless, as shown in Figure 2, strong nonlinear effects should only affect sdBs in the shortest-period orbits which would be tidally synchronized quickly. We therefore suspect that nonlinear effects do not greatly change the synchronization periods we calculated, which matches the observed trends fairly well. Future works should look into the nonlinear effects and their potential influence on sdB tidal spin-up.", '5.6. Tidal Heating': 'Tides not only transfer angular momentum between the star and the orbit, but can also dissipate energy through tidal heating, an effect we ignored in our modeling. The rate of tidal heating can be estimated as:  \n˙ E ∼ ˙ J (Ω orb -Ω spin ) = τ tide (Ω orb -Ω spin ) . (18)  \nFor close sdB binaries, Ω orb ∼ 2 π/ 0 . 2 days ≈ 0 . 4 mHz, and the strongest tidal heating occurs when Ω spin = 0. Therefore, off-resonance tidal torques between 10 32 -10 35 g cm 2 s -2 (see Figure 2) can only produce a tidal heating luminosity of 10 -5 -10 -2 L ⊙ , much less than the nuclear luminosity of a helium-core-burning sdB (typically ∼ 10 L ⊙ ). While on-resonance torques could be 2 -3 orders of magnitude larger, nonlinear dissipation likely makes the actual torque much weaker. Hence we do not expect tidal heating to ever surpass the stellar luminosity.  \nHowever, we note that tides do not dissipate energy uniformly inside the star. For gravity waves with radiative damping, the majority of tidal heating should occur near the stellar surface. Even a small amount of heat may change the stellar structure near the surface, which could in turn affect the dissipation processes in that region. Future works should investigate whether this scenario is important for the tidal evolution of sdB binaries.', '5.7. Other Caveats': 'There are some other caveats with our methods. We calculate the oscillation modes based on pre-calculated nonrotating stellar models. As the star becomes spun up, the modes may start to behave differently. This matters the most when Ω spin becomes comparable with  \nthe Brunt Vaisala frequency N c near the convective core boundary, where the oscillations are couple with the tidal potential. Rotation is also important when Ω spin becomes comparable with ω α , where rotation effects on mode eigenfunctions cannot be ignored.  \nWe checked that the former case is never important for sdB binaries, as N c is typically greater than 10 mHz, much higher than the rotation rates of even synchronized sdBs. For the latter case, as the mode contributing mostly to the tidal torque is the one with ω α ≈ ω f ≡ m (Ω orb -Ω spin ), the rotational effects on modes could only be significant when ω α ∼ Ω orb ∼ Ω spin , i.e., when the system is close to synchronization. However, the tidal torques become small as synchronization is approached, so we do not expect these rotational effects to change our results much.  \nWe ignored the tidal dissipation in the companion star in our spin-orbit evolution. This is valid for WD companions whose dissipation is negligible compared to sdB stars. However, for main sequence companions, tidal dissipation in the companion can occur. Since the orbital moment of inertia of the binary is always much greater than the spin moment of inertia of the companion, this process synchronizes the companion star without changing the orbital period significantly. The dissipation within the companion can thus be ignored if the mass and angular momentum of the binary are conserved.  \nAn exception occurs if the companion star is undergoing magnetic braking. This may remove a significant amount of orbital angular momentum throughout the evolution. The coupled evolution of tides and magnetic braking have only been computed recently in some cases (see, e.g., Sun et al. 2024). Future works should investigate whether this is an important process with more detailed modeling of evolution including magnetic braking.  \nWe set the initial orbital periods of our binary systems as a unified P orb = 60days, while in reality close-in sdB binaries can be born from two separate channels, namely a prior mass-transfer phase, or a common-envelope ejection. In general, sdB binaries born from these channels could have different initial spin periods, but the detailed outcome of these processes is highly uncertain. Nevertheless, we comment that, as long as tides are responsible for most of the sdB angular momentum for synchronized systems, our calculations should be insensitive to the initial orbital setup.  \nWe did not include mass-transfer in the sdB evolution phase, even though it could happen for sdB binaries born at orbital periods less than 2 -3 hours (Bauer & Kupfer 2021). Mass-transfer may remove the sdB outer  \nhydrogen envelope, and hence change the strength of the tides (see discussion in Section 3.1.3). However, we suspect that tidal synchronization will remain efficient for these stars since they fill their Roche lobes.', '6. CONCLUSION': "In this manuscript, we investigated the tidal spin-up of close-in subdwarf B (sdB) binaries. We considered the dissipation of tidally excited gravity waves in the envelopes of sdB stars, and calculated the tidal torques by directly computing the amplitudes of tidally driven oscillation modes in sdB stellar models. We integrated the coupled spin-orbit evolution of these binaries and calculated the resulting sdB rotation rates.  \nWe showed that in contrast to the usual assumption that gravity waves are efficiently damped near the surface ('Zahn's traveling wave limit'), these waves can actually be less damped, and can reflect back to form standing waves in the radiative envelope of sdB stars. The resulting tidal torque is then significantly less than Zahn's theory predicted, and has a complicated resonant dependence on the frequency of the tidal force. At longer periods, the waves are more highly damped and the tidal torque approaches Zahn's limit.  \nFor binaries containing a 0 . 47 M ⊙ canonical sdB, our models predict the system will be tidally synchronized if the orbit is less than ∼ 0 . 2 days. For those with a 0 . 37 M ⊙ low-mass sdB, this tidal synchronization period becomes ∼ 0 . 15 days. These values are very similar to the observed spin rates of sdB binaries (Figure 3), which are tidally synchronized at orbital periods less than ∼ 0 . 2 days. The tidal synchronization timescale has weak dependence on the companion star mass, and is mostly determined by the orbital period.  \nWe investigated how the amount of hydrogen in the sdB envelope could affect the strength of the tidal torque. Since sdBs with more hydrogen have larger radii, and the tidal torque magnitude could scale with the stellar radius as τ ∝ R 6 , tidal torques may be stronger for stars with more hydrogen. However, the existence of unstable oscillations for our sdB models with thicker hydrogen envelopes complicate the calculation of tidal torques.  \nWhen tidally synchronized sdBs evolve into carbonoxygen white dwarfs, we estimate their rotation periods to be between 4 to 17 minutes (if tidal effects after the core-helium-burning phase can be neglected), which corresponds to spin rates roughy a hundred times faster than typical white dwarfs. We pointed out that resonance locking cannot happen in the tidal spin-up phase of sdB binaries, and discussed the limitations of our mode decomposition method to calculate tidal  \ntorques. Differential rotation, rotational effects on oscillations, and tidal heating may also be important. Future works should investigate the above scenarios, and compare them to the growing numbers of rotation rate measurements for sdBs in close binaries.  \nThe agreement between our models and measurements for sdB binaries is very encouraging for the prospect of reliable tidal synchronization predictions. In particular, we expect the physics of tidal spin-up in sdBs to be very similar to that of more massive helium stars in close binaries (Ma & Fuller 2023), which are progenitors of gravitational wave sources and exotic supernovae. We believe the results of this paper increase the credibility of predictions for black hole spins presented in that work.", 'ACKNOWLEDGEMENTS': 'We thank Holly Preece and the anonymous referee for constructive comments on the manuscript. We thank Evan Bauer, Ylva Gotberg, Rich Townsend, Emily Hu, Peter Scherbak, Reed Essick, and Yuri Levin for helpful discussions. This work is partially supported by NASA through grant 20-XRP20 2-0147. L.M. is thankful for the Max Planck Institute for Astrophysics and the Kavli Foundation, who supported the 2023 Kavli Summer Program where many of the above discussions took place.  \nSoftware: MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019), GYRE (Townsend & Teitler 2013; Townsend et al. 2018; Goldstein & Townsend 2020) , SciPy (Virtanen et al. 2020), matplotlib (Hunter 2007), numpy (Harris et al. 2020)', 'A. FIXING MESA PROFILES': 'To calculate the tidal response, we solve for internal stellar oscillations with the GYRE stellar oscillation code (Townsend & Teitler 2013; Townsend et al. 2018; Goldstein & Townsend 2020). The code reads stellar snapshots from MESA as unperturbed background profiles of density, pressure, etc., and then solves the linear perturbation equations of stellar oscillations. The oscillation equations GYRE aims to solve are simplified by assuming hydrostatic equilibrium and mass conservation of the background stellar profile, hence the MESA snapshots provided to GYRE should satisfy the following equations:  \ndP dr = -ρg , (A1)  \ndM r dr = 4 πr 2 ρ, (A2)  \nwhere P, ρ, M r and g = GM r /r 2 are the pressure, density, enclosed mass and gravity inside the star. Further, the Brunt-Vaisala frequency profile inside the stellar model should also satisfy the following equation by definition:  \nN 2 ≡ g ( 1 Γ 1 d ln P dr -d ln ρ dr ) = g ( -g c 2 s -d ln ρ dr ) , (A3)  \nwhere we made use of c 2 s ≡ Γ 1 P/ρ . We rewrite the above equation into the following form:  \nd ln ρ dr = -N 2 g -g c 2 s , (A4)  \nand plot the ratios between the LHS and RHS of Equations A1, A2 and A4 for one of our MESA stellar snapshots in the upper panels of Figure 6. While they should all be unity, we notice that in the original MESA profile, the ratio between the LHS and RHS for equation A4 departs significantly from 1 at the density discontinuity near the convective core boundary at 0 . 025 R ⊙ . In the radiative envelope, this ratio also departs from unity by a few percent at some radii. The GYRE oscillation solutions solved by assuming Equation A4 are hence problematic. In practice, we find that this inconsistency often causes the mode solutions to change drastically between successive MESA snapshots, while in principle we expect them to vary gradually as the modes evolve.  \nWe hence need to fix the stellar profiles provided by MESA to get correct oscillation solutions. As g mode oscillations are most sensitive to the Brunt-Vaisala frequency profile, we aim to keep the value of N as output by MESA, and adjust the density and pressure profiles to satisfy Equations A1, A2 and A4. Hence, we rewrite these equations into  \n2  \nFigure 6. Upper : The ratios between the LHS and RHS of Equations A1, A2 and A4 for an example MESA snapshot, which should all be unity. We see that in the original stellar profile calculated from MESA, Equation A4 is sometimes not satisfied, with the largest discrepancy happening at the convective core boundary with a density discontinuity. We fix the P, M r and ρ profiles by the method described in Appendix A, and the fixed profiles satisfy Equation A4 better. Lower : The original and fixed P, M r , ρ and Γ 1 profiles of the example MESA snapshot. The pressure, density and enclosed mass are almost identical to original MESA profiles, so the change of stellar structure after the fixing process is negligible. Γ 1 is different by a small amount in some regions of the star, which means that the energy processes in the fixed stellar model are in general not consistent.  \n<!-- image -->  \nthe following matrix form:  \nd dr    P ln ρ M r    =    -GM r ρ r 2 -N 2 r 2 GM r -GM r r 2 c 2 s 4 πr 2 ρ    . (A5)  \nThe equations then become a first order ordinary differential equation of the form d y /dr = f ( y , r ; N 2 , c 2 s ) where y ( r ) = [ P ( r ) , ln ρ ( r ) , M r ( r )] is an unknown function to be solved numerically, and N 2 and c 2 s are the ODE parameters that can be fitted from the original stellar profile. To solve for y , we need a set of boundary conditions for P, ρ and M r , which are given by the following physical limits:  \nP ( R star ) = P original ( R star ) ; (A6)  \nρ ( ϵ ) = ρ original (0) ; (A7)  \nM r ( ϵ ) = 4 πρ original (0) 3 ϵ 3 , (A8)  \nHere, ϵ is the radial coordinate of a point very near the center of the star, just above r = 0 to avoid the singularity at r = 0 which can cause numerical problems. Once the solutions are found, the values at r = 0 are acquired by P (0) = P ( ϵ ), ρ (0) = ρ ( ϵ ) and M r (0) = 0. We choose ϵ to be the spatial coordinate of the innermost grid in the MESA model, and we can then solve Equation A5 with the boundary conditions described by Equations A6, A7 and A8.  \nWe use the integrate.solve bvp function in the SciPy python package (Virtanen et al. 2020) to solve for the fixed P, M r and ρ profiles, using the original MESA profiles as our initial guess. We show the comparison between our fixed and original MESA profiles in Figure 6. The fixed P, M r and ρ profiles are almost identical to the original density profile, but they more accurately satisfy Equation A4 in the radiative envelope. While there is still some inconsistency at the convective core boundary due to the density discontinuity, it is likely unimportant as gravity waves are evanescent in the convective core.  \nWe note that, as the adiabatic index Γ 1 ≡ ρc 2 s /P is now calculated from the fixed P and ρ and the original sound speed profile, it can be different from the original Γ 1 by a small amount inside the star. Physically, this means the gas  \nequation of state is now slightly artificial, and future works should investigate a more self-consistent way to deal with this problem. Nevertheless, we find that, after fixing the MESA profiles, GYRE is able to get oscillation solutions that vary gradually across nearby MESA snapshots.', 'REFERENCES': '-. 2001, A&A, 366, 840, doi: 10.1051/0004-6361:20000245 Wolz, M., Kupfer, T., Drechsel, H., et al. 2018, Open Astronomy, 27, 80, doi: 10.1515/astro-2018-0011 Wu, S. C., Dewberry, J. W., & Fuller, J. 2024, ApJ, 963, 34, doi: 10.3847/1538-4357/ad1e54  \nZahn, J. P. 1975, A&A, 41, 329 -. 1977, A&A, 57, 383 Zhang, X., Chen, X., & Han, Z. 2009, A&A, 504, L13, doi: 10.1051/0004-6361/200912336'}