On the Diversity of Fallback Rates from Tidal Disruption Events with Accurate Stellar Structure

The tidal disruption of stars by supermassive black holes (SMBHs) can be used to probe the SMBH mass function, the properties of individual stars, and stellar dynamics in galactic nuclei. Upcoming missions will detect thousands of tidal disruption events (TDEs), and accurate theoretical modeling is required to interpret the data with precision. Here we analyze the influence of more realistic stellar structure on the outcome of TDEs; in particular, we compare the fallback rates—being the rate at which tidally disrupted debris returns to the black hole—from progenitors generated with the stellar evolution code mesa to and γ = 5/3 polytropes. We find that mesa-generated density profiles yield qualitatively different fallback rates as compared to polytropic approximations, and that only the fallback curves from low-mass (1 M⊙ or less), zero-age main-sequence stars are well fit by either a or 5/3 polytrope. Stellar age has a strong affect on the shape of the fallback curve, and can produce characteristic timescales (e.g., the time to the peak of the fallback rate) that greatly differ from the polytropic values. We use these differences to assess the degree to which the inferred black hole mass from the observed light curve can deviate from the true value, and find that the discrepancy can be at the order of magnitude level. Accurate stellar structure also leads to a substantial variation in the critical impact parameter at which the star is fully disrupted, and can increase the susceptibility of the debris stream to fragmentation under its own self-gravity. These results suggest that detailed modeling is required to accurately interpret observed light curves of TDEs.


Introduction
The tidal destruction of a star by a supermassive black hole (SMBH) generates a stream of stellar debris that, over timescales of months to years, feeds the SMBH and generates a luminous, observable signature (Hills 1975;Lacy et al. 1982;Rees 1988). These tidal disruption events (TDEs) therefore offer one of the few means to directly probe the inner regions of otherwise-quiescent galaxies, and dozens have now been observed (e.g., Gezari et al. 2012;Arcavi et al. 2014;Chornock et al. 2014;Blagorodnova et al. 2017;Hung et al. 2017;van Velzen et al. 2019; see Komossa 2015 for a review). However, our ability to confidently use the observed flares from TDEs to study, for example, black hole (BH) demographics depends critically on our physical understanding of the disruption process, and in particular the way in which the properties of the progenitor star translate to a corresponding feeding rate of the SMBH.
Along these lines, Lodato et al. (2009) simulated the full disruption of polytropes, which offer simple, yet physical descriptions of the density profiles of stellar interiors (e.g., Hansen et al. 2004), and also provided a nearly analytical means of calculating the fallback rate from a star with a given density profile (using the impulse, or "frozen-in," approximation; see also Stone et al. 2013 and Section 3 below). They found that, while the fallback rate always approached the expected, t −5/3 scaling at late times, the peak value of the fallback and the time to peak depended on the stellar structure. Guillochon & Ramirez-Ruiz (2013) expanded this work by also studying the effect of varying the point of closest approach of the stellar center of mass to the BH, and found that denser stars more frequently leave bound "cores" that either resist the tidal shear altogether throughout pericenter passage, or reform following the full disruption of the star. The existence of these cores then modifies not only the early-time fallback, but also causes the latetime fallback to deviate from t −5/3 , and these features are a direct result of the stellar structure (in combination with the variation in the stellar pericenter).
Since then, a number of other authors have analyzed the tidal disruption of polytropes, with the aim of assessing one or another aspect of the tidal disruption process (e.g. Hayasaki et al. 2013Hayasaki et al. , 2016Coughlin & Nixon 2015;Shiokawa et al. 2015;Bonnerot et al. 2016;Coughlin et al. 2016Coughlin et al. , 2017Bonnerot et al. 2017;Guillochon & McCourt 2017;Mainetti et al. 2017;Wu et al. 2018;Golightly et al. 2019). However, only relatively few authors have analyzed the disruption of a star with a density profile other than a polytrope. We are aware of the following studies: MacLeod et al. (2012), who investigated the disruption of giant stars by particularly massive SMBHs; Law-Smith et al. (2017), who simulated the disruption of white dwarfs with extended, hydrogen envelopes; Gallegos-Garcia et al. (2018), who analytically calculated the fallback of metal-rich material from the cores of evolved stars; and Goicovic et al. (2019), who analyzed the extent to which a more realistic stellar structure affects the stellar "disruptability." Among the questions that remain concerning stellar structure is the degree to which more realistic stellar density profiles (i.e., those calculated with a stellar evolution code that accounts for more realistic opacities, metallicity gradients, and energy transport) affects the fallback rate compared to polytropes. More realistic stellar profiles could impose additional variability on the fallback rate, or conceivably alter characteristic timescales associated with the fallback (e.g., the time to the peak of the fallback curve). Such timescales are used to place constraints on BH properties (e.g., Guillochon & Ramirez-Ruiz 2013;Mockler et al. 2019), and hence it is necessary to understand the effects that stellar structure can have in modifying them.
In this Letter we analyze the disruption of stars with stellar structure computed with the code MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018, primarily to understand the influence that such structure can have on the fallback of debris to the SMBH. In Section 2, we first describe and present the stellar profiles calculated with MESA, and in Section 3 we use the impulse approximation-as described in Lodato et al. (2009)-to calculate the fallback rate onto the BH from those stars; we show that, compared to polytropes that are matched to the same stellar mass and radius of a given MESA model, there are certain combinations of initial stellar mass and age that yield notably different fallback curves. In Section 4 we present the results of numerical simulations of disruptions of the MESA models, and compare those results to disruptions of polytropes (again, matched to the same stellar mass and radius). We discuss the implications of our simulations for estimating the BH mass from observed light curves in Section 5, and show that polytropic approximations can lead to mass-estimate discrepancies at the order of magnitude level. We summarize and conclude in Section 6.

