Partial tidal disruption events: The elixir of life

In our Galactic Center, about 10,000 to 100,000 stars are estimated to have survived tidal disruption events, resulting in partially disrupted remnants. These events occur when a supermassive black hole (SMBH) tidally interacts with a star, but not enough to completely disrupt the star. We use the 1D stellar evolution code Kepler and the 3D smoothed particle hydrodynamics code Phantom to model the tidal disruption of 1, 3, and 10 solar mass stars at zero-age (ZAMS), middle-age (MAMS), and terminal-age main-sequence (TAMS). We map the disruption remnants into Kepler in order to understand their post-distribution evolution. We find distinct characteristics in the remnants, including increased radius, rapid core rotation, and differential rotation in the envelope. The remnants undergo composition mixing that affects their stellar evolution. Whereas the remnants formed by disruption of ZAMS models evolve similarly to unperturbed models of the same mass, for MAMS and TAMS stars, the remnants have higher luminosity and effective temperature. Potential observational signatures include peculiarities in nitrogen and carbon abundances, higher luminosity, rapid rotation, faster evolution, and unique tracks in the Hertzsprung-Russell diagram


INTRODUCTION
argued that the Galactic Centre of the Milky Way harbours about 10 4 -10 5 stars that survived a tidal disruption event (TDE) with the supermasive black hole (SMBH) at its centre.Manukian et al. (2013) estimated a similar number of stars.The upcoming Vera C. Rubin Legacy Survey of Space and Time will result in an explosion in the detection of partial TDEs (Hambleton et al. 2023).Understanding partial TDEs hence is essential.
Most of the stars in our NSC are old (age > 5 Gyr), late-type giant stars and helium-burning stars on the horizontal branch with mass between 0.5 M ⊙ to 4 M ⊙ .A few asymptotic giant branch stars with temperatures less than 2,700 K along with a few red supergiant stars have also been observed.About 200 young (3-8 Myr) Wolf-Rayet, O and B type stars have been detected within 0.5 pc.Possible explanations of their existence include short events or weak starbursts about 6 Myr ago or the migration of a star cluster to the central parsec (see reviews by Genzel et al. 2010;Neumayer et al. 2020).Within 0.04 pc, a cluster of B-stars orbiting Sgr A* on highly eccentric and inclined orbits have been detected, called the S-star cluster.Their existence in such proximity to the SMBH is puzzling, with proposed explanations including the Hills mechanism (Hills 1988) and migration (Genzel et al. 2010).About six G-objects -having characteristics associated with dust and gas clouds, but dynamics similar to stellar-mass objects (Ciurlo et al. 2020), have also been detected within 0.04 pc (Gillessen et al. 2012;Phifer et al. 2013;Witzel et al. 2017;Ciurlo et al. 2020).The cause of their formation remains a mystery.
The dynamics of the stars present in NSC are influenced by the SMBH's radius of influence ( h =   • / 2 ), where  • is the mass of SMBH and  is the velocity dispersion.Gravitational relaxation processes -perturbations by massive objects (Perets et al. 2007), and/or two-body stellar relaxation interactions (Magorrian & Tremaine 1999) result in stars entering the loss cone.The loss cone is defined as a region in the phase space where the orbits have the closest approach to the black hole (pericentre;  p ) less than the tidal radius ( t ) given by  t = 0.47 au where  * is the stellar mass and  * is the radius of the star (Magorrian & Tremaine 1999).The fate of these stars can include tidal disruption by the SMBH, direct plunge into SMBH, capture by SMBH and inspiral, or ejection (Hills 1975;Alexander 2017).Not all stars that enter the loss cone are disrupted by the SMBH, as they can move in and out of the loss cone several times per orbit (Magorrian & Tremaine 1999).Direct -body simulations suggest that almost all stars entering the tidal radius have eccentricities close to unity, with about half of the stars having  > 1, implying marginally hyperbolic orbits (Zhong et al. 2014).Alexander (2017) argued that the stars scattered into the loss cone are on hyperbolic orbits with the specific orbital energy a small fraction of the binding energy.Hence, these orbits can be approximated as parabolic orbits.Hayasaki et al. (2018) found, using -body simulations, that stars that are disrupted by the SMBH are rarely on an eccentric, precisely parabolic or hyperbolic orbit.Most of the particles followed marginally eccentric or marginally hyperbolic orbits.Zhong et al. (2023) also reported a similar result.
We focus on partial tidal disruption events (TDEs) in this paper.A star is tidally disrupted when it approaches the tidal radius, a distance at which the gravitational forces of the SMBH can overcome the self-interaction forces of a star (Hills 1975;Rees 1988;Evans & Kochanek 1989).TDEs were theorised in the 1970s (Hills 1975) and first observationally detected in the 1990s by the ROSAT All Sky Survey in X-rays (Donley et al. 2002).
During a TDE, a star can either undergo a full disruption or leave behind a remnant -a partial TDE.The strength of a TDE is quantified by the penetration factor  ≡  t / p (Lodato et al. 2009;Guillochon & Ramirez-Ruiz 2013).The boundary between full and partial disruption remains uncertain due to several factors such as the star's structure, SMBH spin and general relativistic effects (Gezari 2021).
A few dozen TDEs have been detected in UV (Martin et al. 2005), optical (van Velzen et al. 2011), soft and hard X-rays since the first detection (see Gezari 2021 and references within).A TDE rate of 10 −4 -10 −5 per galaxy per year (Magorrian & Tremaine 1999;Brockamp et al. 2011;Stone & Metzger 2016) has been estimated but this can be influenced by the spin of the SMBH, the nature of the galaxy (e.g.post-starburst galaxies have a higher rate Stone & van Velzen 2016) and the presence of binary SMBH (see review by (Gezari 2021)).In the past few years, a few partial TDEs have been inferred based on the decay rate of their lightcurves (Blagorodnova et al. 2017;Nicholl et al. 2020).Recently, repeating partial TDEs have also been observed (Payne et al. 2021;Miniutti et al. 2023;Huang et al. 2023;Somalwar et al. 2023;Wevers et al. 2023).Stone et al. (2020) argued that partial TDEs should be more common than full disruptions due to the larger cross-section.Chen & Shen (2021) found that partial TDEs become dominant for  • ≥ 10 6 M ⊙ .Using -body simulations, Zhong et al. (2022) suggested that partial TDEs are 75% more likely than analytical estimates based on the encounter cross-section, as a single star can undergo repeated partial TDEs.Recently, Bortolas et al. (2023) found that partial TDEs are more likely than full TDEs by a factor of ∼ 10.
Previous studies have delved into understanding the properties of remnants formed during partial TDEs.Alexander & Livio (2001) analytically explored remnants of a  = 1.5 polytrope and a solartype star.Both models were initially on slightly hyperbolic orbits in Newtonian gravity.They argued that stars would experience spinup from the SMBH resulting in rotationally induced mixing.The remnants would be bluer and more luminous due to the stripping of material compared with a main-sequence star of the same mass.Hydrodynamical simulations by Antonini et al. (2011) indicate that stars bound to SMBH can undergo multiple TDEs and result in depletion of the outermost stellar region, leaving behind a hot central core.Manukian et al. (2013) used simulations of partial disruptions of solar-type stars on parabolic orbits, finding that stars can have velocities of several hundred km s −1 , preventing them from escaping the galactic centre.Goicovic et al. (2019) simulated the tidal disruption of a 1 M ⊙ zero-age main-sequence (ZAMS) star and found that the surviving core has positive angular momentum in outer layers, while negative in the innermost region.
Guillochon & Ramirez-Ruiz ( 2013) explored the boundary between partial and full TDEs for polytropic stellar models of  = 4 /3 and  = 5 /3 using an adaptive-mesh grid-based hydrodynamics code.They found a critical  of 1.85 for  = 4 /3 and 0.9 for  = 5 /3.Goicovic et al. (2019) determined a critical  ∼ 2.5 for their 1 M ⊙ ZAMS model.Law-Smith et al. (2020) performed a detailed study with a range of stellar masses for stellar models at different ages.They found that critical  is related to the mean ratio of the central and mean density of a star ( c / ρ).Ryu et al. (2020c) utilised a fully general relativistic framework for a range of stellar masses at middle age main-sequence (MAMS), discovering remnants with velocities of few hundred km/s for strong disruptions (mass of remnant ≈ 40 to 60% of the initial star mass).Bound remnants were all on highly eccentric orbits.They found remnants of angular frequencies close to break-up in the outer layers.
In this paper, we study how the stellar evolution of disrupted remnants proceeds post disruption.To address this question, we evolve stars in a 1D stellar evolution code, Kepler (Weaver et al. 1978), disrupt them using a 3D smoothed particle hydrodynamics code, Phantom (Price et al. 2018), considering general relativistic effects (Liptai & Price 2019) with Kerr metric (Kerr 1963) of zero spin, and map them back into 1D to continue their lives.We explore the disruption of stars with masses of 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ at three different stages of their lives for a range of penetration factors.A key finding is that disrupted stars are 'forever young', staying on the main-sequence for longer than the age of the Universe.For stars, partial disruptions are the long-sought Holy Grail, the fountain of youth and the elixir of life.
This paper is structured as follows: In Section 2 we describe the method used to simulate TDEs and map the resulting remnants back into Kepler to continue their evolution.In Section 3 we present our results.In Section 4 we discuss our results and compare them with current literature.We summarise in Section 5.

