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

Spin evolution of stellar-mass Black Holes (sBHs) embedded in AGN accretion disks is an important process relevant to production of gravitaional waves from binary Black Hole (BBH) merger events through the AGN channel. Since embedded sBHs are surrounded by circum-stellar disks (CSDs), the rotation of CSD gas flows determine 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 of CSD orientation occurs at a transition eccentricity that scales nearly proportional with local sound speed. This bifurcation in the CSD flow and thereafter spin-up direction of SBHs leads to formation of a population of nearly anti-aligned sBHs and should be incorporated in future population models of sBH and BBH evolutions.


INTRODUCTION
Active galactic nucleus (AGN) disks have emerged as rich factories for producing massive stars and their remnant stellar-mass black holes (sBH), neutron stars (NSs), including some of the detected binary black holes (BBHs), neutron stars mergers, accretion-induced collapse of NSs (e.g., Artymowicz et al. 1993;McKernan et al. 2012; Bartos et al. 2017;Stone et al. 2017;Leigh et al. 2018;Gröbner et al. 2020;Davies & Lin 2020;Tagawa et al. 2020b;Kimura et al. 2021;Perna et al. 2021;Wang et al. 2021a;Zhu et al. 2021;Li et al. 2021b,c). Due to the very dense environment surrounding the embedded objects in AGN disk, this "AGN channel" is very promising for producing the recent detection of the heaviest BBH merger event GW190521 with a possible electromagnetic counterpart (Abbott et al. 2020;Graham et al. 2020), which can possibly differentiate it from other merger channels (e.g., Rodriguez et al. 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;Wang et al. 2021b;Tagawa et al. 2020aTagawa et al. , 2021a. 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 AGN disk environment, it would require a moderate standard deviation of initial spin |a| < 0.5 for embedded sBHs to produce the observed χ eff distribution.
In contrast, Jermyn et al. (2021) postulated that during 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 common formation of sBHs with high spin aligned with the AGN disk, born from fast spinning 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-Virgo 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 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 in the spin evolution of disk-embedded 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 co-rotating 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 classical theory of companiondisk interaction suggests the formation of a prograde disk around the companion with a circular orbit, due to the steady tidal perturbation (e.g. Korycansky & Papaloizou 1996;Lubow et al. 1999), a companion with a moderate orbital eccentricity may disrupt large-scale 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 circum-planetary 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, since it's 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 anti-aligned) with circular orbiters. Merging of misaligned (anti-aligned) pairs of sBHs may produce low-χ eff events consistent with most LIGO-Virgo detections. This Letter is organized as follows: In §2, we introduce the numerical setup of our global disk simulation and in §3 we analyze results of circular v.s. eccentric diskembedded angular momentum accretion, focusing on how eccentricity changes the circumstellar disk (CSD) flow's direction from prograde to retrograde. In §4 we discuss the influence of eccentricity on the spin evolution of individual massive stars and sBHs in AGN disks. In §5 we give a brief summary and discussion about 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 semi-major axis of r 0 and eccentricity of e. Different e are explored to map the transition of CSD flow from the prograde to regtrograde 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 SMBH.
We expect companions to be embedded in turbulent, possibly star-forming regions of an AGN disk with marginal gravitational stability (Sirko & Goodman 2003;Goodman 2003;Thompson et al. 2005). However, for simplicity we neglect 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 gravito-turbulence 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 Σ 0 = Σ 0 r 2 0 /M SMBH is only a normalization factor in the gap-opening 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 (Nayakshin & Cuadra 2007;Levin 2007). It's 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 sim-ulation (Li et al. 2021c), which mimics a more realistic case h 0 = 0.01, q = (1 ∼ 5) × 10 −6 and α = 0.01 ∼ 0.1 in terms of the gap opening parameter K = q 2 h −5 0 α −1 (Kanagawa et al. 2015a(Kanagawa et al. ,b, 2018. With those parameters, the companion's Bondi radius for the fiducial h 0 is GM c /c 2 s = 0.012 r 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 averaged 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) to introduce a small sinkhole with a radius of r s being same as the softening scale centered around the embedded companion (which is effectively 3 central grids connected in azimuth). The average removal timescale of surface density within r s is chosen to be the dynamical time (5Ω K ) −1 , where Ω K is the Keplerian frequency at the companion's location. 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 Keplerian circum-companion flow pattern Chen et al. 2021). 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's important we capture flow patterns enabled by R c R B with a small in our study. Our choice with a small ratio of r s /R B = /R B = 0.036(h 0 /0.05) 2 (q/3 × 10 −5 ) −2/3 1 ensures that the trend of flow patterns around the companion are 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 hyper-resolution of n r × n φ = 3072 × 16396 with ∼ 30 grids per R B .