Stellar Profiles
Using the stellar evolution code MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018, we evolved a 0.3 M e , 1 M e , and 3 M e zero-age main-sequence (ZAMS) star to the end of the main sequence (MS), being the time at which the hydrogen mass fraction in the core dropped below 0.1%. For each star, we used all of the default values for the standard inputs (e.g., each pre-MS model adopted Solar metallicity, the stars were all non-rotating, there was no mass loss in the form of winds) within MESA, version 10398.
We took snapshots of the density of each star at the ZAMS, the terminal-age main sequence (TAMS), and at one time in between ZAMS and TAMS when the hydrogen mass fraction in the core just fell below 0.2; we denote the latter by a "middle-age main sequence" star (MAMS). 4 Figure 1 shows the density profile of each of these stars, and demonstrates, perhaps not surprisingly, that stellar evolution produces vastly different density structures over the lifetime of a given star.
In the following two sections, we describe two different approaches to modeling the tidal disruption of these stars by an SMBH.

The Impulse Approximation
A useful methodology for analyzing the fallback of debris from a tidal disruption event is the impulse approximation, which posits that the tidal field of the BH acts impulsively as the center of mass of the star reaches the tidal radius. Therefore, prior to reaching the tidal radius the star retains perfect hydrostatic balance, and thereafter the star is "destroyed," meaning that each gas parcel follows its own ballistic orbit in the gravitational field of the BH. Within this approximation and to lowest order in the tidal potential, the binding energy of a given fluid parcel is only a function of the projected distance of that fluid parcel from the BH onto the line connecting the stellar center of mass and the BH. This energy dependence implies that perpendicular "slices" of the star return simultaneously to the BH, which allows one to construct a fallback rate Ṁ -being the rate at which bound material returns to the BH-that accounts for the stellar density profile ρ; the result is (Lodato et al. 2009;Gallegos-Garcia et al. 2018;Golightly et al. 2019) where M å is the stellar mass, τ=t/T mb with the return time of the most bound debris, M the BH mass and R å the stellar radius, η=R/R å with R the spherical distance from the center of the star, and r p = ) is the average stellar density. Note that the integral in Equation (1) is negative for all times t<T mb , and hence the physical fallback rate is zero until the most bound debris element returns to the BH.
From Equation (1), within the impulse approximation the only timescale relevant to the problem is the return time of the most bound debris; hence the time at which the fallback reaches any characteristic value (e.g., the time to peak, the time between half-peaks, the time to reach µ a -+ M t 5 3 with α>0) is also a constant multiple-which depends on the stellar structure-of this timescale. If we denote the dimensionless time to any characteristic fallback value by τ c , then by definition the corresponding physical time t c is This simple expression demonstrates the relative importance that stellar structure can have when inferring BH properties from observations 5 : assume that for a given tidal disruption event with known t c , we know the mass and radius of the disrupted star. Then, for the same event if we ascribe to the disrupted star two different density profiles, ρ 1 and ρ 2 , we will infer two different BH masses for the same event, M 1 and M 2 , their ratio being Thus, relatively small differences in the stellar structure are capable of producing more discrepant BH masses owing to the squared dependence in this expression. Figure 2 illustrates the fallback rate onto the BH (in units of Solar masses per year as a function of time in years) from the frozen-in approximation, where the disrupted star has a mass and radius equal to those of the 1 M e , TAMS star (see Table 1) and the BH has a mass of 10 6 M e . Each curve adopts a different stellar density profile, being a γ=5/3 polytrope (blue, dotteddashed), a γ=1.35 polytrope (green, dashed), and the profile resulting from MESA (red, solid). It is clear from this figure that, despite the fact that the bulk stellar properties are the same, the density profile has a marked effect on the shape of the curve. To highlight the differences induced by stellar structure, the points on each curve give the time to different, characteristic values of the fallback curve, being the time to peak (t max ), the time to reach half the peak (which occurs on the rise and the decay of the curve, denoted, respectively, by t half,1 and t half,2 ), and the time to reach µ -M t 4 3 , t 4/3 . While certain characteristic times are visibly comparable for each model (e.g., t max ), other timescales are more noticeably discrepant (e.g., t 4/3 ).
Because of these discrepancies, for a TDE with a given physical timescale, the BH mass required to yield a fallback curve with that timescale will differ according to the stellar structure model that one employs, and the magnitude of the difference in mass will be larger or smaller depending on the characteristic timescale itself. For example, using Equation (4) and letting ρ 1 describe the density profile of the 1 M e , TAMS MESA progenitor, we find M M 2.3 2 1  if one models the density profile by a γ=1.35 polytrope and uses the difference between the time to peak and the first time to half peak (i.e., the timescale t c =t max -t half,1 6 ); physically, one can interpret this result by saying that a γ=1.35 polytrope requires a BH mass ∼2.3 times larger to reproduce the time to peak from the 1 M e , TAMS MESA progenitor. On the other hand, if we use the same timescale but model the star as a γ=5/3 polytrope, then we find M M 14 2 1  -because a γ=5/3 polytrope peaks considerably earlier than the real star, a much larger BH mass is required in order to extend the time to maximum fallback. Figure 3 directly demonstrates how changing the mass of the BH yields commensurate characteristic timescales: the solid red curve shows the same fallback rate as in Figure 2 for the 1 M e , TAMS stellar profile computed with MESA and a 10 6 M e SMBH (note that this figure is on a linear-linear scale). The dashed, green line, and the dotted-dashed blue line show the fallback rates for a γ=1.35 and γ=5/3 polytrope, respectively, again with the same stellar properties as those of the MESA progenitor. However, here we varied the mass of the BH according to the value required to yield the same time to peak, being M=2.3×10 6 for the γ=1.35 polytrope and =Ḿ M 1.4 10 7  for the γ=5/3 polytrope. The time between the first half-peak (marked with a †) and the peak fallback rate (marked with a å) is the same in each case. This shows that, by varying the BH mass appropriately, one can reconstruct the same physical timescale for different stellar properties.
One can repeat this procedure for the nine different models presented in Table 1 and thereby assess the degree to which a polytropic density profile reproduces the fallback curve-and correspondingly the inferred BH mass-of the one obtained with the MESA model. However, the impulse approximation ignores crucial physics of the disruption process that also alter Figure 2. Fallback rate calculated using the frozen-in approximation, in Solar masses per year as a function of time in years, for the 1 M e , TAMS star, when the density profile is modeled as a γ=5/3 polytrope (blue, dotted-dashed) and a γ=1.35 polytrope (green, dashed); the solid red curve shows the fallback rate calculated from the MESA progenitor. Here the BH mass was set to 10 6 M e . The different points show characteristic times in the fallback rate, including the time taken to reach the peak (t max , asterisks), to reach µ -M t 4 3 (t 4/3 , bullets), and the time taken to reach half the peak fallback rate (which occurs on both the rise and the decay of the curve; respectively t half,1 and t half,2 , crosses). It is apparent that, while each one of these stars possesses the same mass and radius, the density profile also plays a large role in generating differences between these characteristic timescales. Note.Each disruption has β=3, and hence the pericenter distance of each star to the SMBH is obtained by dividing the tidal radius (fourth column) by 3. From left to right, the star name, stellar mass M å in solar masses, stellar radius R å in solar radii, and tidal radius In the next section, we employ hydrodynamical simulations to obtain more realistic fallback rates.

