Simulating the tidal disruption of stars by stellar-mass black holes using moving-mesh hydrodynamics

In the centers of dense star clusters, close encounters between stars and compact objects are likely to occur. We study tidal disruption events of main-sequence (MS) stars by stellar-mass black holes (termed $\mu$TDEs), which can shed light on the processes occurring in these clusters, including being an avenue in the mass growth of stellar-mass BHs. Using the moving-mesh hydrodynamics code \texttt{AREPO}, we perform a suite of hydrodynamics simulations of partial $\mu$TDEs of realistic, \texttt{MESA}-generated MS stars by varying the initial mass of the star ($0.5\,{\rm M}_{\rm \odot}$ and $1\,{\rm M}_{\rm \odot}$), the age of the star (zero-age, middle-age and terminal-age), the mass of the black hole ($10\,{\rm M}_{\rm \odot}$ and $40\,{\rm M}_{\rm \odot}$) and the impact parameter (yielding almost no mass loss to full disruption). We then examine the dependence of the masses, spins, and orbital parameters of the partially disrupted remnant on the initial encounter parameters. We find that the mass lost from a star decreases exponentially with increasing distance of approach and that a $1\,{\rm M}_{\rm \odot}$ star loses lesser mass than a $0.5\,{\rm M}_{\rm \odot}$. Moreover, a more evolved star is less susceptible to mass loss. Tidal torques at the closest approach spin up the remnant by factors of $10^2$--$10^4$ depending on the impact parameter. The remnant star can be bound (eccentric) or unbound (hyperbolic) to the black hole: hyperbolic orbits occur when the star's central density concentration is relatively low and the black hole-star mass ratio is high, which is the case for the disruption of a $0.5\,{\rm M}_{\rm \odot}$ star. Finally, we provide best-fit analytical formulae for a range of parameters that can be incorporated into cluster codes to model star-black hole interaction more accurately.

