Spin Evolution of Stellar-mass Black Holes Embedded in AGN Disks: Orbital Eccentricity Produces Retrograde Circumstellar Flows

The spin evolution of stellar-mass black holes (sBHs) embedded in AGN accretion disks is an important process relevant to the production of gravitational waves from binary BH (BBH) merger events through the AGN channel. Because embedded sBHs are surrounded by circumstellar disks (CSDs), the rotation of CSD gas flows determines the direction of the angular momentum it accretes. In this Letter, we use global 2D hydrodynamic simulations to show that while a disk-embedded sBH on a circular orbit transforms the initial retrograde Keplerian shear of the background accretion disk into a prograde CSD flow, as in the classical picture of companion-disk interaction theory, moderate orbital eccentricity could disrupt the steady-state tidal perturbation and preserve a retrograde CSD flow around the sBH. This switch in CSD orientation occurs at a transition eccentricity that scales nearly proportional to the local sound speed. This bifurcation in the CSD flow and thereafter spin-up direction of SBHs leads to the formation of a population of nearly antialigned sBHs and should be incorporated in future population models of sBH and BBH evolutions.

The effective spin parameter χ eff measured from a merger event, which is the mass-weighted average of the sBH spins projected along the BBH orbital angular momentum axis, can help constrain compact-object merger pathways (e.g., Gerosa et al. 2018;Bavera et al. 2020;Tagawa et al. 2020aTagawa et al. , 2021Wang et al. 2021b). The χ eff distribution inferred from the observed GW events prefers low values on average (The LIGO Scientific Collaboration et al. 2021), which suggests low natal sBH spins or misalignment between the binary orbital angular momentum and sBH spins (Farr et al. 2017). Tagawa et al. (2020a) inferred that in the AGN disk environment, it would require a moderate standard deviation of the initial spin |a| < 0.5 for embedded sBHs to produce the observed χ eff distribution.
In contrast, Jermyn et al. (2021) postulated that during the stellar evolution phase, disk-embedded stars typically grow to large masses and significantly spin up to critical rotation as they continuously receive injections of angular momentum from disk gas. This process may lead to the common formation scenario of sBHs with high spin aligned with the AGN disk, born from fastspinning stars (Shapiro & Shibata 2002;Dennison et al. 2019). Subsequently, sBHs' rapid rotation may be upheld as they continue to accrete gas from the AGN disk. This scenario may be inconsistent with the general low-χ eff distribution in LIGO events. Although potentially efficient angular momentum transfer mechanisms may lead to low-spin remnants during the core collapse (Qin et al. 2018;Fuller & Ma 2019), sBHs with a negligible spin parameter (a ∼ 0) can still gain angular momentum through gas accretion and quickly become aligned with each other (also with the disk rotation). Such a process would lead to high-spin merger events unless dynamical encounters between BBHs and other stars or sBHs can tilt their orbital planes significantly away from the disk plane (Tagawa et al. 2020a). For the purpose of studying gravitational-wave production through the AGN channel, it is essential to understand the details of the spin evolution of diskembedded sBHs.
In this work, we show that a finite orbital eccentricity e  0.05 of embedded massive star/sBH orbits around the supermassive black hole (SMBH) can significantly influence their spin evolution. In the corotating frame centered on a companion, the unperturbed background Keplerian motion of the global disk flows in the retrograde direction with respect to its orbital motion around the SMBH. While the classical theory of companion-disk interaction suggests the formation of a prograde disk around the companion with a circular orbit, due to steady tidal perturbation (e.g., Korycansky & Papaloizou 1996;Lubow et al. 1999), a companion with a moderate orbital eccentricity may disrupt largescale steady flow structures and restore background shear. Recent work by Bailey et al. (2021) showed that eccentric planets embedded in a protoplanetary disk are capable of reversing the rotation of circumplanetary gas flows from prograde (in the case of a circular planet) to retrograde, albeit they did not discuss the criterion and implication of this effect in details. We point out this phenomenon is also generic for low-mass companions (massive stars, sBHs) embedded in an AGN disk, as it is not uncommon for them to obtain eccentricities due to birth kicks (Lousto et al. 2012) and other dynamical interactions (Zhang et al. 2014;Secunda et al. 2019Secunda et al. , 2020. In particular, moderately eccentric orbiters endure spin-up in the retrograde direction and become misaligned (or nearly antialigned) with circular orbiters. Merging of misaligned (antialigned) pairs of sBHs may produce low-χ eff events consistent with most LIGO detections.
This Letter is organized as follows: In Section 2, we introduce the numerical setup of our global disk simulation, and in Section 3 we analyze results of circular versus eccentric-disk-embedded angular momentum accretion, focusing on how eccentricity changes the direction of the flow of the circumstellar disk (CSD) from prograde to retrograde. In Section 4 we discuss the influence of eccentricity on the spin evolution of individual massive stars and sBHs in AGN disks. In Section 5 we give a brief summary and discussion of our work.

Numerical Setup
In this section, we apply fiducial numerical simulations to study the effect of orbital eccentricity on the accretion and spin evolution of embedded massive stars/sBHs (generally referred to as companions hereafter) in AGN disks with LA-COMPASS (Li et al. 2005(Li et al. , 2009). The embedded companion orbits around the central SMBH with its semimajor axis of r 0 and eccentricity of e. Different e are explored to map the transition of CSD flow from the prograde to retrograde configuration. For our initial conditions, we choose an axisymmetric and locally isothermal 2D accretion disk model, with a constant aspect ratio h 0 ≡ H/r over the whole disk and a constant α-viscosity with α = 0.001 (Shakura & Sunyaev 1973). The gas surface density profile is initialized as a power-law profile Σ(r) ∝ r −0.5 with a total disk mass of 10 −3 M SMBH , where M SMBH is the mass of the SMBH.
We expect companions to be embedded in turbulent, possibly star-forming regions of an AGN disk with marginal gravitational stability (Goodman 2003;Sirko & Goodman 2003;Thompson et al. 2005). However, for simplicity, we ignore self-gravity throughout the simulation domain such that in a steady state, a constant accretion rate is given by the α prescription. By construction, our disk model is assumed to be laminar, albeit we comment on how realistic gravitoturbulence might affect our results in Section 5. In addition, when the feedback from the accreting companion is ignored, the magnitude of the local surface density Σ 0 at r 0 in code units S = S r M 0 0 0 2 SMBH is only a normalization factor in the gapopening and mass-accretion processes, and it does not affect scale-free results (Kanagawa et al. 2015b;Tanigawa & Tanaka 2016;Li et al. 2021a).
For typical companion and SMBH masses, the mass ratio q = M c /M SMBH  10 −5 , where M c is the mass of the embedded companion, and the scale height of the accretion disk at r 0 is h 0 ∼ 0.01 (Levin 2007;Nayakshin & Cuadra 2007). It is unfeasible to resolve circumstellar flow for global simulations with such a small CSD, so we do the rescale problem with h 0 = 0.05, q = 3 × 10 −5 , and α = 0.001 for our fiducial simulation (Li et al. 2021c), which mimics the more realistic case with h 0 = 0.01, q = (1 ∼ 5) × 10 −6 and α = 0.01 ∼ 0.1 in terms of the gap-opening parameter a = --K q h 2 0 5 1 (Kanagawa et al. 2015a(Kanagawa et al. , 2015b(Kanagawa et al. , 2018. With those parameters, the companion's Bondi radius for the fiducial h 0 is = GM c r 0.012 c s 2 0 , where c s = h 0 v K is the isothermal sound speed, and v K is the local Keplerian velocity. This radius is mildly smaller than both the average companion Hill radius R H = (q/3) 1/3 r 0 = 0.022 r 0 and the scale height H.
The companion's direct potential term is softened as where δr is the gas fluid's distance from the companion and ò = 0.02 R H is the softening parameter applied to yield numerical convergence and resolve the CSD around the companion. In order to simulate companion accretion, we follow Li et al. (2021a) (see also ;Kley 1999;Tanigawa & Watanabe 2002;D'Angelo et al. 2003;Dobbs-Dixon et al. 2007) in introducing a small sinkhole with a radius of r s , which is the same as the softening scale ò centered around the embedded companion (which is effectively three central grids connected in azimuth). The average removal timescale of surface density within r s is chosen to be the dynamical time W -K 1 at r 0 . As confirmed by previous works (Tanigawa & Watanabe 2002;Li et al. 2021a), the accretion rate converges to the same quasi-steady rate for small enough r s < 0.1R H independent of specific values of r s and removal rate.
In principle, the values of both r s and ò should be chosen to mimic the companion's physical (or "atmospheric") size R c within which the gas flow becomes irrelevant. A low-mass companion's gravity can only strongly modify the background gas flow within R B , but if ò ∼ R B , there would be no direct modification because the direct potential gradient within R B is already reduced to nearly zero. We checked several previous simulations of companion-disk interaction with ò ∼ R B , usually adopted for convenience when studying global effects of companion-disk interactions, and found no apparent circumcompanion flow pattern (Chen et al. 2021;Li et al. 2021a). In the AGN context, the ratio of R c /R B is given by where T is the accretion disk temperature around r = r 0 , typically 10 4 -10 5 K, and μ is the gas molecular weight. Only in the extreme case of red giants may R c grow to be comparable to R B , so it is important we capture flow patterns enabled by R c = R B with a small ò in our study. Our choice with a small ratio of ensures that the trend of flow patterns around the companion is not significantly impacted by the softening scale r s and ò.
A damping method is used to avoid artificial wave reflection (de Val-Borro et al. 2006) and sets both the inner and outer boundary to the initial conditions. Our simulation domain extends from 0.5 r 0 to 2 r 0 and has a uniform hyperresolution of n r × n f = 3072 × 16,396 with ∼30 grids per R B .

Circumstellar Flow Patterns of Eccentric Companions
In our fiducial models, we adopt h 0 = 0.05 and experiment with a circular orbit as well as a wide range of eccentricities e (see Figure 4 for all parameters). The snapshots for the gas surface density at pericenter and apocenter are shown in Figure 1. Weak spiral arms can be seen in the vicinity of the companion (see Zhu & Zhang 2022 for a more detailed discussion on the dependence of global spiral structure on eccentricity). Illustrative comparison of CSD flow patterns, at pericenter and apocenter, for exemplary cases e = 0.02, 0.025, and 0.05, are shown in Figure 2. The gas velocity vectors shown as white arrows are overlaid with the specific angular momentum of gas distribution in the comoving frame of the companion. All results shown in this section are collected at ∼500 orbits when the flow field has already reached semisteady state periodic in orbital phase. 5 In circular and nearly circular cases, the unperturbed initial Keplerian shear of the background rotation of the accretion disk is effectively retrograde with respect to the companion (blue everywhere) in its close vicinity. However, this kinematic distribution is different from the steady-state midplane flow pattern of gas within the horseshoe and CSD regions, which have been extensively studied in the circular cases from gas giants down to sub-Earth-mass planet embryos (e.g., Korycansky & Papaloizou 1996;Lubow et al. 1999;Machida et al. 2008;Tanigawa et al. 2012;Ormel 2013;Fung et al. 2015;Szulágyi et al. 2022) with a central host. In the left panels of Figure 2, we show the flow pattern around a companion with e = 0.02, for eccentricities below which the flow patterns converge with those of the circular case. With the U-turns produced at the azimuthal limits of the companion Bondi radius (shown as black solid circles), the Keplerian shear is deflected by the tidal torque into a prograde CSD flow to connect smoothly with the horseshoe streamlines, manifesting as a stripe of red regions with prograde angular momentum relative to the companion.
This physical mechanism was ignored in Jermyn et al. (2021). In their treatment, the companion spin-up rate is derived by assuming the distribution of gas angular momentum with respect to the companion follows only the averaged initial retrograde Keplerian shear (see their Figure 1). In reality, this so-called "shear-dominated" gas accretion regime is only valid when R c  R B or if other effects (e.g., magnetic torques, large pressure gradient) could influence the gas rotation profile up to some critical radius R T > R B , preventing the flow direction switch around R B . We suggest that most companions (especially sBHs) have R c = R B , and near-circular orbiters always maintain a prograde CSD (red region) instead of accreting directly from the background retrograde shear (blue region).
In contrast, we found that companions with a sufficiently large orbital eccentricity (e  0.025 for our h 0 = 0.05) induce a genuine retrograde flow around themselves. The middle and right panels of Figure 2 show the flow pattern around a companion with e = 0.025 and e = 0.05 at pericenter and apocenter. For e = 0.025, as high velocity fluctuations break the cumulative effect of tidal interaction and disrupt the classical horseshoe region, the CSD orbits connect with the background shear and become effectively retrograde.
There is a subtle explanation for why red regions seem to "reappear" as e reaches larger values (e.g., in the case of e = 0.05, right panel of Figure 2). These regions of prograde angular momentum with respect to the companion are not caused by horseshoe streamlines anymore but by headwinds, as a feature of "background shear flow" that we propose an eccentric orbiter is able to preserve. At around pericenter/ apocenter, the orbital velocity of a highly eccentric companion becomes considerably faster/slower than the initial Keplerian background flow at that radius; therefore, there is a small red region to the left/right of the companion initially (shown in the lower-left corners of flow patterns corresponding to e = 0.05 in Figure 2). The width of this prograde region expands as e increases. However, because the retrograde regions on the opposite side of the companion always have a larger relative velocity magnitude, the background shear flow azimuthally averaged around the companion is still effectively retrograde. After the CSD forms, its flow orientation essentially captures the averaged initial properties and preserves the retrograde spin, but far away there is still a reminiscence of the initial prograde region to the left/right sides of the companion (the red blobs in the right panels of Figure 2), which may penetrate deeper into R B for larger e. However, these regions do not affect the retrograde rotation profile close to the companion.
To quantify average rotation profiles down to grid scale, we show the azimuthally averaged specific angular momentum as a function of distance to the companion δr in Figure 3, averaged from 10 apocenter/pericenter snapshots. The angular momentum profiles for different eccentricities are labeled with different colors, while the Keplerian profile (in orbit around the companion) is shown as the black line for comparison. In the circular and mildly eccentric cases (e  0.02), a classical smooth prograde Keplerian rotation profile is formed within R B , then the profile transforms into the retrograde background shear far away from the companion. In the retrograde cases, the profiles are also smooth for moderate eccentricities (e = 0.025, 0.05), albeit always negative (retrograde). Strong fluctuations around R B starts to appear for e = 0.1 due to penetration of the headwind region, but deep within the Bondi radius, the CSD can still form a retrograde Keplerian rotation profile.
It is worth noting that our transition eccentricity ∼0.025 is smaller than the critical value of ∼0.05 for a similar q and h 0 from 3D simulations of Bailey et al. (2021). This discrepancy is most probably due to the difference in 2D and 3D geometry because time-independent flow fields are generally more prone to disturbances in 2D as perturbations cannot dissipate in the extra vertical direction (in Section 4.1, we show that our massaccretion rate is also more sensitive to eccentricity compared to 3D case; also see Li et al. 2021a).
In reality, the circumcompanion accretion disk is not directly connected to the background flow but fed by meridional flows that transport materials down to its vertical boundaries from large heights (e.g., Tanigawa et al. 2012;Szulágyi et al. 2014), and it is not possible to capture certain important physical processes with 2D simulations, such as how the radial mass flux changes from outflow to inflow with increasing vertical height. However, we are focused on the azimuthal velocity profiles of CSDs in the midplane sufficiently close to the companion, which is shown to settle to Keplerian profiles in 3D as well as 2D simulations. In particular, Figure 6 of Bailey et al. (2021) shows the midplane azimuthally averaged rotation profiles around companions with h 0 and q disk parameters comparable to our fiducial values with various orbital eccentricities; their results demonstrate that for a small enough smoothing length (left panel), the j CSD (δr) profile deep within the Bondi radius also converges to approximately Keplerian, similar to Figure 3. One major difference is the fluctuation of the j CSD profile around R B is very large for high e = 0.1 in our simulation but relatively smooth in Bailey et al. (2021), which may be relevant with complicated 3D flow structures. Nevertheless, because the transition of the Keplerian rotation profile at δr = R B is qualitatively similar, we expect our main conclusion to be unaffected by the absence of vertical flow patterns.

Angular Momentum Budget
It is easy to understand why the j and Ω CSD profiles in the prograde/retrograde cases can be nearly mirror symmetric about the y = 0 line when δr is very small. Within the CSD region where streamlines become axisymmetric with respect to the companion, we may integrate the azimuthal equation of motion for the gas around the companion to obtain where Λ is the injection rate of the angular momentum per unit mass by tidal interaction with the SMBH. We immediately see that when δr is sufficiently small so the tidal interaction term is negligible, both a prograde Ω CSD > 0 and a retrograde Ω CSD < 0 could maintain the same profile of Ω CSD ∝ δr −3/2 when the S  M c profiles around the companion are nearly constant, which we checked to apply for most cases within a region r s < δr < R B  R H . 6 In principle, CSDs can never extend to the Hill radius before it is tidally truncated even if gas thermal motion is subdominant (Martin & Lubow 2011). Additionally, we find that even considering the perturbation term, profiles of |Ω CSD | can still be nearly identical very close to the companion. When e becomes as large as 0.1, however,  M c and Σ becomes quite nonaxisymmetric and the Ω CSD profile becomes substantially nonlinear near the Bondi radius.

Criterion for Flow Pattern Transition
In this subsection, we explore the transition of flow pattern from prograde to retrograde for different h 0 and try to identify the transition eccentricities in Figure 4. For the fiducial h 0 = 0.05, our transition eccentricity is a factor of 2 smaller than in 3D simulations. But because our transition in the inner CSD direction is shown to be qualitatively similar and orbital eccentricity most directly influences midplane hydrodynamics, we expect 2D simulations to capture basic scalings relevant to the transition boundary.
To test how the transition eccentricity depends on h 0 , we run another two sets of simulations with h 0 = 0.03 and 0.07 with q = 3 × 10 −5 . We see from Figure 4 that the critical eccentricity e for the transition from prograde to retrograde gradually increases with increasing h 0 , although the dependence on h 0 is sublinear. This dependence is natural because it is easier to destabilize steady flow patterns when the local sound speed is small. To briefly investigate whether the transition is also relevant to the companion's mass, which exclusively controls the Hill radius of the companion, we have also performed two extra runs with h 0 = 0.05, e = 0.025 but with q = 1 × 10 −5 and q = 1.5 × 10 −4 , respectively. We found that both masses are able to keep a prograde CSD, in contrast to the result for q = 3 × 10 −5 . This shows that the transition boundary, if it does exist for a range of q, may have a complicated and nonmonotonic dependence on the companion mass. This dependence should be investigated in detail in future parameter surveys with robust 3D methods.

Mass-accretion Rates
Without constraints from thermal or kinetic feedback from the companion (i.e., sBH), the hydrodynamical mass-accretion rate of an embedded massive star can be estimated by the Bondi accretion rate~S  M cR H B 0 s B 2 (e.g., Rosenthal et al. 2020) when R B  R H and the Bondi radius separates the background flow with the CSD. Above the corresponding mass ratio, the companion enters the Hill accretion regime ; see Figure 1). Figure 5 shows the measured accretion rates as a function of time in our fiducial cases, measured in units of S W r 0 0 2 K . In all cases, the mass-accretion rates reach a steady state after ∼500 orbits. In the circular case, we see the accretion rate quickly converge to around half of the Bondi rate 7´S W r 1.44 10 K 4 0 0 2 , while for eccentric cases the steady state  M c seems to be excited to considerably larger values, depending on the eccentricity. This scenario may be similar to cases in Kley & Dirksen (2006), Li et al. (2021a), and Tanaka et al. (2022) when the eccentricity of the disk itself is excited by a planet on a fixed orbit, which may generate nonlinear effects that enhance the accretion rate (the surface density within the CSD Σ also deviates from the initial condition Σ 0 ). This effect may be weaker in 3D , see Section 5.1). Bailey et al. (2021) reported that the effective relative change in  M c may be ∝ e/h 0 for eccentric companions (see their Figure 9) in 3D, though they do not have active mass removal in their treatment. In the limit of much larger eccentricity corresponding to supersonic epicyclic speeds, we expect a decrease in  M c dominated by the shrinking of the accretion impact parameter, but further analysis at even larger eccentricity is expected to be carried out and presented elsewhere. A detailed study of this 2D versus 3D discrepancy is beyond the main scope of this paper, and we first adopt the Bondi accretion rate in the estimation below.

Massive Stellar Spin Evolution
Assuming small mass loss, the timescale for the star to spin up to critical rotation can be estimated to be Here J crit is the critical angular momentum, and j is the specific angular momentum per unit accreted mass at the effective accretion radius R c . For circular orbiters, because a prograde Keplerian CSD naturally forms within the Bondi radius and can at most extend down to its physical radius R c , it could maintain a prograde j that on average lies between [ ] GM R GM R , c c c B , depending on where exactly the CSD is truncated. For eccentric cases, the distance within which a Keplerian retrograde flow can be established could still be close to R B if mildly eccentric (e.g., e = 0.025 to ∼0.05) or smaller than R B for a large enough e (e.g., e  0.1).
Around stars where the inner boundary of CSD flow extends down to ∼R c , the spin-up timescale τ s is on the order of the mass-doubling timescale for Bondi accretion   This short timescale suggests that stellar companions in AGN disks can be spun up to extreme rotations rapidly, in qualitative agreement with Jermyn et al. (2021). More importantly, we argue that the spin direction is prograde for circular orbiters and retrograde for moderately eccentric orbiters, which could potentially result in a population of antialigned high-spin massive stars.

sBH Spin Evolution
For sBHs, the accretion rate measured at our sinkhole may not be physical, because the Bondi accretion in disk environments is hyper-Eddington while the accretion onto the embedded sBHs is likely constrained by the Eddington value. The mass-growth timescale is where τ Sal = M • c 2 /L E = 5 × 10 8 yr is the Salpeter (1964) timescale, L E = 1.25 × 10 38 M • /M e erg s −1 is the Eddington luminosity, and M • is specifically used to denote the companion sBH's mass. The accretion efficiency η, which also depends on BH spin, is usually on the order of a few percent for isolated BHs (Jiang et al. 2014(Jiang et al. , 2019 or even much lower due to the strong outflow/jet for the embedded sBHs in AGN disks (Pan & Yang 2021;Tagawa et al. 2022). The gas accreted by the sBH carries the specific angular momentum of the CSD gas at R ISCO so that » j GMR • ISCO (Bardeen 1970;Tagawa et al. 2020a). For sBHs born with negligible natal spin (Fuller & Ma 2019), the accretion of prograde/retrograde gas onto sBHs evolves a • toward ±1 on an asymptotic spin-evolution timescale Here R • is the size of the embedded black hole. Their normalized spin parameter could grow to |a • | ∼ 1 within τ s < τ M , i.e., before their mass grow exponentially. Here the sign of a • is measured with respect to the prograde orbital angular momentum. The collapse of critically rotating stars may also generate sBHs with rapid spins (|a • |  1) (Shapiro & Shibata 2002;Dennison et al. 2019). In that sense, a population of misaligned BHs could form directly from the population of misaligned massive stars. Higher-generation sBHs, as merger products, can also have high natal spin regardless of the individual spins of progenitors (Hofmann et al. 2016). If these rapidly spinning sBHs attain circular orbits, they would be surrounded by CSDs with prograde spins such that a • > 0 sBHs would spin up and those with a • < 0 ones would spin down. The collapse of their main-sequence progenitors and coalescence of their predecessors may also lead to recoil and induce significant eccentricities such that they would be surrounded by CSDs with retrograde spins. In this case, sBHs with a • < 0 would spin up whereas those with a • > 0 would spin down. These effects need to be considered in future population synthesis studies of the AGN channel for BBH mergers.
It is worth mentioning that the spin-down cases pertain to a more general category of configurations where outer CSDs' rotation is inclined relative to central sBHs (with the nearly antiparallel spin being an extreme case). Due to the Lense-Thirring effect, a spinning sBH exerts a torque that warps or breaks the disk (Bardeen & Petterson 1975;Papaloizou & Pringle 1983;Papaloizou & Lin 1995;Papaloizou et al. 1998;Tremaine & Davis 2014). Interior to some critical radius R w the axis normal to the disk plane will still be approximately aligned with the sBH's spin vector. Outside R w , the angular momentum vector of the disk is aligned with that of the gas accreted from the AGN disk subject to effects of sBHs' orbital eccentricity as we have shown above.
When a steady Eddington-limited gas flow from the CSD to the sBH is uninterrupted, the timescale for the sBH's spin vector to align with the outer disk may be approximated (King et al. 2008) by Because |a • | 1 and R • , R ISCO < R w , the alignment timescale is shorter than both the mass-doubling and the spin accretion timescales. This justifies the extreme cases of the parallel or antiparallel configurations discussed above. However, it is not clear whether the outer region of the disk becomes decoupled from the inner warp regions if the Lense-Thirring torque leads to a tearing configuration (Nixon et al. 2012), especially when a • and the CSD rotation is nearly antiparallel.

Summary and Discussion
In this Letter, we performed 2D hydrodynamical simulations of stellar companions embedded in an AGN disk. We found that within their Bondi radius (usually much larger than their physical radii), while circular orbiters accrete gas from a prograde CSD, moderately eccentric orbiters accrete gas from a retrograde CSD. There is a generic transition eccentricity between prograde and retrograde CSDs dependent on disk scale heights and companion mass ratios. This effect has several implications for the spin evolution of disk-embedded companions: 1. Massive stars embedded in AGN disks would quickly spin up to critical rotation within a Bondi accretion timescale, with circular ones being spun up in the prograde direction and eccentric ones being spun up in the retrograde direction. 2. Embedded sBHs may be born from the core collapse of massive stars. If they are born with low spin, the accretion of disk gas will spin them up in the same manner as massive stars. This leads to the formation of a population of misaligned sBHs spun up by prograde and retrograde CSDs. The merging of misaligned pairs of sBHs would contribute to low χ eff GW events. 3. Both core-collapse supernova and hierarchical merger events may also produce sBHs with a high natal spin, a fraction of which would have eccentricities induced by recoil velocities leading to the formation of retrograde CSDs. They would experience spin down if their initial spin is prograde or spin up if their initial spin is also retrograde, in contrast to circular orbiters with prograde CSDs. In cases where sBH spins are misaligned or antialigned with the CSD, the spin-down process may be accelerated by the Lense-Thirring effect.
To facilitate predictions for gravitational-wave production through the AGN channel, prescription for this bifurcation in angular momentum accretion should be incorporated into future population synthesis models of sBH and BBH evolutions.
Investigation of BBHs with their center of mass on eccentric orbits around the SMBH also appeals to subsequent research. In such hierarchical systems, a circumbinary disk (CBD) will form outside the two sBHs' individual CSDs (Baruteau et al. 2011;Li et al. 2021c;Li & Lai 2022). The interaction of CBDs and CSDs with BBHs as well as the background Keplerian flow will contribute to the binaries' orbital angular momentum, whose direction and magnitude influences not only χ eff but also the inspiral timescale and spin of the merger product.
In all of our simulations, we have adopted an isothermal equation of state, which could be unrealistic within the CSD. Higher temperatures around the companion vicinity due to the contraction of gas and/or feedback from the companion may prevent the formation of CSDs (Szulágyi et al. 2016;Szulágyi 2017;Li et al. 2021b) and thus significantly affect the angular momentum accretion. Future work should implement more detailed treatments of radiative and kinetic feedback around massive stars or sBHs.
We have performed our simulations in 2D because our uniform hyperresolution treatment is quite intense with each run taking about 0.1 million CPU hours. We found 2D rotation profiles to be similar to 3D midplane results for marginally gapopening companions in a disk (Bailey et al. 2021; see their Figure 6), albeit we cannot capture certain important physical processes. Subsequent global simulations should expand our limited number of samples for eccentric orbiters and scan a larger and more realistic parameter space across (q, e, h 0 ) with less expensive adaptive mesh numerical schemes.
The background flow in our accretion disk is effectively laminar with a prescribed α. However, in a realistic gravitoturbulence state, α is driven by fluctuations in the gravitational potential and the velocity field of the turbulent background (Gammie 2001). The dominant turbulent eddies has scales of H, which is larger than R B for low-mass companions. If so, CSD rotation may be driven more strongly by turbulent eddies with randomly distributed angular momenta than by either the steady-state flow or the background shear, albeit these patterns might be recovered in a time-averaged sense (Baruteau & Lin 2010, see their Figure 10). Under this condition, the spin evolution of companions might be episodic and dominated by stochastic components. We seek to understand this different flow geometry and propose prescriptions to quantify such processes in a parallel project.
Although our main discussion revolves around sBHs in AGN disks, the generic change of flow direction around diskembedded companions is also relevant in the protoplanetary disk context, at least for planets with masses in the super-Earth regime and physical/atmosphere radii considerably smaller than their Bondi radii. Retrograde circumplanetary flows induced by supersonic orbital eccentricity can strongly influence the evolution of planetary spins and planetary sizes  and may also lead to formation of retrograde satellites. The observational implications of such impacts appeal to further exploration.