METHODS
We used the general relativistic smoothed particle hydrodynamics code Phantom (Price et al. 2018;Liptai & Price 2019), and the stellar evolution code Kepler (Weaver et al. 1978) for our simulations.We mapped Kepler models into Phantom to simulate the tidal disruption, and then remapped into Kepler to continue the evolution.

Mapping from 1D to 3D
We computed stellar models using Kepler for stars with masses of 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ at three different stages of the star's lifecycle.Table 1 lists the parameters of our models.These stages correspond to the Zero-Age Main-Sequence (ZAMS), Middle-Age Main-Sequence (MAMS), and Terminal-Age Main-Sequence (TAMS) of the star.ZAMS, MAMS, and TAMS are defined as the stellar evolution stages when the mass fraction of hydrogen in the core has dropped by 1% (absolute), is half the initial value, and when it has reached 1% (absolute), respectively.Figure 1 shows the evolutionary paths of these Kepler models.We also show the three stages at which the models were analysed.As the star evolves on the main-sequence, the luminosity increases while the effective temperature decreases.Our Kepler models use an approximate nuclear network, APPROX19, that comprises neutrons, 1 H, 3 He, 4 He, 12 C, 14 N, 16 O, 20 Ne, 24 Mg, 28 Si, 32 S, 36 Ar, 40 Ca, 44 Ti, 48 Cr, 52 Fe, 54 Fe, and 56 Ni.This network is adequate to capture the main nuclear burning phases of massive stars and their nuclear energy generation.
To obtain a 3D density profile from the initial 1D density data, we used the stellar profile mapping and relaxation procedure outlined in Lau et al. (2022).We fixed the radial entropy profile of the star during the relaxation.We stopped the relaxation procedure when the ratio of kinetic to potential energy dropped below 10 −7 and the error in density was less than 1%.All stellar models were resolved with 10 6 SPH particles.We employed a specific heat ratio of  = 5 /3 and used an adiabatic equation of state to calculate the energy and temperature from the relaxed pressure and density profiles.
The composition data obtained from Kepler was then interpolated onto the SPH particle positions and saved into a file.We did not simulate nuclear reactions in Phantom, as the simulation duration was significantly shorter than the nuclear-burning timescales even during the peak of the compression.

Disrupting stars in Phantom
Once the star was relaxed, we placed it onto an orbit corresponding to a Newtonian parabolic orbit around a 10 6 M ⊙ SMBH at a distance ∼ 10 times the tidal radius (Eq.1).We later realised that this procedure resulted in slightly hyperbolic orbits due to the difference between the general relativistic and Newtonian energy at 10 r t .While this difference results in a faster velocity of the remnant at infinity, the effect on the tidal disruption of the star is not significant (see Appendix D for a comparison).
We used the following conditions in our simulations in Phantom.A Courant number of 0.3 was used which is also the default value.As we used equal mass particles, we used a proportionality factor of 1.2 that specified the smoothing length in terms of mean local particle spacing.To use the smoothing length and density interchangeably during the time-stepping, we used a tolerance of 10 −4 which is the default value.A tolerance of 10 −7 was used for both position and momentum iterations in the GR timestep integrator (Liptai & Price 2019).We also used a tree accuracy of 0.5 in the self-gravity calculation.We used high-resolution shock-capturing dissipation in our simulations, with the shock viscosity factors  = 1 and  = 2.We also used a shock conductivity of 0.1.We used an adiabatic equation of state with a mean molecular weight of 0.5988 for all the models.Heating from  d work and shocks were assumed to be trapped rather than radiated.We used  = /  as the energy variable (Liptai & Price 2019).We performed simulations in the Kerr (1963) metric with an SMBH with no spin (Liptai & Price 2019).
We performed the initial disruption with an accretion radius set to the last stable orbit at 6  / 2 inside of which particles were deleted from the simulation.After the initial disruption, for computational efficiency, we manually reset the accretion radius such that the material falling onto the SMBH would be removed but not so large that it would remove the material falling back onto the remnant (about 100-1,000  / 2 depending on the simulation).
We assumed that the mass of the accretion disk formed after disruption is negligible compared to the mass of the SMBH and the particles within the accretion radius are accreted by the SMBH.
Table 2 lists the  values, pericentre distance of our simulations.

Mapping back: 3D to 1D
To model the post-disruption stellar evolution, we mapped the Phantom models back into Kepler.To project the 3D data from Phantom back into 1D Kepler models, we binned the star by radius.

What is part of a remnant?
We used the point of maximum density as the centre of the star instead of using the centre of mass of all the particles because the removed particles/debris could move to very large distances, which could result in an inaccurate centre.To determine which particles were part of the debris, we sorted particles by radius with respect to the centre and then calculated the sum of kinetic, potential and internal energy for each particle , where Φ  is the self-gravitional potential between particles calculated by Phantom,   is the velocity of the particle relative to the maximum density particle and   is the particle mass, which is fixed in Phantom.We ignored the black hole as the star was sufficiently far away at the time of analysis.We determined that a particle was bound to the remnant if the total energy was negative and the kinetic energy was less than half of the gravitational potential energy, according to the virial theorem.This ensured that only particles close to the point of central density were considered as part of the remnant.
Figure A2 shows that the remnant stars from stronger disruptions contain streams.By examining the debris temperature as a function of radius, we found that more particles were dispersed horizontally at lower temperatures.We leveraged this to distinguish particles that were part of the stream from those that were part of the hydrostatic star.See details in Appendix A.

Binning algorithm
We calculated the density, temperature, composition, angular velocity, and radial velocity of each SPH particle.To determine temperature, we first computed the mean molecular weight based on each particle's composition and then used it along with specific energy in the adiabatic equation of state to find the particle's temperature.
The radial velocity for each particle was determined using where r was taken with respect to the position of the particle with the maximum density in the remnant.The angular momentum for each particle was calculated using We sorted particles by radius and then fixed the number of bins to a constant,  bins .To achieve optimal mass resolution in the core while limiting the ratio of masses in neighbouring zones to not exceed two, we binned the particles by the following procedure.We initially placed one particle in the first bin.Moving outward in radius, we then doubled the number of SPH particles in each bin until the desired mass resolution of  p / bins •   , where   is the number of SPH particles in the remnant star, and   is the particle mass.We continued adding particles until reaching the desired maximum zone mass of  p / bins •   or the radius is increased by 30%, whichever was achieved first.
Finally, we calculated the required quantities for each bin.Density, temperature, and radial velocity were determined by averaging all particles within each bin.We set a minimum temperature of 1,000 K.
To compute angular velocity, we computed the moment of inertia tensor for each bin, accounting for variations in axes of rotation among the particles, where    were the components of the tensor I.We then obtained the angular velocity.