In this work, we examine TDEs of stars by stellar-mass black holes or SBHs (termed µTDEs).These are less studied and have only garnered interest recently.Although µTDEs have lower observable rates (Perets et al. 2016), they can shed light on the processes occurring in the centers of globular and nuclear star clusters where they are most likely to occur.In particular, they are an important avenue in the mass growth of SBHs to form intermediate-mass black holes (IMBHs, e.g., Stone et al. 2017;Rizzuto et al. 2023).However, these studies make simplistic assumptions about the mass accreted onto a black hole after a TDE, which highlights the need for detailed simulations to model this interaction more accurately for the next generation of globular ⋆ corresponding author cluster simulations.Partial tidal disruptions (PTDEs), which are more likely than full tidal disruptions (FTDEs), can be responsible for the tidal capture, and subsequent encounters, of the remnant star by the BH if the remnant ends up in a bound orbit.The remnant stars of such interactions tend to be spun up from the torque due to the BH and can have peculiar internal structures (e.g., Alexander & Kumar 2001).
µTDEs are also of interest in the context of low-mass Xray binaries (LMXBs) and the newly-discovered Gaia BHs (El-Badry et al. 2023a,b;Chakrabarti et al. 2023).Isolated binary evolution cannot satisfactorily explain the observed rates of binaries with BHs and low-mass stellar companions (e.g., Podsiadlowski et al. 2003;Kiel & Hurley 2006).Hence, the dynamical formation of such binaries in clusters, through PTDEs and subsequent tidal captures, may very well be a crucial process to explain them.
Accurate modeling of post-disruption orbits and the internal structure of remnant stars requires detailed hydrodynamics simulations.The first detailed study on µTDEs was carried out by Perets et al. (2016).They estimated that these events could occur in Milky Way globular clusters at a rate of 10 −6 yr −1 MWGal −1 , and might observationally resemble ultra-long GRB events.Similar rates were obtained for scenarios where a supernova natal kick to a newly formed BH results in a chance encounter with its binary companion.Wang et al. (2021) used smoothed-particle hydrodynamics (SPH) to find the fallback rates onto the BH in the case of partial µTDEs of polytropic stars.Kremer et al. (2022) performed a large suite of SPH simulations of partialand full-disruptions, with varying BH and stellar masses, stellar polytropic indices, and impact parameters.Among other things, they observed that stellar structure plays a crucial role in the boundedness of the remnant stars after partial disruption, with the possibility of less centrally 'concentrated' stars ending up in hyperbolic orbits.Xin et al. (2024) looked at low-eccentricity 'tidal-peeling' events of realistic MS stars, with the star being slowly stripped away in multiple orbits.
Hydrodynamics studies on three-body encounters involving SBHs in globular clusters have also garnered interest, motivated by the significantly larger cross-sections of binary-single tidal interactions.Lopez et al. (2019) studied the interaction of binary BHs and polytropic stars and showed that tidal disruption can alter the spin of one of the BHs.An extensive series of studies was carried out by Ryu et al. (2022Ryu et al. ( , 2023aRyu et al. ( ,b, 2024b)), using both SPH and moving-mesh codes, on the different combinations of close encounters between realistic MS stars and SBHs.
In our study, we use the moving-mesh code AREPO (Springel 2010;Pakmor et al. 2016;Weinberger et al. 2020) to simulate µTDEs of low-mass MS stars with SBHs.It should be noted that Ryu et al. (2023aRyu et al. ( ,b, 2024b) also employed AREPO to model TDEs with remarkable success.We generate realistic stellar models from detailed 1D MESA (Paxton et al. 2011) profiles.TDEs of MESA-generated stars by SMBHs and IMBHs have been investigated in the past (e.g., Gallegos-Garcia et al. 2018;Golightly et al. 2019;Law-Smith et al. 2019, 2020;Ryu et al. 2020a,b,c,d;Kıroglu et al. 2023), but µTDE studies have typically probed polytropic stars.In this paper, we focus on the dependence of post-disruption mass, spin, and orbital parameters on the initial conditions and provide best-fit functions for the same.These fits can be incorporated into N-body or other cluster codes for better treatment of star-BH encounters.
The paper is organized as follows.In Section 2, we detail our simulation methodology, including our grid of initial conditions.We briefly describe the analysis of certain quantities in Section 3 and present our results and best-fit functions in Section 4. Section 5 and Section 6 discuss the implications of our work and conclude, respectively.
To solve the fluid equations, AREPO builds an unstructured Voronoi mesh with cells of varying volume and calculates fluxes between the cells in a finite-volume approach.The Voronoi mesh moves in time according to the approximate bulk velocities of the fluid elements.Gravity is handled in a tree-particle-mesh scheme (see Bagla 2002), with the minimum gravitational softening for the gas (star) cells set as one-tenth of the smallest gas cell for accuracy.AREPO allows for adaptation of spatial resolution according to arbitrary criteria on top of the default adaptivity to density, inherited from the near-Lagrangian nature of the scheme.It also includes particles that interact only gravitationally.We assume Newtonian gravity and do not include magnetic fields in our simulations.

Stellar models
We use the 1D stellar modeling code MESA (Paxton et al. 2011) to generate main-sequence (MS) stellar models with initial masses m ⋆ = 0.5 M ⊙ and 1.0 M ⊙ , and metallicity Z ⋆ = 0.02 (near-solar) at different ages.These masses are typical for globular clusters that predominantly consist of old, low-mass stars (e.g., Salaris & Weiss 2002;De Angeli et al. 2005;Marín-Franch et al. 2009), while the metallicity is higher than those of most clusters, which typically have a tenth (or less) of Solar metallicity (e.g., Harris 2010)).However, the effect of metallicity is not expected to be significant and is quantified in Section 4. MESA solves the stellar structure equations to compute the evolution and provides us with stellar mass-and radial-profile parameters (e.g., densities, temperatures, and chemical abundances) at different evolutionary phases.We convert these 1D profiles to 3D AREPO initial conditions using the scheme adopted by Ohlmann et al. (2017), with the Helmholtz equation of state (Timmes & Swesty 2000).
After generating the star, we relax it for five stellar dynamical timescales t dyn,⋆ = (r 3 ⋆ /Gm ⋆ ) 0.5 , where r ⋆ is the stellar radius and G is the gravitational constant.We then use this relaxed star for the disruption simulations.
Stellar density profiles are a major factor in determining whether a star undergoes a partial tidal disruption (PTDE) or a full tidal disruption (FTDE) for a given periapsis distance (e.g., Ryu et al. 2020c).To that end, we consider MESA models of three evolutionary stages of 1.0 M ⊙ MS stars with varying core hydrogen abundances H c , which correspond to distinct density profiles (see also Gallegos-Garcia et al. 2018;Golightly et al. 2019;Law-Smith et al. 2019;Goicovic et al. 2019).The first is close to the onset of hydrogen burning (H c ≈ 0.70), i.e., zero-age main-sequence (hereafter ZAMS), the second is approximately midway through the main-sequence (H c ≈ 0.34), i.e., middleage main-sequence (hereafter MAMS), and the third is close to the depletion of core hydrogen (H c ≈ 0.00), i.e., terminal-age main-sequence (hereafter TAMS).The total stellar mass remains nearly constant throughout the MS lifetime owing to the insignificant wind mass loss.In the case of the 0.5 M ⊙ stars, we exam-  ine only ZAMS profiles (H c ≈ 0.70).This choice is motivated by the fact that their internal structure barely changes over time (see Figure 2) owing to their large MS lifetimes.One could choose a 0.5 M ⊙ star at an age comparable to a Hubble time (e.g., Ryu et al. 2020a,c) but the results are expected to be very similar.Table 1 lists the stellar parameters of the MS star models for the stellar masses and ages we simulate.
Furthermore, we performed resolution scaling tests to ensure that the PTDE simulations were robust.To this end, we varied the number of fluid (star) cells making up a 1 M ⊙ star -5 × 10 4 , 8 × 10 4 , 2 × 10 5 , 5 × 10 5 , and 8 × 10 5 .We then performed PTDE simulations involving a 1 M ⊙ star and a 10 M ⊙ BH for each resolution.The masses lost due to tidal disruption plateaued close to the resolutions 2 × 10 5 and 5 × 10 5 , with the mass loss val-ues using 5 × 10 5 and 8 × 10 5 cells being essentially the same.Hence, we chose an initial resolution of 5 × 10 5 fluid cells for the stars in this study.It should be noted that the number of cells increases as the simulations progress since AREPO automatically adds extra cells when the density gradient is high.
Figure 1 shows that more evolved 1 M ⊙ MS stars have a denser (left panel) and hotter (right panel) core and a puffier outer layer than less evolved MS stars.Since the three stars have different radii r ⋆ , they also have different t dyn,⋆ and different tidal radii r t = (m BH /m ⋆ ) 1/3 r ⋆ for the same BH mass m BH .Figure 2 shows that 0.5 M ⊙ stars of different ages -ZAMS (H c = 0.70) and ∼ 13.5 Gyr (H c = 0.60) -have very similar density and temperature profiles, which justifies excluding the latter from our simulation runs.It should be noted that the AREPO profiles, espe- Notes.The columns indicate the initial parameters -the stellar mass m ⋆ , the BH mass m BH , the core hydrogen fraction H c , and the impact parameter b.In each row, every comma-separated combination of values in all columns is simulated, to give a total of 58 simulations.
cially the temperature, diverge significantly from the MESA ones near the surface of the star (corresponding to a fractional mass of ∼ 10 −4 -10 −3 ).This amount of mass near the surface that deviates from the MESA model sets the minimum amount of stripped mass in PTDEs that we can resolve.Therefore, we only present the results of PTDEs yielding a mass loss ≳ 10 −4 m ⋆ in this paper.

Black hole
The black holes are initialized in AREPO as non-rotating sink particles that interact gravitationally and grow in mass through accretion.We set the gravitational softening length of the BH to be ten times the minimum softening length of the gas (star) cells.The scheme we use for accretion onto the BH, described here briefly, is the same as the one used in Ryu et al. (2022Ryu et al. ( , 2023aRyu et al. ( ,b, 2024b)).Firstly, the cells around the BH within a radius of 10 4 r g , where gravitational radius r g = Gm BH /c 2 , are identified.Secondly, the accretion onto the BH is calculated using a weighted average of radial flux to account for the higher contribution of gas closer to the BH and vice-versa.Thirdly, mass is extracted from the neighboring cells of the BH, and momentum is inserted in the BH to conserve these quantities.Finally, the BH is spun up from accretion by adding angular momentum from the selected neighboring cells.The radiation feedback from accretion is not taken into account.Although angular momentum is conserved through this scheme, there is no guarantee that the total energy is conserved.However, given that accretion onto the black hole is small in our simulations, the energy lost due to accretion is negligible.
We run simulations with SBH masses of 10 M ⊙ and 40 M ⊙ , which are in proximity to the lower-and higher-mass peaks of observed LIGO BH masses (e.g., Abbott et al. 2021Abbott et al. , 2023)).Given the mass of the BH m BH , the mass of the star m ⋆ , and the radius of the star r ⋆ , the classical tidal radius is defined as

Initial conditions
The distance of the closest approach is characterized by the impact parameter: b ≡ r p /r t 1 , where r p is the periapsis distance of the initially parabolic orbit.We vary the mass of the star: 0.5 M ⊙ and 1 M ⊙ , and the mass of the black hole: 10 M ⊙ and 40 M ⊙ .We then generate a suite of simulations of b varying from 0.25 to 2.00 (1 M ⊙ stars) or 2.50 (0.5 M ⊙ stars), and the three stellar ages (density profiles) to determine the post-disruption orbital, mass and spin parameter.The additional simulations for 0.5 M ⊙ and b = 2.50 are performed to properly analyze the PTDE trends since 0.5 M ⊙ stars are fully disrupted for b ≲ 1.00.
1 β ≡ b −1 is often introduced in the literature for TDEs by SMBHs as the 'penetration factor'.
We also assume that the initial stars are non-rotating and that the star-BH encounters are parabolic.This can be justified for encounters in globular clusters because the velocity dispersions of Milky Way globular clusters are in the range ∼ 1-20 km s −1 (e.g., Baumgardt & Hilker 2018).For instance, given a velocity dispersion σ ≃ 10 km s −1 , |1 − e| ≃ 10 −4 (very close to parabolic) for the parameters considered in this paper.
Each simulation is run for ∼ 70 t dyn,⋆ after periapsis passage.This duration was chosen to ensure that the stellar masses reach steady values for proper analysis.Table 2 provides an overview of the parameters we use in our suite of simulations.

Calculation of bound and unbound mass
To find the center of mass of a self-bound object (i.e., star before disruption, remnant after disruption) in each snapshot, we used an iterative procedure fairly similar to Guillochon & Ramirez-Ruiz (2013), with some differences to improve the accuracy of the identification of bound and unbound cells.For each snapshot, we chose the gas cell with the maximum density as the 'initial' guess for the center of mass (CoM) position.Subsequently, we calculate the specific total energies of all the star cells relative to the star CoM, ε cell,CoM , and the BH, ε cell,BH , as the sum of their relative kinetic, potential and internal energies.We consider a cell to be bound to a self-gravitating object if ε cell,CoM < 0, ε cell,CoM < ε cell,BH .The last condition ensures that the cells bound to the star are 'close' to the CoM.We compute a 'new' CoM position (and velocity) using these bound cells and determine the 'new' specific energy.We iterate this process until there is no change in the relative position of the CoM.Subsequently, we calculate its 'final' bound mass.In the case of an FTDE, we also visually inspect the surface density plots and assign the mass bound to the star to be zero.
The mass bound to the BH is calculated in a similar fashion, though the position of the BH is trivially known from the simulation.Finally, we consider any cells with ε cell,CoM ≥ 0 and ε cell,CoM ≥ 0 unbound from the system.

Calculation of orbital and spin parameters
For each snapshot post-disruption, we compute the instantaneous Keplerian orbital parameters -the semimajor axis a and orbital eccentricity e -using the current values of the mass bound to the star m ⋆ , the mass of the BH (plus the gas mass bound to it) m BH , the positions and velocities of the star's CoM (r CoM , v CoM ) and the BH (r BH , v BH ).For completeness, with r ≡ r CoM − r BH and v ≡ v CoM − v BH , the equations are as follows: Article number, page 4 of 17 Notes.The columns indicate the polytropic index γ, the impact parameter for full disruption b FTDE , the inverse of the density concentration factor ρ −1 conc and the aforementioned parameter bFTDE .The b FTDE values for polytropic stars are from Mainetti et al. (2017).
We determine the spin angular momentum L ⋆ about the star's CoM using the masses m i , relative positions r i and velocities v i of the bound gas cells.In the cases of full disruption, we ignore the orbital and spin parameters.
The 'final' post-disruption mass, orbital, and spin parameters reported in the following section are the means and standard deviations of these quantities during the last ten snapshots of each simulation (the time between consecutive snapshots was chosen to be approximately equal to the initial dynamical timescale of the star).In the cases where a partially disrupted star returns on a second passage within the simulation time, we chose snapshots close to the apoapsis of the first passage for the parameter calculations.

Results
Figure 3 shows selected snapshots of our suite of simulations.It can be visually seen that the final outcomes (FTDE or PTDE) depend on all the parameters that we considered -the stellar mass m ⋆ , the core hydrogen fraction H c , the BH mass m BH , and the impact parameter b.
The following sections detail the quantitative results.It should be noted that, henceforth, initial star parameters are given the subscript ⋆ (e.g.m ⋆ , L ⋆ ), while the post-disruption remnant parameters are given the subscript ⋄ (e.g.m ⋄ , L ⋄ ).Table A.1 provides detailed initial and post-disruption parameters for our 58 simulations.Ryu et al. (2020b) showed that, in the case of SMBHs, the critical impact parameter for full disruption, b FTDE , is inversely proportional to the density concentration factor ρ conc ≡ (ρ c /ρ) 1/3 , where ρ and ρ c are average and central stellar density respectively (see also Law-Smith et al. 2020).A larger value of ρ conc represents a denser core with a puffier outer layer, and viceversa.This indicates that a star with higher core density requires a smaller impact parameter to fully disrupt.Quantitatively, b FTDE can be written as:

Effect of density concentration factor
Here, bFTDE is the critical impact parameter for the full disruption of a hypothetical star with uniform density, i.e., ρ conc = 1 (see also Equation 15of Ryu et al. 2020b).Table 3 displays the values of b FTDE of stars with polytropic indices γ (where P ∝ ρ γ ) 4/3 and 5/3 as calculated by Mainetti et al. (2017).Shown also are the analytically computed ρ conc for these polytropes, and the subsequently evaluated (using Equation 3) bFTDE values.We also estimate bFTDE to be 1.95 ± 0.20 from our suite of simulations by fitting the remnant masses for given initial stellar and BH masses (see the following section for details), and find the values of b FTDE .Despite the significant difference in the masses of the black holes, our estimated values agree with those from Mainetti et al. (2017).Estimates of this factor calculated by Ryu et al. (2020b), Law-Smith et al. 2020 andKıroglu et al. (2023) all lie within the error range.Since the former two studies are on SMBHs and the latter is on IMBHs, we reiterate that bFTDE is almost independent of the mass of the BH for a very large range of BH masses.The consensus between the previous works in spite of different numerical schemes (general relativistic gridbased hydrodynamics, adaptive-mesh refinement, and SPH respectively) and the current work strengthens this point.Therefore, we will use our computed value of bFTDE = 1.95 ± 0.20 for the fit functions in the following sections.
For completeness, Figure 4 (and Table 1) shows the trend ρ −1 conc for MESA 1 M ⊙ MS stars, with our three stellar models highlighted.We see that a more evolved star has a smaller ρ −1 conc (harder to fully disrupt) compared to a less evolved star (easier to fully disrupt) and that this relation is close to linear as a function of stellar age t ⋆ .Figure 5 shows the averaged ρ −1 conc values for a range of masses of MESA-generated MS stars (0.1 M ⊙ -20.0 M ⊙ ) of two metallicities -0.02 (near-Solar) and 0.002 (sub-Solar).We see that higher mass MS stars have denser cores and puffier outer layers (hence lower ρ −1 conc values, which plateau for stars of masses ≳ 1 M ⊙ at ∼ 1.6-2.0),while the opposite is true for lower mass MS stars.Moreover, metallicity affects ρ −1 conc only slightly, with the metal-poor stars being slightly less (more) centrally concentrated when m ⋆ > 1.0 M ⊙ (m ⋆ ≤ 1.0 M ⊙ ).The (approximate) values of ρ −1 conc are crucial to applying our fits (see the following sections).

Post-disruption stellar and BH masses
After a disruption event, a star loses mass, some of which gets bound to (and eventually accreted onto) the BH while the rest of it becomes unbound from the system.Evidently, the closer the approach of the star is to the BH, the larger the mass lost from the star.
Figure 6 shows this trend for the case of a 1 M ⊙ (left) and a 0.5 M ⊙ (right) star that is disrupted by a 10 M ⊙ (top) and a 40 M ⊙ (bottom) BH.The post-disruption fractional mass loss ∆m ⋄ /m ⋆ , where ∆m ⋆ is the difference between the star's initial mass m ⋆ and the remnant mass m ⋄ .It decreases almost exponentially with the impact parameter b.Moreover, a more evolved MS star loses less mass compared to a less evolved star for the same value of b.This is because a more evolved star, with a higher central concentration, requires a stronger tidal force (or smaller b) to strip the same amount of mass when compared to a less evolved star.Consequently, b FTDE , corresponding to ∆m ⋄ /m ⋆ = 1, is also lower for more evolved stars.
Furthermore, we see that roughly half the mass lost from the star stays bound to the BH, although this is not always true, especially when there is little mass loss (see Table A .1).This bound mass decreases through time due to continuous debris interaction resulting in mass unbinding, and hence is an upper limit for the mass that can be accreted onto the BH.Details of the accretion process and feedback from the BH might prevent some of the mass bound to the BH from being accreted, but following this process is beyond the scope of our paper.The other half of the stripped mass is unbound from the BH-star system.We can fit (with standard 1σ errors) the fractional mass loss ∆m ⋆ well by:

Post-disruption stellar spins
An initially non-rotating star gains spin after a tidal encounter due to tidal torques from the BH.Although we do not delve into the details of the spin-up, this results in differential, and not rigid, rotation of the star.It should be noted that, due to spurious motions in our initial AREPO stellar models, the spin angular momenta L ⋆ are not zero but ∼ 1-2 × 10 47 g cm2 s −1 (∼ 2 × 10 46 g cm 2 s −1 ) for 1 M ⊙ (0.5 M ⊙ ) stars.Since these values correspond to negligible equatorial surface velocities 2 of ∼ 10 m s −1 , the stars are considered to be non-rotating.The spin angular momentum of the remnant L ⋄ depends on the internal structure of the original star and the impact parameter.
Figure 7 shows the trend of L ⋄ with respect to the impact parameter b for the case of a 1 M ⊙ (left) and a 0.5 M ⊙ (right) star being partially disrupted by a 10 M ⊙ (top) and a 40 M ⊙ (bottom) BH.In addition, we indicate the approximate, orderof-magnitude estimates of the break-up angular momenta (assuming uniform density and rigid rotation) of the remnant stars, 3 .Here, we crudely estimate the remnant radius as (r ⋄ / R ⊙ ) = (m ⋄ / M ⊙ ) 0.7 and assume the breakup velocity to be the Keplerian velocity at r ⋆ , i.e., v ⋄,break = √ Gm ⋄ /r ⋄ .In reality, the remnant is significantly puffed up, and the break-up velocity would depend on the true moment of inertia of the remnant (and thus, its density profile) and the differential rotation due to spin-up.For our range of simulation parameters, L ⋄ remains below L ⋄,break , although significant mass loss can bring the remnant very close to break-up, e.g., m ⋆ = 0.5 M ⊙ and b ≲ 1.00.For higher values of b (but still less than b ∼ 2), L ⋄ decreases almost exponentially as b increases.On the other hand, when b is less such that a fractional mass loss is ≥ 2 × 10 −1 , the trend reverses until FTDE.This trend reflects the counterbalance between mass loss and spin-up: for stronger PTDEs, although the remnant is spun up more rapidly, the remnant mass is smaller.
The large spin-up for smaller b indicates that the specific angular momentum monotonically increases as b decreases.Motivated by this, we considered the scaled (moment of inertiaadjusted) parameter L ⋄,scaled ≡ L ⋄ /(m ⋄,frac r2 ⋄,frac ).Here, the fractional remnant mass m ⋄,frac ≡ m ⋄ /m ⋆ and the post-disruption estimated fractional radius r⋄,frac / R ⊙ = m 0.7 ⋄,frac is a crude estimate of the radius of an MS star.
We found the following exponential fit (with standard 1σ errors) describes L ⋄,scaled well: From the definition of r⋆,frac before, we have: Here, (∆m ⋄ /m ⋆ ) is estimated using Equation 4. This equation implies that the post-disruption spin depends significantly on the fractional mass loss.Moreover, there is no dependence on m BH .Figure 7 plots these best-fit curves, with the error bars shaded.

Post-disruption orbital parameters
Depending on the stellar structure and the BH mass, a star in an initially parabolic orbit ends up in an eccentric or even a hyperbolic orbit after a PTDE (see also Ryu et al. 2020c;Kremer et al. 2022;Kıroglu et al. 2023).The specific orbital energy of the star-BH system, ε orb , is the sum of the specific energy injected into stellar tides, ε tide , the specific binding energy of the ejected material, ε bind , and the specific 'kick' energy, ε kick (see Kremer et al. 2022).Manukian et al. (2013) estimated the potentially unbinding 'kick' received by a star undergoing PTDE due to an SMBH from asymmetric mass loss and observed that it depends solely on the mass loss and not on the star-SMBH mass ratio.However, this is not true for µTDEs -the mass ratio determines the kick, and hence the boundedness, of the star-SBH system (see also Figure 2 of Kıroglu et al. 2023 for TDEs due to IMBHs).
Figure 8 shows the post-disruption specific orbital energy of the star-BH system, ε orb , for the disruption of a 1 M ⊙ and a 0.5 M ⊙ star by a 10 M ⊙ (top panel) and a 40 M ⊙ (bottom panel) BH.The values are normalized to the square of σ = 10 km s −1 , the typical velocity dispersion of Milky Way globular clusters (e.g., Baumgardt & Hilker 2018).More negative (positive) values represent more bound (unbound) orbits.In general, ε orb be-comes more negative with decreasing b.If the density concentration factor ρ conc is high enough (as in the case of a TAMS 1 M ⊙ star; see Figure 4), this trend continues till FTDE.However, in the case of a ZAMS 1 M ⊙ star (and less extremely in a MAMS star), there is a certain value of b at which ε orb starts to increase with decreasing b.This is due to the momentum kick from asymmetric mass loss, which becomes significant when the mass loss is high.In the case of a PTDE of a 1 M ⊙ ZAMS star by a 40 M ⊙ BH, the post-disruption orbit can be hyperbolic (unbound) for b values close to FTDE.This trend is more extreme for a 0.5 M ⊙ star, with hyperbolic orbits being possible even in the case of PTDE by a 10 M ⊙ BH.Moreover, their velocities at infinity are high enough (∼ 100-300 km s −1 ) to escape the cluster entirely.The two different BH masses result in quite different quantitative orbital parameters.Similar to Kremer et al. 2022, we can write each of the energy components as follows: Here, T 2 and T 3 are terms which depend on the stellar mass, radius, and density structure, BH mass, and distance of approach (Press & Teukolsky 1977;Lee & Ostriker 1986).For the range of parameters considered, with m BH > m ⋆ and r p ∼ r ⋆ , the values of T 2 and T 3 lie in the range 0.01 -0.1.These values are systematically lower for stars with lower ρ conc values.As a result, in 0.5 M ⊙ stars (and ZAMS 1 M ⊙ stars), ε tide is dominated by ε kick , resulting in them becoming unbound from the BH.It should be noted that ε tide depends very strongly on the distance of approach.
In addition, we find that the change in the specific orbital angular momentum of the star-BH system, h orb , is not significant.For orbits very close to being parabolic, which is the case for all our post-disruption orbits (see Figure 9), h orb ≈ 2G(m ⋆ + m BH )r p ≈ 2Gm BH r p .This implies that the postdisruption periapsis distance r p = a(1 − e), where a is the semimajor axis and e is the eccentricity, is nearly the same as the initial periapsis distance r p ≡ br t , which is also seen in PT-DEs by SMBHs (Ryu et al. 2020c).Consequently, given that a = −0.5 G(m ⋆ + m BH )/ε orb , the calculation of e is straightforward.Figure 9 shows the post-disruption e, computed using Equations 1 and 2, as a function of b which follows the trend in ε orb .The dashed line represents parabolic orbits and separates the parameter space into bound (eccentric) and unbound (hyperbolic) orbits.When mass loss is relatively low, a lower value of b generally results in more negative ε orb in the case of 1 M ⊙ star.This trend continues for the TAMS star.However, for the ZAMS and MAMS stars, significant mass losses can reverse this trend, with hyperbolic orbits also being a possibility (e.g., 1 M ⊙ ZAMS star and 40 M ⊙ BH mass with b < 1.00).0.5 M ⊙ stars, on the other hand, can become unbound for larger b, with boundedness also depending on the mass of the BH.

Rates of tidal encounters
The rate of µTDEs per year per Milky Way-like galaxy (similar to Perets et al. 2016) can be expressed as: Here, N clus ∼ 150 MWGal −1 is the number of globular clusters (in the Milky Way, e.g., Harris 2010), N ⋆,c is the number of stars in a cluster core, N BH,c is the number of BHs in a cluster core, Σ = πr 2 p (1 + v 2 p /σ 2 ) is the gravitational focused encounter cross-section with periapsis distance r p = br t (see Section 2) and periapsis velocity v p = G(m BH + m ⋆ )/r p , and σ is the velocity dispersion of the cluster.For globular clusters in the Milky Way Galaxy, the typical values of these parameters are n ⋆,c = N ⋆,c V −1 c ∼ 10 2 -10 7 pc −3 (e.g., Baumgardt & Hilker 2018), N BH,c ∼ 10-100 (e.g., Morscher et al. 2015;Askar et al. 2018), σ ∼ 1-20 km s −1 (e.g., Baumgardt & Hilker 2018).With these values, and assuming v p >> σ in globular clusters, and typical star and BH masses of 1 M ⊙ and 10 M ⊙ respectively, the rate is: This equation implies that PTDEs (with larger r p /r t ) are more frequent than FTDEs.More specifically, from our simulations of TDEs of 1 M ⊙ stars by 10 M ⊙ , encounters with b ∼ 1-2 would result in eccentric bound orbits, which can circularize over time through tides.A more detailed rate calculation would involve integrating Equation 11 over a range of star and black hole masses.The above rates are derived for µTDEs produced in dense stellar systems such as globular clusters and nuclear star clusters.However, other channels also contribute to µTDEs (see discussion in Perets et al. 2016), such as ultra-wide binaries and triples in the field (Michaely & Perets 2016, 2020) and post-natal-kicks leading to close encounter between a newly formed BH and a stellar companion (Michaely et al. 2016).These channels may give rise µTDEs in young environments and in the field.

Summary of fits
Below we present the post-disruption fits detailed in previous sections for use in globular cluster codes.m ⋆→BH refers to the gas mass bound to the black hole after disruption, including the accreted mass.

Masses of remnant star and BH
4.6.2.Spin of remnant star

Validity of fits
It should be noted that the density concentration factor, ρ conc , only captures the full disruption of a star.Similarly defined factors for each radius r, i.e. (ρ(r)/ρ) 1/3 , are required to calculate the distance at which mass beyond r is lost (see Ryu et al. 2020b).
Stars of higher masses tend to have much steeper density profiles when compared to 1 M ⊙ stars, thus Equation 4, which fits the mass loss, is not universal.In particular, the mass loss fit for 2 M ⊙ stars would involve an exponent of ∼ 3, unlike our derived exponent of ∼ 2.1 in Equation 4. Similar arguments follow for the fits of spins and orbits.

Implications
In this section, we discuss the general observational implications and the caveats of our study.µTDEs can potentially produce both unique and peculiar transient events as well as longer-lived systems.

Low-mass X-ray binaries
Low-mass X-ray binaries (LMXBs) are interacting binary systems which consist of compact objects (neutrons stars or black holes) accreting mass from their low-mass (≲ 1 M ⊙ ) stellar companions, thereby producing X-rays.A few hundred of LMXBs have been detected in the Milky Way Galaxy (e.g., Liu et al. 2007;Avakyan et al. 2023), highlighting the prevalence of such systems.However, theoretical (e.g., Portegies Zwart et al. 1997;Kalogera 1999) and population synthesis (e.g., Podsiadlowski et al. 2003;Kiel & Hurley 2006) studies have not been able to match the rates of observed LMBXs through isolated binary evolution alone.The primary reason is that a low-mass companion has insufficient energy to expel a common envelope during the giant phase of the primary star, thereby resulting in the spiraling in of the companion and an eventual merger.Thus, a dynamical channel formation channel for LMXBs, where a low-mass star is tidally captured by a black hole in a globular/nuclear cluster, holds promise (e.g., Voss & Gilfanov 2007;Michaely & Perets 2016;Klencki et al. 2017).Some of our simulations of encounters of low-mass stars with black holes are not close enough for significant mass loss and can result in highly eccentric bound orbits.For example, in the parabolic encounter of a 1 M ⊙ star and a 10 M ⊙ BH, for b ∼ 1.0-1.5 (see Section 4), the fractional mass loss is ∼ 1-10%, and the post-disruption bound orbit has a period of 0.1-1 yr and an eccentricity of ∼ 0.97-0.98.However, tidal dissipation at subsequent periapsis passages can reduce the eccentricity and orbital period (e.g., Klencki et al. 2017) such that the resulting bound system evolves into a compact LMXB.
Another type of system that is potentially formed through dynamical encounters are non-interacting BH-star systems, e.g., El-Badry et al. (2023a,b); Chakrabarti et al. (2023).These systems consist of ∼ 1 M ⊙ stars orbiting ∼ 10 M ⊙ black holes in moderately-eccentric orbits with periods of ∼ 1 yr.These parameters correspond quite well with our simulations with b ∼ 1.0-1.5, which lends credence to dynamical formation.

Intermediate-mass black hole growth
Intermediate-mass black holes (IMBHs; see Greene et al. 2020 for a review), with masses greater than those of stellar-mass black holes (≳ 10 2 M ⊙ ) and less than those of supermassive black holes (≲ 10 6 M ⊙ ), are expected to exist in the centers of dwarf galaxies (e.g., Kunth et al. 1987;Filippenko & Sargent 1989;Reines 2022;Gültekin et al. 2022) and globular clusters (e.g., Davis et al. 2011;Farrell et al. 2014;Pechetti et al. 2022).Among several IMBH formation mechanisms (see Volonteri et al. 2021 for review), one possible avenue is runaway mass growth through tidal capture and disruption, which has been examined analytically (e.g., Stone et al. 2017) and numerically (e.g., Rizzuto et al. 2023;Arca Sedda et al. 2023).One of the most critical factors in this mechanism for BH growth is the amount of debris mass that ultimately accretes onto the BH.For example, Rizzuto et al. (2023) assume that a star approaching a BH within the nominal tidal radius r t is fully destroyed and 50% of the stellar mass is accreted onto the BH.
While these assumptions may be valid for a part of the parameter space, it is important to develop a prescription for the outcomes of TDEs that works over a broader parameter range to more accurately assess the possibility of BH mass growth through TDEs.Our simulations can provide such prescriptions that can improve the treatment of TDEs.Firstly, our simulations show, like previous work on TDEs by supermassive black holes (e.g., Guillochon & Ramirez-Ruiz 2013;Mainetti et al. 2017;Goicovic et al. 2019;Ryu et al. 2020a,b,c,d;Law-Smith et al. 2020), that the periapis distances at which stars are fully or partially destroyed depend on the stellar internal structure.Ryu et al. (2020b) analytically demonstrated that the nominal tidal radius r t is a more relevant quantity for partial disruption events involving stars with m ⋆ ≳ 1 M ⊙ .Only for lower-mass stars, where ρ conc ∼ 1, does r t come close to the genuine full disruption radius4 .Our fitting formula for the fractional mass loss, Equation 4, can provide a better prescription for determining the fate of the star (full or partial disruption) as well as the mass of the remnant.
The amount of debris accreted onto the BH remains highly uncertain.To the zeroth order, if the accretion rate is super-Eddington, the strong radiation pressure gradient would drive strong outflows, hindering continuous and steady accretion (S ądowski et al. 2014).However, the accretion efficiency would be collectively affected by many factors, including magnetic fields, black hole spins, accretion flow structure, and jet formation (e.g., S ądowski et al. 2014;Jiang et al. 2019;Curd & Narayan 2023;Kaaz et al. 2023).Our fitting formulae for the fractional mass loss cannot provide an accurate prescription for the accreted mass but can place constraints on the maximum mass that can be accreted onto the BH.

Fast blue optical transients
Fast blue optical transients (FBOTs) are a class of optical transients characterized by high peak luminosities > 10 43 erg s −1 , rapid rise and decay times of the order of a few days, blue colors (Perets et al. 2011;Drout et al. 2014), and peak blackbody temperatures of a few 10 4 K.Although a majority of these events can be explained by supernovae with low-mass ejecta (Pursiainen et al. 2018) Chrimes et al. 2024), are too bright at their peaks and fade too rapidly to be explained by supernovae.Proposed mechanisms involve compact objects, such as black hole accretion or magnetars formed in core-collapse supernovae (Prentice et al. 2018), mergers between Wolf-Rayet stars and compact objects accompanied by hyper-accretion (Metzger 2022), or tidal dis-ruption events by intermediate-mass BHs (Kuin et al. 2019) or stellar-mass BHs (Kremer et al. 2021).
Some luminous FBOTs also emit in radio or X-ray (e.g., Margutti et al. 2019;Ho et al. 2020) which indicates the presence of a circumstellar medium that may result from partial TDEs (Kremer et al. 2021).Our simulations support this notion by showing the presence of debris near the BH with occasional weak outflows from the BH.In addition, the remnants in our simulations are often bound to the BH and can subsequently undergo multiple partial disruptions, e.g., our simulations of TAMS 1 M ⊙ stars and 10 M ⊙ BHs for b = 0.25 or b = 0.33.While we cannot confirm a direct correlation between FBOTs and µTDEs, multiple disruption events of a star may fulfill some of the astrophysical conditions (e.g., the presence of a gas medium and rapid accretion onto a compact object) necessary to explain FBOTs with radio emission.For such cases, the detection of repeated bursts in FBOTs could provide valuable constraints on their formation mechanisms.