Setup
We use the smoothed particle hydrodynamics (SPH) code PHANTOM (Price et al. 2018) to model the hydrodynamics of the disruption process. Following previous work (e.g., Coughlin & Nixon 2015;Coughlin et al. 2016;Wu et al. 2018;Golightly et al. 2019) we model the central SMBH as a Newtonian point mass at the origin with an accretion radius, inside which particles are removed from the simulation. We include the self-gravity of the star, and we model the stellar pressure using an adiabatic equation of state P=Kρ γ , where K is a conserved quantity for each fluid element but can be spatially dependent within the original star. K is chosen such that the fluid pressure is the total pressure in the MESA calculation, and thus the SPH star has an equilibrium density structure that matches the stellar density structure calculated by MESA. The properties of each star are given in Table 1.
We construct the star with particles placed on a close-packed sphere that is stretched to achieve the desired density distribution. We then relax the particle distribution in isolation with a velocity damping force until the star settles into a numerically relaxed configuration. We plot the density structure at this point against the solutions from MESA in the Appendix (Figure 9). We further checked that these solutions are numerically relaxed by evolving them in isolation for the time taken for the star to reach pericenter in each case, and we found no subsequent evolution of the density profile.
We then placed the relaxed stars on parabolic orbits around a central SMBH, with the initial location at five tidal radii. To ensure full disruption of at least some of the stars we chose an calculated from the mass and radius of each star (see Table 1) and r p is the pericenter distance of the stellar center of mass. Thus, while the β is the same for each simulation (and therefore the average strength of the tidal field at pericenter is the same), the physical pericenter is different from simulation to simulation. As we shall see below, β=3 is still not large enough to completely disrupt all of the MESA-generated stars, but we neglected to go to higher β because of the prohibitively large particle number required to achieve converged results. For each simulation we employed 10 6 particles, though we also ran tests with 10 5 particles and found only small differences (at the noise level) in the fallback.
To understand the impact of stellar structure, additional simulations were also performed with a γ=1.35 polytrope and a γ=5/3 polytrope, with the mass and radius of each polytrope matched to those of the MESA progenitor and the same orbital properties (i.e., the same β). For the polytrope disruptions, we set the adiabatic index equal to Γ=5/3, such that the microphysics and the stiffness of the equation of state is identical to that employed in the MESA calculations. In this way, the only difference between the polytrope disruptions and the MESA star disruptions is the density profile-the stellar mass, radius, and microphysics are identical-and these simulations therefore isolate the imprint that the density profile has on the fallback rate.
In this Letter we are interested in the fallback rate, defined as the rate at which disrupted material returns to pericenter, and we therefore increase the accretion radius to 3 r t once the star has passed through pericenter. The fallback rate of material through this radius will, in general, differ from the true accretion rate, which is the rate at which material passes through the horizon of the BH. In general, the latter requires detailed, high-resolution simulations that accurately model the formation of the accretion flow around the BH. This is not currently computationally feasible for standard TDE parameters, but has been attempted by various authors in restricted cases (see, e.g., Hayasaki et al. 2013Hayasaki et al. , 2016Shiokawa et al. 2015;Bonnerot et al. 2016;Saḑowski et al. 2016). For the relatively high-β simulations considered here, the general relativistic advance of periapsis will be large, which should enhance energy dissipation and the formation of an accretion disk. We therefore expect the fallback rates we find to closely track the true accretion rate onto the BH, but we leave a detailed study of this process to a future investigation (for which our fallback rates could be used as inputs).
Finally, as we noted above, the more highly evolved MESA stars have extremely dense cores, and-even for a β of 3 encounter-those cores survive the tidal interaction with the BH. In these instances, the time step is extremely limited because of the high density and sound speed at the center of the stellar remnant, which makes these simulations prohibitively expensive to run for the duration over which the fallback occurs. Therefore, once the surviving core recedes to a significant distance from the SMBH, we replace all of the particles in the bound core that have a density above the maximum (non-core) stream density by a single sink particle. The position and velocity of the sink is set equal to the center of mass position and velocity of all the particles used to create the sink, and the accretion radius of the sink particle is equal to the maximum distance of these particles from their center of mass (the sink position). For the simulation with the densest core (0.3 M e , TAMS) we inserted the sink at a time of 3.5 days post  Figure 2). The dashed green curve and the dotted-dashed blue curve show the fallback rates for a γ=1.35 and γ=5/3 polytrope, respectively, when the stellar properties (mass and radius) are matched to those of the MESA progenitor. By changing the BH mass to the value shown in the legend, we are able to reproduce the time taken to go from half peak to peak, being t c =0.63 yr. For each curve, this timescale is the time taken to go from the † to the å.
pericenter. The other simulations with bound cores were less computationally expensive and ran to later times. As such we inserted the sink at times (post pericenter) of 15 days (1.0 M e , TAMS), 11 days (3.0 M e , TAMS), 173 days (1.0 M e , MAMS), and 27 days (3.0 M e , MAMS). We tested the robustness of this approach by changing the time at which we replaced the resolved core with the sink, and found negligible changes in the subsequent fallback rate.
One way in which this approximation could adversely affect the fallback rate is if we erroneously included the marginally bound (to the core) radius within the particles that constitute the core. If we were to make this error, the fallback rate onto the BH would abruptly terminate at a late, but ultimately finite time. We have checked that the mass of the sink particle increases very slightly after its formation, indicating that the marginally bound radius is indeed outside of the radius of the sink particle. Figure 4 illustrates the fallback rate onto the BH in Solar masses per year as a function of time in years from six different stellar models, with the specific stellar model shown in the legend. In each panel the solid curve shows the fallback from the star with the MESA density profile, the dashed curve is the γ=5/3 model matched to the MESA stellar mass and radius, and the dotted-dashed curve is the γ=1.35 polytropic model (again, with the same mass and radius as the MESA star). It is evident from the top-left panel and the middle-left panel that the 0.3 M e , ZAMS and the 1.0 M e , ZAMS MESA fallback curves are extremely well-reproduced by γ=5/3 and γ=1.35 polytropes, respectively. This finding indicates, correspondingly, that such stars can be very well-modeled by single polytropes that are gas-pressure and (nearly) radiationpressure dominated (see Figure 5). This result for the 1 M e , ZAMS star was also recovered by Goicovic et al. (2019), who found that their fallback rates (from a 1 M e , ZAMS star generated with MESA) were very similar to those of Guillochon & Ramirez-Ruiz (2013), who used a polytropic approximation.