Kepler simulations
Kepler is a 1D hydrodynamic stellar evolution code and solves the equation of stellar structure by modelling the conservation of energy, mass, and momentum (Weaver et al. 1978).Convection, which is inherently multi-dimensional, is modelled using the Mixing Length Theory (MLT; Böhm-Vitense 1958) with a mixing length of where  is the mixing length parameter,  is the pressure,  is the density, and  is the local gravitational acceleration (Böhm-Vitense 1958).The value of  depends on the properties of the convection zone and is well determined for the surface convection zone of our sun, where we can accurately measure effective temperature, luminosity, and even structure from helioseismology.Both mixing length parameters and opacity scaling were tweaked to give the solar luminosity at the present age of the sun.As convection does not play a significant role in stars which are ≥ 1.3 M ⊙ , a value of  = 1 was used for 3 and 10 M ⊙ models.Kepler uses the Ledoux criterion for convective stability, as well as semiconvection and thermohaline convection (Heger et al. 2005) When mapping models back into Kepler after the binning procedure implemented in Phantom, we first let the model relax until the radius of the star was less than the initial model radius, and the timestep was ≳ 1 year (specifically we chose Δ > 2 × 10 7 s).We did not perform this step if the initial radius of the model was less than the radius of the Kepler model that was disrupted.Whenever timestep exceeded 10 6 s, we allowed fine rezoning of the outer layers to re-establish an atmosphere that can not be modelled in Phantom, as it does not have sufficient resolution.
As the models were not in thermal equilibrium in Kepler after mapping, we estimated an approximate Kelvin-Helmholtz timescale for the models to settle into thermal equilibrium, using approximate main-sequence values for Figure C1 shows the evolution of a 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ ZAMS models that were mapped into Phantom, and then mapped back into Kepler.The evolution is similar to the original Kepler models, hence, showing the accuracy of our mapping back procedure.We note that for a few models, we had to ignore the outermost five bins to successfully map them into Kepler as otherwise, these models crashed.These models were 1 M ⊙ TAMS,  = 1.08; 3 M ⊙ ZAMS,  = 1.29; 10 M ⊙ ZAMS,  = 0.99; 10 M ⊙ MAMS,  = 0.98 and  = 1.95; 10 M ⊙ TAMS,  = 0.94.We also compared the remnants with Kepler models of the same mass.

What happens to the star?
Figure 2 shows snapshots of the Phantom simulations of the disruption of a 3 M ⊙ MAMS model for a range of penetration factors,  (top to bottom), showing the column density perpendicular to the orbital plane, with snapshots superimposed at every 3 h within each panel up to 8 days since the beginning of the simulation.For an initial  = 0.8 the star can be seen to have lost almost no mass, whereas for  = 2.39, we can observe that the star loses about 80% of its mass.Post-pericentre passage, we see the formation of streams around the remnants.Eventually, some of this material falls onto the SMBH.Several authors have tried to determine the power-law fits of this accreting material onto the SMBH (Guillochon & Ramirez-Ruiz 2013;Goicovic et al. 2019;Golightly et al. 2019;Ryu et al. 2020c).In this paper, we focus on the evolution of the remnants themselves.
Table 1 lists the central density and the radius of the stars we disrupted in Phantom.As a star evolves from ZAMS to TAMS, it becomes more dense and its radius increases.This affects how long it can survive the encounter.For very deep encounters, the star can be completely ripped apart, fall back, and eventually reassemble itself into a remnant.We saw this behaviour for some of our closest encounters with the SMBH.An example is the 10 M ⊙ MAMS,  = 2.15 remnant of 0.24 M ⊙ .The central density decreases for the first 14 days, but afterwards, it increases and becomes constant as a remnant is formed (see Appendix B).This model has  p = 115.80R ⊙ and exhibits complete disruption and then re-collapse.We performed simulations for different stars where they approached the SMBH at even lower  p values, but did not undergo complete disruption and then re-collapse.Guillochon & Ramirez-Ruiz (2013) and Nixon & Coughlin (2022) had used Newtonian physics but also found that the streams can re-collapse to form a remnant for their polytropic models.Nixon & Coughlin (2022) argued that this is a result of gravitational instability in the stream.

Stellar oscillations
Figure 3 shows the central density of the 3 M ⊙ MAMS remnants as a function of time for different penetration factors ().We see that the density of the remnant starts to increase as it approaches the pericentre, becoming maximum around the pericentre, and then decreasing as the star moves away from the SMBH.It eventually settles to a constant value.Even though the  = 0.5 model does not lose any mass, the star experiences a slight compression which eventually dissipates.
We find that the encounter induces oscillations in the remnant star.The period of this ringing is close to that of the fundamental mode of the star, ( osc = √︁ 3/(32 ρ) ).For example, we looked at the  = 1.79 case for 3 M ⊙ MAMS star.We considered the material bound within 2 R ⊙ as part of the remnant.This resulted in a oscillation time of 0.04 days, whereas  osc ∼ 0.07 days.We observe a similar trend in the central temperature as a function of time.
Figure 4. Absolute accretion rate onto the remnant as a function of time since the beginning of the simulation, calculated by considering the difference in mass of the remnant every hour.The grey line represents the time at which we perform analysis for our models.We selected times of 5 days, 8 days, and 10 days for 1 M ⊙ , 3 M ⊙ , and M ⊙ models, respectively for our analysis.For 10 M ⊙ ,  = 2.15 model, we chose the snapshot at 60 days.

Accretion onto the remnant
Figure 4 shows the mass accretion rate onto the remnant as a function of time, determined by calculating the difference in the mass of the remnant over 12 hours intervals.For most of the encounters shown in Figure 4, the mass accretion rate can be seen to become negligible after several days.However, we observed a persistent, albeit small, stream of material continuing to accrete over longer timescales.This stream moves away from the remnant as the simulation progresses.
For example, we found some material still accreting onto the remnant at 100 days for  = 1.79 and a 3 M ⊙ MAMS model.

Mass of remnants
Table 2 lists the remnant masses for each model.As anticipated, TAMS model stars can survive higher  encounters compared with the ZAMS models, owing to their more centrally concentrated structure.Recall that for these stars, the same  corresponds to more distant encounters compared to ZAMS models because these stars have initially lower average densities.For  ≲ 1, stars exhibit varying degrees of mass losses, ranging from cases with no mass loss ( ≲ 0.5) to others with losses ≤ 10%.
Figure 5 shows the mass of the remnant as a function of the penetration factor, .The remnant mass is defined as described in Section 2.3.1.Additionally, we plot polynomial fits given by where , ,  are the parameters determined by fitting, as lines connecting our data points.Table 4 gives the fitted parameters as well as the root mean squared error in the fit.Ryu et al. (2020a) determined a functional form for mass of remnants in the Newtonian limit as  rem / ★ = 1 − ( p /R t ) −3 where R t is the physical tidal radius -defined as the maximum distance at which star would experience full disruption.Since  t does not account for the internal structure of the star, it can survive the disruption at distances  p <  t , aligning with our findings.Ryu et al. (2020b) provides R t / g values of 22.5 ± 1.2, 33.9 ± 2.0, and 52.1 ± 3.1 for their 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ MAMS models, where  g = 2  • / 2 .They also list the  t / g as 47.5, 79.8, and 123, respectively.Utilising these values, we calculated the  rem for a range of  values based on the  of our models.The results are shown as shaded regions in Figure 5. Notably, our 1 M ⊙ and 10 M ⊙ MAMS models fall within these shaded regions, but our 3 M ⊙ remnant masses surpass the expected range.We also observed that the remnant mass calculated by this functional formula is negative for higher  values, even in instances where remnants are still present for our simulations.Such an example is 3 M ⊙ MAMS disrupted for a  = 2.39.This inconsistency may stem from differences in the internal structures of our models.

Critical 𝛽 (𝛽 c )
We determined the critical , beyond which our stellar models would experience complete disruption ( c ), using the polynomial fits.The roots of these functions were determined, corresponding to  values where the remnant mass becomes zero.Values and uncertainties of  c are listed in Table 1.These uncertainties correspond to the coarseness of the simulations run.The fitting function determined a  c = 3.4 ± 0.5 for the 1 M ⊙ TAMS model.We found that a remnant can form ,  = 3.45 albeit of only 0.01 M ⊙ .
Law-Smith et al. ( 2020) and Ryu et al. (2020b) argued that  c is related to the ratio of critical density to the mean density.
Figure 6 shows the critical  as a function of the ratio of central density to the mean density of the initial stellar model.We conducted fits between  c and (  / ρ) and determined the following fit functions: = 0.39 × (  / ρ) 1/2.6 for 10 M ⊙ .(11c) than its initial radius for  > 1.For encounters with pericentre more distant than the tidal radius (e.g.,  = 0.54) we find only a slight increase in the radius of the remnant compared to the initial star, which is due to lower stripping of the star.

Properties of the remnants
Figure 8 shows density as a function of radius for the same set of models shown in Figure 7.The first column presents a comparison of the density distribution for the models mapped into Phantom without tidal disruption.The binning algorithm is observed to maintain the expected density profile.The central density of the remnant can be seen to decrease as the star passes closer to the SMBH (increasing ; left to right).For instance, the remnant from the  = 1.94 simulation of 1 M ⊙ exhibits a post-disruption central density of 2.02 g cm −3 , whereas for  = 0.54, the central density is 150.84 g cm −3 .In very weak encounters like  = 1.08, the remnant loses some material (∼ 3%), but the central density remains similar to the original model.We have also plotted the radially-binned density (Section 2.3.2) on the same plot for comparison, showing that our remnants are generally spherical except in the outermost layers where material continues to accrete onto the remnant.