Ultra-long gamma-ray bursts and X-ray flares
The exact appearance of transients resulting from µTDEs is highly uncertain and depends on various assumptions.Perets et al. (2016) suggested that, if the accretion onto the BH is efficient and gives rise to jets, there is a possibility of producing ultra-long gamma-ray bursts (GRBs) and/or X-ray flares (XRFs).
Assuming a proportional relation between the material accreted to the compact object and the luminosity produced by the jet, we can generally divide µTDEs into two regimes: (1) t min ≫ t acc and (2) t min ≪ t acc , where t min is the typical fallback time of the debris following the disruption, and t acc is the typical viscous time of the accretion disk that can form around the BH (see Perets et al. 2016 for more details).
t min ≫ t acc : When the mass of the star is much smaller than that of the BH, the accretion evolution is dominated by the fallback rate, i.e., the light curve should generally follow the regular power-law (e.g.t −5/3 power-law for a full disruption).t min ≪ t acc : When the masses of the star and the BH are comparable, the fallback material is expected to accumulate and form a disk on the fallback time, which then drains on the longer viscous time maintaining a low accretion rate.We would expect the flaring to begin only once the material is accreted onto the compact object.Therefore, we expect four stages in the light curve evolution: (1) A fast rise of the accretion flare once the disk material is processed and evolves to accrete on the compact object.
(2) Accretion from the disk until the accumulated early fallback material is drained; if we assume a steady state accretion until drainage, one might expect a relatively flat light curve.(3) Once the disk drains the accumulated early fall-back material the light curve should drop steeply.(4) The continuous fall-back of material would govern the accretion rate at times longer than the viscous time, and the accretion rate should then follow the t −5/3 rate (or a different power-law).The exact light curve at the early stages is difficult to predict, but we do note a non-trivial signature of the µTDEs relating to the early and late stages.
The expected properties of µTDEs are consistent with and might explain the origins of ultra-long GRBs (Levan et al. 2014): long-lived (∼ 10 4 s) GRBs which show an initial plateau followed by a rapid decay.µTDEs may also explain the origin of some Swift TDE candidates (Bloom et al. 2011;Burrows et al. 2011;Cenko et al. 2012;Brown et al. 2015), suggested to be produced through a TDE by a supermassive BH.The typical timescales for the latter are longer than the observed 10 5 s, challenging the currently suggested origin, but quite consistent with a µTDE scenario.
µTDEs producing ultra-long GRBs are also expected to produce afterglow emission if/when the strong outflows and jets launched by accretion propagate into the surrounding medium and the resulting shock interaction produces high energy particles which give rise to synchrotron and inverse Compton radiation.Such observational signature of µTDE has been little explored, though analogous observation and modeling have been suggested for TDEs by massive black holes (e.g., Giannios & Metzger 2011;Zauderer et al. 2011).Detection of afterglows from µTDEs would provide information on the energetics and dynamics of the outflow, robust independent evidence for the jet collimation, as well as identification of the host galaxy.