Results
For every other stellar model, however, there are notable differences between the fallback curves from the MESA and the polytropic models. Specifically, we see that employing the MESA density profile over either polytropic model systematically shifts the return time of the most bound debris to earlier times, the time to peak to earlier times, and the magnitude of the peak rate itself is also increased. Furthermore, because the total mass is the same, the larger peak in the accretion rate and the earlier time-to-peak imply that the MESA curves must fall below the polytropic ones at some point, and this is indeed recovered in each case. It is also apparent that the MESA models conform to a power-law decline at an earlier time than do the polytropic models.
Every polytropic star is completely disrupted by the SMBH for these β=3 encounters, which agrees with the results of Guillochon & Ramirez-Ruiz (2013) and Mainetti et al. (2017), who demonstrated that the critical β for the full disruption of a γ=5/3 polytrope is b 0.9  , while that for a γ=4/3 polytrope is b 2  . Each ZAMS MESA progenitor is also fully disrupted. However, more highly evolved stars start to yield partial TDEs, and the MESA 1 M e MAMS and 3.0 M e MAMS leave stellar cores at the location of the marginally bound orbit. We see that, for the case of the 3.0 M e , MAMS progenitor, the presence of the core has the affect of modifying the late-time fallback rate, which declines approximately as ∝ t −9/4 instead of t − 5/3 . This result is in agreement with the analytic predictions of .
The TAMS MESA progenitors all possess extremely dense cores that survive the tidal encounter, each of which modifies the late-time fallback rate onto the BH, as shown in Figure 6 (the time taken for the 3.0 M e progenitor to go from MAMS to TAMS is very short, and hence the TAMS fallback curve appears nearly identical to the bottom-right panel of Figure 4; we therefore opted not to show this fallback rate). We also emphasize that the lifetime of the 0.3 M e star is well in excess of the age of the universe, and hence this star cannot be disrupted by an SMBH (at least not any time soon). However, the density profile of this star could conceivably be achieved by a more massive progenitor, which would evolve to the TAMS within the age of the universe, with different initial properties (e.g., metallicity, rotation).
As we noted above, the 1 M e MAMS MESA progenitor also possesses a bound stellar core at the center of mass of the tidally disrupted debris stream. In this case, however, the star is initially completely disrupted, and the core recollapses out of the stream at a time significantly after the stellar center of mass passes through pericenter. For this reason, the mass of the surviving core is only a small fraction of the initial progenitor ( 15%  ), and consequently the fallback curve shows little evidence of the gravitational influence of the core over ∼1 yr and appears to approach a decline ∝ t −5/3 . Figure 7 shows the fallback from the 1 M e MAMS MESA progenitor out to 10 yr post-disruption and demonstrates, however, that the presence of the core does start to affect the fallback at later times. In particular, we see that while the first year shows little evidence of the existence of a bound core, there is a clear break in the fallback curve at a time of 1-2 yr where the power law of the fallback rate transitions from ∝ t −5/3 to one that is better matched by ∝ t −9/4 . Interestingly, this time at which a break in the power law is exhibited is very close to the time predicted by the analytical model in Coughlin & Nixon (2019; see their Figure 2), and coincidentally also occurs around the same time at which the fallback rate drops below the Eddington limit of the SMBH (dashed black line, assuming a radiative efficiency of 10% and an electronscattering opacity of 0.34 cm 2 g −1 ).