Rotational properties
When a star is tidally disrupted by a SMBH, one might expect that the outer layers rotate faster than the core due to the stronger impact of tides.This, however, is not what happens in our simulations.Figure 9 shows the time evolution of the angular velocity (  ), computed in cylindrical coordinates in a cross-section slice in the orbital plane for our 3 M ⊙ MAMS model for a penetration factor  = 1.79.At 16 h before pericentre (top left panel), we see the formation of two distinct regions in the radially stretched star, with one region rotating in the opposite direction.These vortices transfer angular momentum from outer regions to the central region of the star, resulting in a remnant rotating faster in the central region (∼ radius of the initial star) while the fluffy outer layers rotating slower ( = 24 h).The time  at which these vortices disappear varies by model, being longer for closer encounters.Goicovic et al. (2019) found something similar but their vortices were persistent structures, and were still present at 42 h.
Figure 10 shows the angular velocity as a function of spherical radius.Initially, our models exhibited no rotation, as shown in the first column.We have also plotted the binned angular velocity (solid lines) and break-up velocity (dashed lines) given by for comparison. enc is the enclosed mass at radius .For disruptions where the star loses almost no mass (e.g.,  = 0.54 model for 1 M ⊙ MAMS; second panel from left in the top row), the remnant is rigidly rotating with a few outer zones having higher rotation.None of the regions rotating rigidly have rotation close to the break-up velocity.
For  ≳ 1 (rightmost panels), the differentially rotating layers are close to the break-up velocities.In summary, models that lose ≥ 40% of their initial mass result in slower-rotating remnants compared to models that lose less mass.
Figure 11 shows the formation of a circum-stellar envelope in the case of 3 M ⊙ MAMS,  = 1.79 simulation.The top panel shows what is identified as a remnant based on our energy condition in the face-on view (- plane), while the bottom panel shows the edge-on view (- plane) over 100 days.As time progresses, the outer layers of the remnant expand and result in the formation of an envelope that is bound to the remnant.
Figure 12 shows the evolution of angular velocity with time.We found that as time progresses, the region that was rotating rigidly starts to rotate slower, and the differentially rotating zones, rotate faster.This could be due to angular momentum transfer from the core to the outer region, and as the material is still falling back, it could also provide angular momentum to the outer layers.This could affect the evolution of the stars, but as we did not implement angular velocity in our mapped back models, we can not quantify the possible difference.

Escape velocity of remnants
Table 2 lists the initial and final escape velocities derived from our simulations.For the initial escape velocity, we followed the trajectory of a single point mass using a code that solves the geodesic equations (Liptai & Price 2019).We started these simulations at ∼ 10  t .We progressed these simulations until the geodesic reached a distance > 60  t for 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ models.We calculated the velocity by determining Newtonian kinetic and potential specific energy ( =  kinetic +  potential ) using  =

√
2.This approach remains valid at this stage as geodesics are sufficiently far away from the black hole, rendering the GR corrections negligible.
To determine the final escape velocity of the remnant, we initially calculated the centre of mass position and velocity relative to the SMBH at 5 days, 8 days, and 10 days since the start of the simulations.
As the encounter strength increases, the star's orbit undergoes more extreme changes.We notice that TAMS models slow down while the ZAMS and the MAMS models tend to speed up.For the same value of  the more evolved TAMS stars have a larger periapsis distance which decreases the effects of general relativity.Whether the postencounter remnant is faster or slower than the incoming star, however, may depend on the detail of the tidal interaction, which could be affected by the interaction of the dynamical tide with oscillation resonances excited within the star (Zahn 1975;Mardling & Aarseth 2001).
We also found that a few remnants formed by disrupting 10 M ⊙ TAMS model had negative energy.We calculated the semi-major axis and eccentricity of these orbits by using  = − ( • +  rem )/(2) and  = 1 −  p /.The results are given in Table 2 and Table 3.
We found that the velocity increases with time.Hence, there is some variability in the presented results.For example, for 3 M ⊙ ,  = 4.5 scenario, the remnant escape velocity is 384 km s −1 at 8 days, but at 20 days, the velocity is 425 km s −1 .The remnant mass and radius continue to change slightly with time, impacting the centre of mass calculation and the remnant's velocity.

Composition mixing during disruption events
Figure 13 shows the 14 N abundance perpendicular to the orbital plane for 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ MAMS models (top to bottom) for a range of penetration factors (across), at 4.8, 7.3, and 9 days postdisruption.We see that as the strength of the encounter increases, there is more mixing inside the star as the disruption can affect the deeper layers of the model.We note that in the disruptions where mass loss is ≲ 10% of the initial model's mass, there is no mixing inside the star.As the star approaches the SMBH at a closer distance, more material is mixed.Moreover, the streams have more material stripped from the central regions as the disruption encounter becomes more violent.
Figure 14 shows the 4 He and 14 N mass-fractions as a function of the enclosed mass of remnants formed by disrupting a 3 M ⊙ MAMS model with  = 1.49(top) and  = 2.39 (bottom).The green line illustrates the radially-binned mass fraction of the 3 M ⊙ before disruption, the blue line represents the remnant formed after disruption, and the red line shows the mass fraction obtained from a Kepler model of the same mass.Kepler models are evolved until the central 4 He mass fraction matches that of the remnant.We see that the  = 1.49,where the model lost ∼ 13% mass, the mass fraction of the remnant is slightly higher than the before disruption mass-fraction for mass enclosed > 0.7 M ⊙ due to mixing.The Kepler model of the same mass has overall lower mass fractions of 4 He and 14 N.
The remnant formed by disrupting this star,  = 2.39 results in a mass loss of ∼ 83%.There is more mixing in the remnant.We can see 4 He enrichment in the envelope.The Kepler model of the same mass has lower mass-fraction than the remnant. 14N content in Kepler model is 0.38 times lower than that of the remnant.

Evolution of remnants
Figure 15 shows the evolution of all the models following the disruption event after their respective Kelvin-Helmholtz time.The simulations were carried out until the core ignited core helium burning.Models with mass < 0.5 M ⊙ posed challenges for re-evolution in Kepler and therefore are not shown in the plot.
From ZAMS models (left panel of Figure 15), we see that the remnants formed by disruption evolve on a track that is nearly identical to an undisrupted star with the same mass as the remnant.For example, the disruption of the 10 M ⊙ ZAMS model with  = 1.48 and the 3 M ⊙ ZAMS model with  = 1.28 end up with masses of 2.58 M ⊙ and 2.66 M ⊙ , respectively.Their tracks (purple dashed line and blue dashed lines lying close to each other at log(/ ⊙ ) ∼ 2; see legend) are nearly identical (although the slightly lower mass model exhibits correspondingly lower luminosity and effective temperature, as expected).Similarly, the remnant of the 3 M ⊙ ZAMS model with  = 1.78 and the remant of the 1 M ⊙ ZAMS model with  = 1.28 have masses of 0.64 M ⊙ and 0.58 M ⊙ , respectively.They follow similar evolution tracks (blue and orange lines lying close to each other in the bottom right of the diagram).This is because the composition of the star is unchanged by the disruption process.
For MAMS models (middle panel of Figure 15), the situation is different.The remnants formed from disrupting a higher mass star have higher luminosity and effective temperature in comparison to undisrupted stars of the same mass.For example, the remnant of the disruption of a 10 M ⊙ model with  = 1.95 (lowest purple line in the middle panel), has higher luminosity than a 3 M ⊙ star (solid blue line in middle panel) despite having lower mass (2.85 M ⊙ ).Hence, in this case, we have produced an over-luminous star (see Section 4).Another example is disruption of a 3 M ⊙ MAMS model with  = 2.39 which results in a remannt with higher luminosity than the disruption of a 1 M ⊙ MAMS model with  = 1.29 despite having a lower mass.
For TAMS models (right-most panel of Figure 15), the situation is similar to the MAMS models.For example, the remnant of 10 M ⊙ TAMS model with  = 4.23 (last purple line lying in the region of blue lines) has higher luminosity than a 1 M ⊙ star which never encountered an SMBH.This remnant is also more luminous compared to a remnant of 1.26 M ⊙ formed by disrupting 3 M ⊙ TAMS model with  = 3 (dotted blue line after the last purple line).The disruption of a 3 M ⊙ TAMS model with  = 4.5 yields a remnant that has higher luminosity than the 1 M ⊙ model with  = 2.16 despite the former having a smaller mass (blue and orange lines in the right corner).
We also found that disrupted stars have longer lifetimes.Figure 16 shows the luminosity as a function of time for remnants of 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ MAMS models.The remnants can live longer on the main-sequence as the mass of the remnant decreases.
A relationship between the luminosity of the star and its mass has been noted for main-sequence stars (Kuiper 1938).Figure 17 shows this relationship for our remnants at the Kelvin-Helmholtz time.As we can see, the remnants do not follow this relationship, showing that these stars are unlike undisrupted stars.