Circum-stellar 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 eccentricity 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 co-moving frame of the companion . All results shown in this Section are collected at ∼ 500 orbits when the flow field has already reached a semi-steady state periodic in orbital phase. 1 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;Ormel 2013;Fung et al. 2015;Tanigawa et al. 2012;Szulágyi et al. 2022) with a center 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 that 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 as to connect smoothly with the horseshoe streamlines, manifesting as a stipe of red regions with prograde angular momentum relative to the companion.
This physical mechanism was neglected in Jermyn et al. (2021). In their treatment, the companion spinup 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 "sheardominated" gas accretion regime is only valid when R c R B or if other effects (e.g. magnetic torques, large pressure gradient) could influences 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 shows 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 accumulative effect of tidal interaction and disrupt the classical horseshow region, the CSD orbits connect with the background shear and becomes effectively retrograde.
There is a subtle explanation for why red regions seem to "reappear" as e reach 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 corre- sponding to e = 0.05 in Figure 2). The width of this prograde region expands as e increases. However, since 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 side 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, since timeindependent flow fields are generally more prone to disturbance in 2D as perturbations cannot dissipate in the extra vertical direction (in §4.1 we show that our mass accretion rate is also more sensitive to eccentricity compared to 3D case, also see Li et al. (2021a)).
In reality, the circum-companion 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's 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 azimuthallyaveraged rotation profiles around companions with h 0 , q disk parameters comparable to our fiducial values with various orbital eccentricity, their results demonstrate that for 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 that the fluctuation of 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, since the transition of 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 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 a same profile of Ω CSD ∝ δr −3/2 whenṀ c /Σ profiles around the companion are nearly constants, which we checked to apply for most cases within a region r s < δr < R B R H 2 . In principle, CSDs can never extend to the Hill radius before it is tidally truncated even if gas thermal motion is sub-dominant (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,Ṁ c and Σ becomes quite non-axisymmetric and the Ω CSD profile becomes substantially non-linear 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 since our transition of 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, 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 sub-linear. This dependence is natural since it's easier to destabilize steady flow patterns when 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 simulations with h 0 = 0.05, e = 0.025 but for 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 non-monotonic 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Ṁ B ∼ Σ 0 c s R 2 B /H (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 Σ 0 r 2 0 Ω 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 3 1.44 × 10 −4 Σ 0 r 2 0 Ω K , while for eccentric cases the steady stateṀ 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); Tanaka et al. (2021) 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 CSD Σ also deviates from the initial condition Σ 0 ). This effect may be weaker in 3D , see §5.1). Bailey et al. (2021) reported that the effective relative change inṀ 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Ṁ c dominated by the shrinking of 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 v.s. 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, since 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 c R c , √ GM c R B ], depending on where exactly the CSD is truncated. For eccentric cases, the distance within which a Keplerian retrograde flow can be established could be still close to R B if mildly eccentric (e.g., e = 0.025 ∼ 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 mass doubling timescale for Bondi accretion (5) This short timescale suggests that stellar companions in AGN disks can be spinned 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 anti-aligned 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 environment 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 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 black hole spin, is usually on the order of a few percent for isolated black holes (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. 2021b). The gas accreted by the sBH carries the specific angular momentum of the CSD gas at R ISCO so that j ≈ √ GM • R 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 • towards ±1 on an asymptotic spinevolution 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 are inclined relative to central sBHs (with nearly anti-parallel spin being an extreme case). Due to the Lense-Thirring effect, a spinning sBH exerts a torque which 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 time scale for sBH's spin vector to align with the outer disk may be approximated (King et al. 2008) by Since |a • | ≤ 1 and R • , R ISCO < R w , the alignment timescale is shorter than both the mass doubling and the spin accretion timescale. This justifies the extreme cases of the parallel or anti-parallel configurations discussed above. However, it is not clear whether the outer region of the disk becomes decoupled to 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 anti-parallel.

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: • 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.
• Embedded sBHs may born from core-collapse of massive stars. If they are born with low spin, accretion of disk gas will spin them up in the same manner as massive stars. This leads to formation of a population of misaligned sBHs spun up by prograde and retrograde CSDs. Merging of misaligned pairs of sBHs would contribute to subsequent low χ eff GW events.
• Both core-collapse supernova and hierarchical merger events may also produce sBHs with high natal spin, a fraction of which would have eccentricities induced by recoil velocities leading to 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 CBD 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 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 effect 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 since our uniform hyper-resolution 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 gap-opening 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 realistic gravito-turbulence 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, the 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 (Dubois et al. 2014). 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 disk-embedded 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 circum-planetary flows induced by super-sonic orbital eccentricity can strongly influence the evolution of planetary spins, planetary sizes , and may also lead to formation of retrograde satellites. The observational implications of such impacts appeal to further exploration.