Remnants
Partial µTDEs of stars by stellar BHs can spin up the stars significantly due to strong tidal interactions (see Fig. 7 and Alexander & Kumar 2001).Though magnetic breaking can potentially spin down such stars over time, observations of highly spun-up low-mass old stars in globular clusters could provide potential signatures for partial µTDEs.In addition, the remnant undergoes violent chemical mixing during the first pericenter passage and has higher entropy than an ordinary star of the same mass and age (Ryu et al. 2020c(Ryu et al. , 2023a)).If the mass loss is significant, a unique chemical composition profile may persist even after the remnant has relaxed to a stable state.This suggests that the remnant could exhibit unique asteroseismic signatures that distinguish it from ordinary stars (Bellinger et al. 2023, in prep).

Caveats
As mentioned in Section 2, we exclude the effects of relativity and magnetic fields in our study.The Newtonian approximation is justified for TDEs by SBHs because most of the debris stays outside r p , which is approximately four orders of magnitude greater than the gravitational radius r g ≡ Gm BH /c 2 5 .However, relativistic effects and proper treatment of radiation are more important for the accretion flow, which is beyond the scope of our paper.The effect of magnetic fields on the dynamics of the remnant is unclear, which we will leave for our future work.A final point to note is that, unlike globular clusters, nuclear star clusters have typical velocity dispersions of σ ≳ 100 km s −1 Figer et al. 2003;Schödel et al. 2009) depending on the galaxy, and subsequently, |1 − e| ≳ 10 −2 .Thus, encounters are not always parabolic, and for a more thorough study, hyperbolic encounters must be taken into account.