Nitrogen enrichment
Figure 18 shows the 12 C and 14 N abundance at the surface of the remnant as a function of penetration factor at the Kelvin-Helmholtz time.ZAMS models have the least 14 N enrichment and some enrichment does take place for  > 1.The remnants of more evolved initial stars have the most amount of enrichment.The enrichment becomes significant for deeper encounters.This would lead to significant changes in the observable [C/N] ratio and provide a signature for detecting TDE survivor stars.

Comparison with non-disrupted stars of same mass
Figure 19 shows the H-R diagrams of 3 M ⊙ remnants at ZAMS, MAMS and TAMS (across).Solid lines correspond to models that were disrupted in Phantom and then mapped back into Kepler.The dashed lines correspond to the Kepler models of the same mass which never encountered an SMBH.We can see that the ZAMS models follow similar evolution, but for MAMS and TAMS, the disrupted models have higher luminosity and higher effective temperature.The difference in evolution becomes significant as we move towards lower mass remnants.

Comparison with stars that had their outer layers removed
We wanted to understand how different the remnants are from the original star with simply the outer layers removed.To achieve this, we used the mapped files of our remnants but used the interpolated composition from the original undisrupted star for the remnant mass (referred to as stripped star).
Figure 20 shows the H-R diagram of 10 M ⊙ models.The dashed lines represent the stripped models, while solid lines are the remnants.We show the evolution after 40,000 years.We also show the Kelvin-Helmholtz time of each model (starred points).We see that the more a star is disrupted, the lower its luminosity is compared with a stripped star.For example,  = 3.29 remnant has lower luminosity than the stripped star of the same mass.In the encounters where the star lost only ∼ 45% of its mass, the evolution of the remnant is similar to the stripped star.We found a similar behaviour for the MAMS models.
Fallback causes a H-rich surface layer, unlike just stripping the outer layers.Yet we retain the pre-disruption composition in the   = 2.39 (bottom).The green line illustrates the radially-binned mass fraction of the 3 M ⊙ before disruption, the blue line represents the remnant formed after disruption, and the red line shows the mass fraction obtained from a Kepler model of the same mass.Kepler models are evolved until the central 4 He mass fraction matches that of the remnant.Notably, the remnant retaining ∼ 87% of the initial model's mass experiences slight composition mixing.The remnant retaining ∼ 21% exhibits higher composition mixing.Additionally, the Kepler models display lower mass fractions at lower enclosed masses for both these models. 14N content in Kepler model is 0.38 times lower than that of the remnant formed with  = 2.39.core.Figure 21 shows mapped 4 He from Phantom to Kepler postdisruption as a function of enclosed mass for 10 M ⊙ TAMS models for a range of penetration factors.This shows that disrupted stars are unlike stripped stars for more deeply penetrating encounters.

DISCUSSION
Our study has several inherent limitations.Firstly, we did not incorporate rotation into the mapping within Kepler.Additionally, we did not disrupt stars present in the Nuclear Star Cluster of the Milky Way.Furthermore, the boundary distinguishing partial and full TDEs The evolution of remnants of similar mass formed from ZAMS models is not significantly different.However, MAMS and TAMS remnants exhibit higher luminosity and temperature for those formed from a higher-mass star.requires further exploration, necessitating an in-depth understanding of how the spin of the star and SMBH influence such encounters.

Properties of remnants
Coughlin & Nixon (2022) argued that the the tidal stripping of the star would take place for  ≳ 0.6 which is in agreement with our numerical results.Partial tidal disruption events manifest as violent occurrences, inducing oscillations in the remnant star close to the fundamental mode, accompanied by an increase in central density near the pericentre (Figure 3).Guillochon & Ramirez-Ruiz (2013) also identified post-pericentre oscillations in their tidal disruption simulations of polytropic models.The small increase in density before pericentre passage is not obvious in their plots since they use a log scale.Nixon & Coughlin (2022), utilising polytropic models, found that the central density of their star increased near the pericentre.They argued that the first spike is a result of vertical oscillation, while the second spike is caused by in-plane caustics.
TDEs lead to the stripping of stars, resulting in mass moss, as has been reported by several authors (Alexander 2005;Lodato et al. 2009;Guillochon & Ramirez-Ruiz 2013;Manukian et al. 2013;Goicovic et al. 2019;Ryu et al. 2020a,b,c).The amount of mass lost increases with the impact parameter () as shown in Figure 5. Sometimes a star might be completely destroyed but material can fall back onto itself to form a remnant as has also been found by (Guillochon & Ramirez-Ruiz 2013;Nixon & Coughlin 2022).Evolved stars, such as TAMS can survive disruption closer to the SMBH compared to ZAMS stars, consistent with the findings of Law-Smith et al. (2020).Remnants exhibit lower densities, temperatures and larger radii after the encounter (Figures 7 and 8).In contrast to the prediction by Alexander & Livio (2001) that the radius of the remnants would be smaller, our results suggest otherwise.Additionally, they had predicted that disrupted stars would experience chemical mixing, which we also observed in our models (Figures 13 and 14).
We observed the formation of vortices that can redistribute angular momentum and result in the central regions rotating faster (Figure 9).This observation aligns with Goicovic et al. (2019), who disrupted a 1 M ⊙ ZAMS MESA model (Paxton et al. 2011), although their vortices persisted for the duration of their simulations (44 hours).
Our remnants exhibit rigid body rotation in the central regions (close to the radius of the initial model), and differential rotation in outer layers (Figure 10).This is unlike Goicovic et al. (2019) who found that their core did not exhibit rigid rotation.The discrepancy might be attributed to their models having an initial rigid rotation profile, whereas we assumed initially non-rotating stars.Sacchi & Lodato (2019) disrupted a 1 M ⊙ ,  = 5 /3 polytrope with  = 1.Unlike our approach, they provided the initial model with initial rigid body rotation.In their simulations, they found that a disk can form around 3 days after disruption, tightly bound to the core and exhibiting Keplerian rotation.This is similar to our observation, although we do not see the formation of a disk in our models.We speculate that the discrepancy might be due to their model having an initial rotation profile close to differential rotation.Further work is required to resolve this difference.
Nixon et al. ( 2021) also found a disk in for  = 1.6, 1 M ⊙ ZAMS model, although they did not specify if it was bound to the star.The cause of the discrepancy with our results remains uncertain.They argued that the disk formed around the core is spun off from it.In our simulations, the differentially rotating region belongs to the outer layers of the initial star that were stripped away.As the simulation progresses, material from the outer layers of the rigidly rotating core can become part of these differentially rotating layers; moreover, the streams can fall back onto the star.Ryu et al. (2020c) also calculated the angular velocity, but unlike our models, their models rotate faster in the outer layers.They found that for the same fractional mass loss, the models that have higher initial mass get closest to break-up velocity.We do not observe such a trend.For example,  = 1.29 for a 1 M ⊙ MAMS, and  = 2.39 for a 3 M ⊙ MAMS, form similar mass remnants, but the remnant formed from the latter model is not closer to break-up velocity than the former models' remnant.
By evolving some of our models to long times, we found that the size of the remnant increases (Figure 11 and Figure 12).Moreover, the angular momentum is redistributed with time, with differentially rotating layers rotating faster while the core rotates slower.
Further investigation is necessary to determine the cause of discrepancies in the rotation profiles within different simulation codes.