Implications for BH Mass Estimates
It is apparent from Figures 4 and 6 that, for the majority of progenitors, non-polytropic stellar structure generates substantial differences in the fallback rate onto the BH. Notably, the MESA profiles yield earlier times to peak and larger peak fallback rates, and they more rapidly approach a power-law falloff as compared to the polytropes. To use these differences to estimate the corresponding differences in the inferred BH mass that would arise by assuming a given density profile, we must have a mapping between a characteristic timescale (e.g., the time to peak), the BH mass, and the properties of the star. As shown in Section 3, when the fallback rate is computed with the impulse approximation, this mapping arises through the timescale p µ ( )  which is the return time of the most bound debris. Additional dependence on the stellar structure modifies the fallback rate through a dimensionless function of time normalized by T mb , and that dimensionless function can be calculated from the (assumed-unaltered) density profile when the star is at the tidal radius. Thus, the dependence on the BH mass arises only as ∝ M 1/2 , and in this approximation, fallback curves are simply scaled in time and magnitude by M .
By comparing Figure 2 to Figures 4 and 6, wee see that the frozen-in approximation does not accurately reproduce many of the features of the numerically obtained fallback rates. In addition to the fact that the time to peak is shorter and the peak itself is higher in the numerical simulations (by a factor 10), the ordering of the curves is actually inverted between the two The specific star is shown by the name in the legend, and panels on the left side show the fallback from stars at ZAMS, while those on the right are more highly evolved. It is apparent from the top-left and middle-left panels that the fallback curves from the 0.3 M e , ZAMS and the 1 M e , ZAMS progenitors are very well-reproduced by γ=5/3 and γ= 1.35 polytropes, respectively. Every other fallback curve from a MESA-generated density profile, however, shows significant deviations from the polytropic approximations. We also see that the 3.0 M e , MAMS follows ∝ t −9/4 at late times, which results from the presence of a bound core that survives the encounter ; no bound core is left when the star is modeled as a polytrope). The 0.3 M e , MAMS MESA star also shows enhanced variability in the fallback rate, which arises from the fact that the stream-unlike the polytropic models for the same MESA star mass and radius-has fragmented vigorously into small-scale clumps.
approaches: while Figure 2 shows that the γ=5/3 polytrope peaks earlier than γ=1.35 polytrope, which itself peaks earlier than the MESA model, Figures 4 and 6 demonstrate that the γ=5/3 polytrope always reaches a peak after the γ=1.35 polytrope. 7 Moreover, for every case except the 0.3 M e ZAMS and 0.3 M e MAMS progenitors, the MESA model peaks earlier than the γ=1.35 polytrope. The numerically obtained return time of the most bound debris also differs for each density profile, whereas, under the frozenin approximation, this timescale-for the same BH mass-is only affected by the stellar mass and radius (which, for a given star, are identical by construction).
These discrepancies indicate that the impulse approximation does not include enough physics to accurately capture the bulk features of the fallback rate. As discussed at length in Coughlin et al. (2016), it is likely that the most crucial physical ingredient lacking from the impulse approximation is the self-gravity of the debris stream, as the stellar center of mass rapidly recedes outside of the tidal sphere of the BH. At this point, the stellar density is comparable to the "BH density," being r~M r • t 3 , and the selfgravity of the stream is capable of competing against the shear of the BH. The self-gravity of the stream induces density waves that traverse the stream radially, and these waves serve to generate more pronounced "shoulders" near the extremities of the stream and flatten the dm/dr∝ρ curve from the polytropic one that follows from the frozen-in approximation. It is, in fact, because of this nearly flat dm/dr generated by self-gravity that the fallback curves more rapidly approach the t −5/3 decline (or the t −9/4 decline for the partial TDEs). The higher central density of the γ=1.35 polytrope also generates more vigorous density waves, which correspondingly produce a flatter density distribution and give rise to an earlier time-to-peak as compared to a γ=5/3 polytrope.
Nonetheless, it is likely that after some amount of time following the disruption, the mass distribution is approximately frozen-in, meaning that self-gravity has smoothed out any density perturbations and the stream is long enough that the time-dependent potential due to self-gravity is small. 8 In this case, the energies of gas parcels comprising the stream are still Keplerian in the potential of the BH, and the energies themselves are simply established at some later time; indeed, these arguments were used and verified by a direct evaluation of the energy distribution at different times by Lodato et al. (2009) andRamirez-Ruiz (2013) to calculate fallback rates to the BH after only a small fraction of the return time of the most bound debris had been directly simulated. 9 Moreover, if for a given β and stellar progenitor the density profile at the time the energy is frozen-in is independent of the BH mass, which the simulations of Wu et al. (2018) verify 10 (see their Figure 1), then it follows that any physical timescale in a TDE can be written where f å,β is a function that depends only on the stellar properties and β. We therefore see that, for the same star and the same orbital parameters, we recover the same result as we did in Section 3: for two TDEs with identical orbital and stellar properties and physical timescales t 1 and t 2 , we can satisfy  7 This effect can also be seen in Figures 2 and 10 of Lodato et al. (2009), though this inversion was not noted by those authors. 8 However, the arguments of Coughlin et al. (2016) and the simulations of Coughlin & Nixon (2015) suggest that the stream is weakly gravitationally unstable, and hence the freezing of the energy distribution is only valid over an integrated region of the stream that contains many clumps that form out of the instability. 9 When the tidal encounter leaves a surviving core behind, a Keplerian energy distribution is no longer upheld; however, as shown by , when the self-gravity of the stream itself no longer significantly modifies the density distribution along the stream, one can make a change of variables when calculating the fallback rate that shows that Equation (6), and hence Equation (7), still holds. 10 This is also a reasonable expectation, as β measures the tidal strength of the BH; thus, for encounters with the same β, it follows that the energy distribution should be roughly fixed at the same time after self-gravity (which depends only on the stellar properties) has modified the density distribution, and the absolute value of the BH mass should not matter. This assumption breaks down, however, once the orbital timescale becomes shorter than the time over which self-gravity acts to modify the density distribution, which occurs for very small BH masses (where even the tidal approximation itself starts to break down).
we would require a BH mass of 10 5  to reproduce the observed timescale. Table 2 gives the ratio M 2 /M 1 required to shift the timescale of the polytropic star to the timescale reproduced by the MESAstar disruption. The timescale itself is shown in the top row of the table, where t max -t half,1 is the time to go from first-half-max to the peak fallback rate, t half,2 -t max is the time taken to fall by a factor of two below the peak fallback rate, and t half,2 -t half,1 is the FWHM of the fallback curve. The stellar progenitor is given in the left-most column of the table, and the value in the left (right) of each cell is the ratio M 2 /M 1 required to yield the MESA-generated timescale by modeling the star as a γ=1.35 (γ=5/3) polytrope. For example, if we were to model the 0.3 M e star as a γ=1.35 polytrope, then we would require a BH mass of M 2 =5×M 1 to reproduce the timescale t max -t half,1 found from the disruption of the MESA model, and the number M 2 /M 1 =5 is shown in the top-left cell of the table.
We see from this table that, for stellar density profiles that are well-reproduced by polytropes (the 0.3 M e ZAMS and the 1 M e ZAMS stars, as shown in Figure 5), the mass ratio required to reproduce the MESA-generated timescale with a polytropic one is very close to unity. In these instances, the inferred BH mass estimate is not far off the true, underlying value. However, for more massive progenitors and more highly evolved stars, the extremely dense core of the MESA progenitor shifts each timescale earlier, which consequently requires a significantly smaller BH mass to yield the same characteristic timescale with a polytropic model.
As a direct demonstration of this effect, Figure 8 illustrates the fallback rate from the disruption of the 3 M e , ZAMS MESA model by an SMBH of mass M 1 =10 6 M e , which is the same curve shown in the bottom-left panel of Figure 4. The dotteddashed curve is the fallback curve from the disruption of a γ=1.35 polytrope, the stellar mass and radius identical to those of the MESA star. In this case, however, we set the mass of the disrupting SMBH to M 2 =2×10 5 , which is, as seen in Table 2, the value predicted to equate the time to go from the first-half-max to the peak between the two models, and we aligned the first time to half-max of the polytrope fallback curve to that of the MESA star (i.e., we initially "see" both disruptions at the same time, that time being the first time to half-max). We see that this polytropic model provides an extremely good fit to the "data" obtained from the fallback curve of the MESA model, but at the expense of incorrectly inferring the SMBH mass by nearly an order of magnitude.
Of course, our approach here to "modeling" the light curve of a TDE is overly simplistic, as one does not necessarily have any prior information about the nature of the progenitor or the BH, and one must use a combination of timescales to recover the best-fit model parameters (as is done in, for example, Guillochon et al. 2018). However, these results do demonstrate In each case the MESA density profile yields a bound core that survives the encounter, while the polytropes do not; this results in a late-time power-law falloff that declines approximately as t −9/4 (dotted-dotted-dashed line) and is significantly steeper than t −5/3 (dotted line). The Eddington luminosity of the BH, assuming a radiative efficiency of 10% and an electron-scattering-dominated opacity of 0.34 cm 2 g −1 , is shown by the long-dashed black line. The dotted line shows the scaling ∝ t −5/3 (predicted to be the asymptotic power law followed by the fallback if there is no surviving core), the dotted-dotted-dashed line gives the scaling ∝ t −9/4 (predicted by Coughlin & Nixon 2019 to be the asymptotic power-law decline of the fallback if there is a surviving core), and the long-dashed line gives the Eddington luminosity of the hole if the radiative efficiency is 10% and the opacity is set to the electronscattering opacity of 0.34 cm 2 g −1 . Here the stream possesses a gravitationally bound core (with a mass of ∼15% of the progenitor star) at the location of the marginally bound radius that reforms out of the stream after the star is initially completely disrupted. We see that for roughly the first year after the return of the most bound debris there is little evidence of the existence of the core on the fallback, and the fallback curve appears to asymptote to a t −5/3 decline. However, around roughly 1 yr, there is a noticeable break in the falloff, and the curve steepens to a decline that is well matched by the power law ∝ t −9/4 . that care needs to be taken to ensure that the template fallback curves used to interpret observed data sets contain a sufficiently broad range of stellar density profiles (e.g., accounting for stellar age, and non-polytropic density profiles) to ensure that the inferred parameters are accurate and that the error bars that result from the data fitting are appropriate.