Summary and conclusion
We performed a grid of 58 hydrodynamics simulations, using the moving-mesh code AREPO, of partial (PTDE) and full (FTDE) tidal disruptions of main-sequence (MS) stars with stellar-mass black holes (SBH) on initially-parabolic orbits.Our varied parameters include stellar mass m ⋆ (0.5 M ⊙ and 1.0 M ⊙ ), SBH mass (10 M ⊙ and 40 M ⊙ ), core hydrogen abundance H c (a proxy for MS age t ⋆ ) and impact parameter b ≡ r p /r t .Our stellar models are initialized in accordance with 1D detailed non-rotating stellar profiles from MESA.We also define a star's density concentration parameter ρ conc ≡ (ρ c /ρ) 1/3 .Our main results are summarized below: -The mass of a post-disruption stellar remnant m ⋄ decreases with decreasing b.A higher ρ conc (more centrally concentrated) results in a lesser mass loss for the same b.m ⋄ depends only on b and ρ conc , and not on m BH .Roughly half of the disrupted mass remains bound to the black hole, while the other half is unbound from the system.-The spin angular momentum of a stellar remnant L ⋄ increases after periapsis passage.For large b (when the mass loss is small), L ⋄ increases with decreasing b, but this trend reverses for very small b (close to FTDE) due to significant mass loss.This reversal occurs for lower b when ρ conc is higher.The final spin is also independent of m BH .For low b values, the rotational velocities can be very close to breakup velocities.-Unlike mass and spin, the orbit of the remnant star also depends on m BH .For large b, eccentricity e and semimajor axis a decrease with decreasing b keeping a constant periapsis distance, i.e., the remnant tends to be more 'bound'.Depending on m BH and ρ conc , this no longer holds for approaches close to FTDEe can increase and the remnant can become unbound.In general, a higher m BH /m ⋆ ratio and a lower ρ conc tend to be unbinding.-We provide relatively simple and accurate fitting formulae for the masses and spins of the remnant stars.These analytical fits are listed in Section 4.6 for easy access.
We discussed the implications that µTDEs can have on lowmass X-ray binary formation, intermediate-mass black hole formation, and transient searches.Given the limited range of stellar masses and ages that we covered in this work, we will extend our suite of simulations in the future to cover a wider parameter space, including massive stars.It is necessary to explore this regime to model stellar dynamics in dense stellar clusters.