Critical 𝛽 for total disruption
We calculated the critical  ( c ) values as 1.56, 2.04, 3.40 for 1 M ⊙ ZAMS, MAMS and TAMS models as described in Section 3.5.For 3 M ⊙ model,  c was determined as 1.91, 2.59, and 4.8 for ZAMS, MAMS and TAMS stages.In the case of 10 M ⊙ , we estimated it as 1.59, 2.17, and 4.49 ZAMS, MAMS and TAMS stages, respectively.Guillochon & Ramirez-Ruiz (2013) were the first to explore  c for polytropes with  values of 4 /3 and 5 /3 as 1.85 and 0.9, respectively.As polytropic models have lower central densities, the star would not be able to survive encounters closer to the SMBH, compared with a star of higher central density.
Law-Smith et al. ( 2020) utilised MESA stellar models and disrupted them in Newtonian physics.For their 1 M ⊙ ZAMS model, they determined a  c of 1.8 ± 0.1.Their model had a radius of 0.9 R ⊙ and  c / ρ of 42, which is similar to our model.Hence, the difference in  c might stem from using general relativistic effects in our simulations.Their TAMS model was obtained at 8.4 Gyr and  We also plot the evolution of stars which would have only been stripped by the SMBH.All models are plotted 40, 000 years since the start of Kepler evolution.We show the Kelvin-Helmholtz time as star points for each model.We see that stripped stars are more luminous compared with our remnants.
had a radius about 5% larger and  c / ρ about 75% higher than our model.They determined a  c of 7 for this model, about twice of what we found.This can be attributed to the difference in the radius and density of the stellar models.Their 3 M ⊙ ZAMS and TAMS models had radius of 1.89 R ⊙ and 3.32 R ⊙ with  c / ρ of 73 and 1198.Their ZAMS  c was about 6% higher than ours.This could be either a result of a difference in  c / ρ or general relativistic effects. c / ρ for our TAMS model was 24% lower.This could be a possible explanation for the difference in  c .Goicovic et al. (2019) used 1 M ⊙ ZAMS model from MESA, and determined that the full disruption takes place after  ≥ 2. Hence, there is a slight difference in the literature results.Ryu et al. (2020c) used real stellar models that had spent about half of their lives on main-sequence and disrupted them in a fully general relativistic framework.There is a small discrepancy between their  c values and ours for the 1 M ⊙ and 10 M ⊙ MAMS models.This could be due to us using hyperbolic orbits for the simulations as for deep encounters, parabolic orbits lose more mass than the hyperbolic cases as discussed in Section D. As Ryu et al. (2020c) took their models from MESA while we utilised Kepler, this could also be a reason for the difference.The small discrepancies could also be a result of different simulation methods (Mainetti et al. 2017).
We also explored the correlation between  c and ( c / ρ).The results are given in Eq. 11.Law-Smith et al. (2020) performed the fit based on the distinction in density ratio, proposing a function of the form  c = 0.5 × ( c / ρ) 1 /3 for ( c / ρ) ≲ 500 and  c = 0.39 × ( c / ρ) 1 /2.3 for ( c / ρ) ≳ 500.Our 1 M ⊙ TAMS and 10 M ⊙ TAMS models, however, have  c / ρ of 569.76 and 563.72, respectively.As models below a mass of 1.3 M ⊙ have convective cores while those above this mass have radiative cores, fitting the same function to  c would overlook the evolutionary structure of these stars.
Law-Smith et al. ( 2020) used Newtonian physics in their models and observed remnant formation at higher  values than the  c that we determined.Additionally, the  c / ρ values of their models differ from ours, owing to slight differences in stellar evolution codes.Ryu et al. (2020b) derived a semi-analytical formula of the form  c = 2.14 × ( c / ρ) 1 /3 .Coughlin & Nixon (2022) determined a a functional form of  c = 0.62 × (  / ρ) 1 /3 and found good agreement with previous numerical simulations.But Nixon & Coughlin (2022) had found that the star can survive the disruption encounter for  >  c for their  = 5 /3 polytrope,  = 16, while previous studies had determined a value of ∼ 0.9 (Guillochon & Ramirez-Ruiz 2013;Mainetti et al. 2017).The cause of this discrepancy is unclear.

Evolution of remnants
Hydrogen burning in the low-mass stars is dominated by the PP chain, however, for more massive stars the CNO cycle starts to dominate hydrogen burning.At first, 12 C is converted to 14 N very swiftly at the beginning of burning.This is the CN cycle.There is also a branch that converts 16 O to 14 N to participate in the CN cycle, however, this conversion by the ON branch requires higher temperatures.It is usually only encountered toward the end of the evolution of a 1 solar mass star and still proceeds much slower than the initial conversion of 12 C even in higher-mass stars.Usually, not much of that occurs at the early evolution stage we refer to as ZAMS but will be present increasingly more with stellar age and mass, so, e.g., more present in the MAMS 10 M ⊙ star than in the MAMS 3 M ⊙ star.Hydrogen burning by the CNO cycle -dominated by the CN branch -hence is enhanced when 16 O has been converted to 14 N, which would not occur in the early evolution of low-mass stars.Hence, if one made a low-mass star from the core of the high-mass star, the resulting larger 14 N can affect the evolution as a larger CN isotope abundance can burn hydrogen faster, i.e., a lower core temperature is needed to produce the same luminosity.
The other key factor is the enrichment of helium.For chemically homogeneous stars, the luminosity is proportional to the fourth power of the mean molecular weight,  (Kippenhahn & Weigert 1990), i.e.,  4 .The change in mean molecular weight during hydrogen burning is dominated by the change in helium mass fraction -the change in the CNO isotopes is negligible in comparison.Hence the helium-rich core of a stripped evolved star is more luminous than a non-stripped star of the same mass and same core helium content, as it will have more helium in the outer layers and hence a higher average mean molecular weight.The higher helium abundance in the outer layers also leads to a lower electron scattering opacity, however, our H-R diagram shows that the stars do become redder, e.g., well-seen for the 10 M ⊙ models.Alexander & Livio 2001;Alexander 2005 had predicted that a partial tidal disruption would lead to remnants of main-sequence stars exhibiting higher luminosities for a given mass due to higher mean molecular weight and a larger convective region.Our findings support this proposition, as evidenced by a comparison of our remnants with stars of the same mass, revealing higher luminosities and temperatures (Figure 19).However the remnants are less luminous compared to stars where only the outer layers are removed.This is because the remnants retain the pre-disruption composition in the core and have a H-rich surface layer (Figure 20 and Figure 21).
Figure 15 shows that as the remnant mass decreases, it traverses the H-R diagram, towards lower effective temperature and luminosity.This movement is not smooth for TAMS models.Moreover, the observed correlation between mass loss and extended main-sequence lifetimes suggests that stars experiencing greater mass loss can reside on the main-sequence for a longer duration -a true 'elixir of life' (Figure 16).
We also found that disruption can result in depletion of 12 C and enrichment of 14 N at Kelvin-Helmholtz time as was predicted by Alexander & Livio (2001) (Figure 18).

Bound or unbound?
Stars with velocities greater than the local escape velocity of the galaxy can be either classified as hyper-velocity stars (HVSs) or hyper-runnaway stars (HRSs), where the former have a galactic centre origin (Hills 1988;Brown 2015;Marchetti et al. 2018).An ejection rate of 10 −2 -10 −4 yr −1 has been estimated for HVSs (Brown 2015;Evans et al. 2022).Detection of HVSs can help us to understand the galactic centre environment.
Since the first detection by Brown et al. (2005), 20 HVSs have been confirmed with velocities of 300-1,700 km s −1 (Brown et al. 2014;Brown 2015;Koposov et al. 2020), and more than 500 candidates have been identified (Boubert et al. 2018;Luna et al. 2024).Most of these stars are B-type stars with masses ranging from 2.5-4 M ⊙ (Luna et al. 2024).Only Koposov et al. ( 2020)'s A-type HVS with a velocity of ∼ 1,700 km s −1 has an orbit that points towards a galactic centre origin and a result of Hill's mechanism (Hills 1988).The current astrometry is not precise enough to trace other HVSs to a galactic centre origin (Evans et al. 2021).
Our remnants formed by disrupting stars on hyperbolic orbits result in velocities at infinity of 200-2,000 km s −1 .This is different from the result of Ryu et al. (2020c) who had found that all their remnants formed from disrupting 1 M ⊙ , 3 M ⊙ and 10 M ⊙ MAMS stars were bound to the SMBH on highly eccentric orbits.They used parabolic orbits while we used hyperbolic orbits for our paper.If we compare our results from the zero energy orbits discussed in Section D, we note that only  = 1.49disruption results in an unbound remnant with a velocity ∼ 600 km s −1 .All the subsequently bound models follow highly eccentric orbits around the SMBH post-disruption.Some of our remnants formed from disrupting 10 M ⊙ models were also bound to the SMBH on highly eccentric orbits.These could have repeated disruptions with the SMBH (Antonini et al. 2011;Ryu et al. 2020c).Manukian et al. 2013;Gafton et al. 2015 had found that single stars on parabolic orbits before partial disruption could have velocities of a few hundred km s −1 after interaction with a SMBH, due to asymmetry in the mass loss.This can explain our remnant formed by disrupting a star on a zero-energy orbit with  = 1.49.
Disruptions on hyperbolic orbits could result in HVSs.Though Ryu et al. (2020c) had argued that unbound remnants can experience a change in their angular momentum from gravitational encounters which can alter their trajectories.With a logarithmic potential, post sphere of influence, the unbound remnants could reach a turning point, and end up in the galactic centre but with a small possibility of another disruption with the SMBH.