Discussions and Conclusions
In this Letter we presented simulations of the tidal disruption of stars encountering an SMBH. We modeled the tidally disrupted stars with the same bulk properties (mass and radius) using three different prescriptions: (1) a γ=5/3, polytropic density profile, (2) a γ=1.35(≈4/3), polytropic density profile, and (3) a density profile calculated from the MESA stellar evolution code. Our simulated disruptions included stars with masses of 0.3 M e , 1.0 M e , and 3.0 M e , each of which was evolved to the ZAMS, the TAMS (where the hydrogen mass fraction in the core fell below 0.1%), and the "middle-age main sequence," which we defined to be the time at which the hydrogen mass fraction in the core fell below 20%. We therefore simulated a total of 27 disruptions (nine different stars, each star modeled with a MESA density profile and two different polytropic profiles).
In each simulation we maintained the same physics: we employed a polytropic equation of state where the Lagrangian entropy was fixed and set to ensure the isolated star was in hydrostatic balance (see Figure 9); we fixed the adiabatic index in each simulation to 5/3, which ensures that the dynamics of the stream's self-gravity (see Coughlin & Nixon 2015) differs only by the mass-entropy distribution along the stream; and we fixed the tidal effects from the BH on each star by employing a b º = r r 3 t p for each simulation, where r p (r t ) is the pericenter (tidal) radius of the star (see also Table 1). Our simulations therefore isolate the impact of the stellar density profile calculated from MESA when compared to those calculated by γ=5/3 and γ=1.35 polytropes.
In general we find that there are significant differences in the simulated fallback rates for stars with different masses and different ages, and further that in most cases these fallback rates deviate significantly from predictions made using polytropes. The exceptions, which come as no surprise as their structures are accurately modeled by polytropes (see Figure 5), are the 0.3 M e ZAMS star-which is well-modeled by a γ=5/3 polytrope-and the 1.0 M e ZAMS star-which is wellmodeled by a γ=1.35 polytrope. At both MAMS and TAMS we find that neither polytrope provides an acceptable description for the fallback curve. Similarly the fallback rates from the 3.0 M e stars are all significantly different to the fallback rates from either polytrope.
There are also differences found in the overall dynamics of the disruption event. In several cases, most notably for the 0.3 M e MAMS star, the debris stream is significantly more self-gravitating for the MESA star than for the polytropes. This results in the fallback curve exhibiting more variability on the power-law decay (see the top-right panel of Figure 4). We also find that several of the MESA stars are not fully disrupted, even with an impact parameter of β=3, and this arises from the more centrally concentrated nature of the MESA stars. 11 The difficulty of fully disrupting real stars, and particularly more highly evolved stars, implies that a greater fraction of the events we observe will be partial, rather than full, disruptions (though we note that stars spend more of their lives near the ZAMS, where the MESA profiles still yield full disruptions for β=3). It has recently been shown  that TDEs that leave a bound core have a fallback rate whose power-law index asymptotes to ≈−9/4 rather than the usual −5/3. In each case that leaves a bound core, our simulations recover this result. In future work we will explore whether simulations are also capable of recovering, e.g., the time at which the power-law slope changes to this value as a function of the mass of the core that survives the encounter. We find (see also Guillochon & Ramirez-Ruiz 2013) that in some cases, here for the 1.0 M e MAMS MESA star, that the core can be initially fully disrupted but reform after leaving the tidal radius. We attribute this to the velocity field imparted in the stellar debris by the BH tides, which can at later times cause the stream to converge along its width and augment the density within the stream (Coughlin et al. 2016;Steinberg et al. 2019).
In Section 5 we described the impact of using realistic stellar models in simulations of TDEs on the inference of the BH mass from observed TDE light curves. We showed that for many types of stars, and particularly those that are more massive at the ZAMS or more highly evolved, the shape of the fallback curves from simulations that employ polytropic stellar models can lead to large errors in the estimated BH mass. Therefore, employing accurate models for the imprint of stellar structure on the fallback rate to produce templates for TDE light curves appears essential for accurately inferring the BH mass. It could also be that degeneracies between the parameters in TDEs (e.g., stellar mass, stellar age, stellar spin, stellar metallicity, BH mass, BH binarity, BH spin, impact parameter, inclination and phase angles of the stellar orbit, and the accretion dynamics on small scales) mean that an unambiguous estimation of system parameters (or at least an understanding of the true level of error within those estimates) requires accurate and detailed modeling of the stars that are being disrupted.