Fig. 2 :
Fig.2: Similar profiles as Figure1, but for a 0.5 M ⊙ MS star.A 0.5 M ⊙ star of age ∼ 13.5 Gyr with H c = 0.60 (purple curves) is very similar to a ZAMS star with H c = 0.70 (red curves), albeit with a slightly larger radius.

Fig. 3 :
Fig. 3: Grid of log density slices of stars ∼ 20 t dyn,⋆ after undergoing µTDEs, with the BH at the center of each panel.Each row is for a different m BH (top and bottom half for 10 M ⊙ and 40 M ⊙ respectively), m ⋆ or H c , whereas the columns represent increasing b(from left to right).The dashed and solid circles around the BH, most visible in the rightmost panels, denote the tidal radius, r t , and the periapsis distance, r p , respectively.Note that the spatial range of each panel is not the same.

Fig. 4 :Fig. 5 :
Fig.4: MESA-derived inverses of density concentration factors ρ −1 conc for a 1 M ⊙ MS star as a function of stellar age.The three dashed vertical lines represent the ages of three chosen MS models in our study, corresponding to the profiles in Figure1and the values in Table1.The dependence with age is very close to linear.

Fig. 6 :
Fig.6: Post-disruption fractional mass loss, ∆m ⋄ /m ⋆ , for a star of mass 1 M ⊙ (left) and 0.5 M ⊙ (right), due to a BH of mass 10 M ⊙ (top) and 40 M ⊙ (bottom), as a function of impact parameter b and MS (H c abundance).The dashed line at the top signifies the full disruption of the star.The fractional mass losses are largely independent of BH mass.Mass loss is higher for lower b values, and this decreases roughly exponentially (best-fit lines and standard error regions shown) with increasing b.A less evolved star (higher H c ) loses more mass for the same value of b.A 0.5 M ⊙ star has a higher mass loss than a 1 M ⊙ star for the same b.