Possible observational candidates
Due to difficulty in observing stars in the galactic centre, there are no present partial TDE star candidates.Here, we propose a few stellar objects that might be a result of an encounter with a SMBH.About six G-objects have been detected within 0.04 pc of our galactic centre (Gillessen et al. 2012;Phifer et al. 2013;Witzel et al. 2017;Ciurlo et al. 2020).Ciurlo et al. (2020) argued that Br is the main identified feature of these objects which results from ionisation of dust and is independent of the mass of the possible central object which could be a star.There is a lack of  ′ -band emissions but it does not rule out the possibility of a remnant embedded within these ionised envelopes.A few possible explanations have been proposed for explaining the formation of these objects, such as gas or dust clouds (Pfuhl et al. 2015), a star with stellar wind (Scoville & Burkert 2013), or a star having an evaporating protoplanetary disk around a young star (Murray-Clay & Loeb 2012), or a product of stellar binary collision (Prodan et al. 2015).As the dynamics are similar to a stellar object, we propose that partial tidal disruption of a star could have resulted in such objects.As our simulations evolve, we see that the star increases in size (object identified as being a star based on Section 2.3.1),forming a bound envelope around the dense central region.As we looked at some of our models only about 100 days, the exact comparison between the sizes of G-objects and our remnants is not possible.We also did not consider cooling in our models which could have an impact on the structure of these outer layers.Furthermore, our stars had no rotation in the beginning.All these possibilities need to be further explored.Alexander & Livio (2001) had identified IRS 7, an M supergiant in the galactic centre as a candidate for a partial TDE.It was found to have abundances that were consistent with the dredge-up of CNO cycle products but required deep mixing in excess compared to standard models of red giant stars (Carr et al. 2000).But recently, Guerço et al. (2022) showed that the chemical abundances are similar to rotating stellar models.We found tidal disruption can cause mixing in the star.Moreover, the remnants can rotate faster compared to the original star (see Section 3.7).It would be interesting to compare the effects of rotation (Heger et al. 2000;Meynet & Maeder 2000) induced by disruption on the evolution of a disrupted star, and see if they can explain IRS 7 or not.Schiavon et al. (2017) found N-enriched stars in the inner region of the Milky Way galaxy.Some of these stars might be products of tidal disruption from the SMBH.As shown in Figure 18, remnants of partial TDEs can have Nitrogen enrichment with Carbon depletion at Kelvin-Helmholtz time after the disruption.But the detected Nenriched stars are uniformly spread out and move in a way that cannot be distinguished from other stars within the same volume (Schiavon et al. 2017), kinematical analysis is required to understand if this can be a formation mechanism for such stars.

CONCLUSION
In this paper, we obtained 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ models at zeroage, middle-age and terminal-age main-sequence for a range of  values using the 1D stellar evolution code Kepler.Then we disrupted these models in 3D general relativistic hydrodynamics using Phantom, and finally mapped these back into Kepler to understand the post Kelvin-Helmholtz evolution.
Our key results are (i) 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ stars can all survive deep encounters ( ≳ 1), and result in formation of remnants.
(ii) These remnants rotate rigidly in the core and exhibit differential rotation outside.This differentially rotating region forms a large circumstellar envelope.
(iii) Around pericentre, there is an increase in central density, but as stripped material moves away from the star, the density decreases.In the end, this leaves behind remnants with lower densities and temperatures.
(iv) Remnants can live on the main-sequence for much longer than the original star they were disrupted from.Thus partial disruption of a star supplies the "elixir of life"!(v) Remnants formed from different mass ZAMS stars evolve similarly to undisrupted stars of a lower mass.But for MAMS and TAMS, the remnants of higher mass stars have higher luminosities and temperatures than an undisrupted star of the same mass.
(vi) We find that partially disrupted stars can show surface enrichment of 14 N, with more enrichment occurring for deeper encounters.
(vii) Remnants are more luminous compared with stars of the same mass that were never disrupted by a SMBH.
(viii) Remnants are less luminous compared with stars where only the outer layers were removed.This is because the remnants retain a H-rich outer layer.but the model does not look spherical.Hence, we used the snapshot at 60 days for analysis.

APPENDIX C: CONTROL CASE: KEPLER MAPPING BACK
To test that our Kepler mapping is robust, we compared the evolution of 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ ZAMS models that were mapped into Phantom but not disrupted.Figure C1 shows the evolution of models in Kepler before they were mapped into Phantom (orange line), and models that were mapped back into Kepler from Phantom without disruption (blue line).We can see that the 1 M ⊙ mapped back model leaves the main-sequence at about at 99.8% of the original model's time.When the models reaches MAMS, the mapped back model has an luminosity of 0.97 L ⊙ which is the same as the luminosity of the original model, but there is a difference in effective temperature.The mapped back model has an effective temperature of 5,773 K while the original model has an effective temperature of 5,765 K.This corresponds to an error of 0.13%.Hence, there are a few discrepancies.We found small differences in the mapping of the 3 M ⊙ model.For 10 M ⊙ , mapped back model, we see oscillations but this effect is not significant enough for us to discard this mapping.Overall, our mapped-back models are in good agreement with the original model's evolution, albeit small differences do exist.
Figure C2 shows the luminosity and effective temperature as a function of time of 3 M ⊙ TAMS with  = 2.5.The model was mapped back at 8 days, 50 days, 100 days.We see that all models follow similar trends in luminosity and effective temperature as a function of time.Hence, the time at which models are mapped back does not have a significant effect on stellar evolution.

APPENDIX D: IMPLEMENTING 𝐸 = 0 ORBITS IN PHANTOM
Here, we explain the appropriate way for implementing zero energy orbits in Phantom.We used the geodesic code (Liptai & Price 2019) to evolve 100 test particles for a range of Newtonian pericentre distances (47.6-196.8R ⊙ ).but computed from a large distance of 2.4 × 10 5 R ⊙ where the Newtonian approximation is valid.We then calculated the position and velocity of the particle at our starting distance of ten times the tidal radius for 1 M ⊙ ZAMS model, once all the simulations were finished.This was achieved by interpolating to find the starting position and velocity closest to our desired initial separation in the geodesic trajectory.
The resulting position and velocity were used for setting up orbit   We analysed snapshots at 8 days, 50 days, and 100 days from the start of the simulation shown as blue, orange and green lines.The results indicate good agreement in the star's evolution, suggesting that the time of analysis for the remnant snapshot does not significantly impact the outcome.
Table D2.This table lists the semi-major axis and eccentricity of the bound models presented in Table D1.Bound remnants were formed by disrupting the 1 M ⊙ ZAMS model on  = 0 orbits.Column 1 gives the impact parameter .Column 2 lists the semi-major axis.Column 3 gives the eccentricity (), specifically log 10 (1 − ).We assumed Keplerian orbits for our analysis.The top panel shows a hyperbolic orbit and the bottom panel shows a zero energy orbit.Both simulations were evolved 8 days and each snapshot is at a 3 day interval.We see that the zero energy orbit is slightly more perturbed than the hyperbolic orbit.Moreover, the hyperbolic orbit remnant retains 2.5 times higher mass than the zero energy orbit.Table D1 lists the simulation results for both hyperbolic and zero energy orbits run for the same  for comparison.Column 1 gives the orbit type.Column 2 gives the  value.Column 3 lists the time of comparison of models for the same .Columns 4 and 5 list the mass of the remnant and final escape velocity.We see that the models have similar mass at lower  values, but as the star gets closer to the SMBH, the remnant mass is lower for zero energy orbit compared with hyperbolic orbit for the same .Moreover, all remnants of zero energy orbits are bound except for the remnant of  = 1.49which has a velocity of 634 km s −1 .This velocity might be the 'kick' caused by the streams on the remnant as argued by Manukian et al. (2013).showing the column density perpendicular to the orbital plane.All simulations were run for 8 days and each snapshot is at a 3 hour interval.We note that orbits are similar, though zero energy orbit is slightly more perturbed, and has more mass being accreted onto the SMBH.

Figure 1 .
Figure 1.Hertzsprung-Russell (H-R) diagram of the 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ stars we evolved in Kepler.We also show the three different stages (ZAMS, MAMS, and TAMS) of the stars for which the models were obtained.

Figure 2 .Figure 3 .
Figure 2. Snapshots of 3 M ⊙ MAMS model disrupted at different impact parameter (; top to bottom), showing the column density perpendicular to the orbital plane.The grey circle corresponds to the tidal radius.All simulations were evolved for 8 days and each snapshot was at a 3 h interval.The star loses more mass as it approaches the SMBH at a closer distance.Streams are formed around the remnant which fall back onto the remnant.

Figure 5 .
Figure5.Mass of remnant as a function of penetration factor, .The lines connecting the points represent the fitting functions that were determined for each model.We also plot the remnant masses based on the Newtonian functional form given byRyu et al. (2020a) as shaded regions for their MAMS models.More evolved stars (MAMS, TAMS) can survive the encounter closer to the SMBH compared to the ZAMS because as stars evolve, the central density and radius increase.

Figure 7 Figure 6 .
Figure7shows the column density perpendicular to the orbital plane showing only the material identified as the stellar remnants, as per Section 2.3.1, for 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ MAMS models (top to bottom) for different initial  values (left to right).The models are shown at 4.8 days, 7.3 days, and 9 days post-pericentre (i.e., 5 days, 8 days, and 10 days since the start of the simulation).The first column shows the initial undisrupted star.The remnants exhibit dense cores (∼ radius of the initial star) with fluffy outer envelopes extending to ∼ 10-100 R ⊙ for the most penetrating encounters, a result of material stripping leading to expanded radii.The radius of the remnant also increases with increasing .The remnant's radius increases with higher , e.g., the 1 M ⊙ model's remnant is nearly 10 times larger

Figure 7 .Figure 8 .Figure 9 .Figure 10 .
Figure 7. Remnants of 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ MAMS models at 4.8 days, 7.3 days, and 9 days post-pericentre, respectively (top to bottom).Except, the 10 M ⊙ ,  = 2.15 is analysed at 59 days post-pericentre.We show column density perpendicular to the orbital plane, showing only material identified as belonging to the remnant (i.e.bound and not belonging to the streams).Each panel is 200 R ⊙ × 200 R ⊙ .Each row shows the models at different penetration factors, with the first column showing the models before the disruption.The radius of the remnants increases as the star approaches closer to the SMBH.The remnants bear a striking similarity to G-objects in the galactic centre (Ciurlo et al. 2020).

Figure 12 .
Figure 11. 3 M ⊙ ,  = 1.79 model remnant at different times (across).Each panel is 300 R ⊙ × 300 R ⊙ .We show the face-on (top) and edge-on (bottom) view.As time increases, the size of the remnant increases too.There is the formation of differentially rotating circum-stellar envelope which is bound to the rigidly rotating central region.

Figure 13 .
Figure 13. 14N composition cross section slices in the orbital plane of 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ MAMS models at 4.8 days, 7.3 days, and 9 days since disruption, respectively (top to bottom).Each row shows the models at different penetration factors, with the first column showing the models before the disruption.Stronger disruption events result in more material being removed from the centre of the star, which eventually ends up in streams.Some material in the streams falls back onto the remnant.

Figure 15 .
Figure15.The leftmost, central and rightmost panels show the H-R diagrams for ZAMS, MAMS and TAMS of remnants formed by disrupting 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ at a range of  values after the Kelvin-Helmholtz time.The evolution of remnants of similar mass formed from ZAMS models is not significantly different.However, MAMS and TAMS remnants exhibit higher luminosity and temperature for those formed from a higher-mass star.

Figure 16 .
Figure16.Luminosity as a function of time for MAMS models.The leftmost, central and rightmost panels correspond to remnants of 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ stars respectively.The models that lose more mass tend to have longer main-sequence lifetimes.

Figure 18 .
Figure18. 12C abundance (top) and 14 N abundance (bottom) at the surface of the remnant as function of  at a Kelvin-Helmholtz time after the TDE.With increasing , 14 N is enriched, whereas 12 C is depleted, leading to significant changes in the observable [C/N] ratio.

Figure 19 .
Figure 19.H-R diagram of 3 M ⊙ remnants for ZAMS, MAMS and TAMS stages shown in leftmost, central and rightmost panels respectively.Solid lines correspond to the evolution of remnants after the Kelvin-Helmholtz time.Dashed lines correspond to post hydrogen ignition evolution of models evolved in Kepler having the same mass as the remnant models.ZAMS remnants have similar evolution to the Kepler models.While MAMS and TAMS remnants have higher luminosity and effective temperature compared with Kepler model of the same mass.

Figure 20 .
Figure 20.H-R diagram of 10 M ⊙ TAMS remnants for a range of  values.We also plot the evolution of stars which would have only been stripped by the SMBH.All models are plotted 40, 000 years since the start of Kepler evolution.We show the Kelvin-Helmholtz time as star points for each model.We see that stripped stars are more luminous compared with our remnants.

Figure 21 .
Figure 21.Mapped 4 He from Phantom to Kepler post-disruption as function of enclosed mass for 10 M ⊙ TAMS models for a range of penetration factors.Pre-disruption composition is retained at the centre of the remnant.

Figure A1 .Figure A2 .
Figure A1.Temperature as function of radius of 3 M ⊙ MAMS disrupted with  = 1.5 at 8 days.We performed a cut in temperature on particles that were found to be bound to the remnant.The red line shows the temperature determined based the 3  condition.Then we refine this cut resulting in the final cut used for analysis.

Figure B1 .
Figure B1.Maximum density as function of time for 10 M ⊙ MAMS,  = 2.15 model.The central density decreases until ∼ 15 days and then increases, becoming constant at ∼ 30 days.

Figure C1 .
Figure C1.H-R diagrams of the original model compared with the mapped back model for 1 M ⊙ , 3 M ⊙ , and 10 M ⊙ (left to right) ZAMS models.The Kepler model is shown in orange and the mapped back model is shown in blue.We see that our mapped-back models are in good agreement with the original models.
Figure C2.H-R diagram of 3 M ⊙ TAMS model disrupted with  = 2.5.We analysed snapshots at 8 days, 50 days, and 100 days from the start of the simulation shown as blue, orange and green lines.The results indicate good agreement in the star's evolution, suggesting that the time of analysis for the remnant snapshot does not significantly impact the outcome.
in Phantom.FigureD1shows the snapshots of 1 M ⊙ ,  = 1.49, showing column density perpendicular to the orbital plane.

Figure D1 .
Figure D1.Snapshots of 1 M ⊙ ZAMS disrupted,  = 1.49.The top panel shows the hyperbolic orbit and the bottom panel corresponds to the zero energy orbit,showing the column density perpendicular to the orbital plane.All simulations were run for 8 days and each snapshot is at a 3 hour interval.We note that orbits are similar, though zero energy orbit is slightly more perturbed, and has more mass being accreted onto the SMBH.

Table 1 .
Properties of the nine Kepler models used in this paper.Columns 1 and 2 give initial mass and evolutionary stages, respectively.ZAMS, MAMS, and TAMS are zero-age, middle-age, and terminal-age main-sequence.Column 3 gives the age at which the tidal disruption occurs.Columns 4 and 5 list the central density and radius of the model respectively.Column 6 gives the ratio of the central density to the mean density.The last column lists the critical , where no remnant would be formed calculated by Section 3.5.The uncertainties are based on the coarseness of the simulations run.

Table 2 .
Properties of Phantom simulations.Column 1 lists the model evolution stage; Column 2, the impact parameter, ; Column 3, the mass of the remnant; Column 4 and 5, the star's initial and the remnant's final (error ±100 km/s) escape velocities; Column 6, the Newtonian pericentre,  p =  t /; and Column 7, the Kelvin-Helmholtz time of the remnant.

Table 3 .
Semi-major axis and eccentricity of the bound models presented in Table2.All bound remnants were formed by disrupting the 10 M ⊙ TAMS model.Column 1 gives the impact parameter , Column 2, the semi-major axis, and Column 3, the eccentricity, .These Keplerian orbit elements assume Newtonian dynamics only.

Table 4 .
Column 1 gives the mass of the stellar model.Column 2 lists the stage of evolution.Columns 3, 4, and 5 list the parameters of the polynomial fit.Column 6 gives the root mean squared error.

Table D1 .
This table lists the simulation properties of hyperbolic and zero energy orbits for 1 M ⊙ ZAMS models disrupted with different  values.The first column gives the orbit type that was set in Phantom.Column 2 lists the penetration factor.Column 3 lists the time of comparison between different orbits.Column 4 lists the mass of the remnant.Column 5 gives the estimated final velocity of the remnant at infinity based on its total energy.