Fig. 7 :
Fig. 7: Post-disruption spin angular momentum, L ⋄ , for an initially non-rotating star of mass 1 M ⊙ (left) and 0.5 M ⊙ (right), due to a BH of mass 10 M ⊙ (top) and 40 M ⊙ (bottom), as a function of impact parameter b and MS age (H c abundance).The dashed lines indicate the (approximate) estimated break-up angular momenta of the remnants L ⋄,break respectively.The spin angular momenta are largely independent of BH mass.It decreases almost exponentially for sufficiently high b values, and it drops for very low b values when the remnant mass is small (best-fit lines and standard error regions shown).The drop-off of angular momentum occurs at larger b for a less evolved star.

Fig. 8 :
Fig.8: Post-disruption normalized specific orbital energies of the star-BH system, ε orb , for a star of mass 1 M ⊙ (left) and 0.5 M ⊙ (right) initially in a parabolic orbit, due to a BH of mass 10 M ⊙ (top) and 40 M ⊙ (bottom), as a function of impact parameter b and MS age (H c abundance).The dashed line represents parabolic orbits and separates the parameter space into bound (eccentric) and unbound (hyperbolic) orbits.When mass loss is relatively low, a lower value of b generally results in more negative ε orb in the case of 1 M ⊙ star.This trend continues for the TAMS star.However, for the ZAMS and MAMS stars, significant mass losses can reverse this trend, with hyperbolic orbits also being a possibility (e.g., 1 M ⊙ ZAMS star and 40 M ⊙ BH mass with b < 1.00).0.5 M ⊙ stars, on the other hand, can become unbound for larger b, with boundedness also depending on the mass of the BH.

Fig. 9 :
Fig.9: Post-disruption orbital eccentricities, e, for a star of mass 1 M ⊙ (left) and 0.5 M ⊙ (right) initially in a parabolic orbit, due to a BH of mass 10 M ⊙ (top) and 40 M ⊙ (bottom), as a function of impact parameter b and MS age (H c abundance).The dashed line represents parabolic orbits and separates the parameter space into bound (eccentric) and unbound (hyperbolic) orbits.The reasoning for the trends is the same as that for Figure8.

Table 2 :
Initial parameters of our suite of simulations.

Table 3 :
Evaluated values of the impact parameter for the full disruption of a star with uniform density, bFTDE = b FTDE /ρ conc .