Eddington envelopes: The fate of stars on parabolic orbits tidally disrupted by supermassive black holes

In the classical picture of tidal disruption events (TDEs), the debris from the tidal disruption of a star on a parabolic orbit by a supermassive black hole (SMBH) rapidly circularises to form an accretion disc via relativistic apsidal precession (Rees, 1988). The predicted mass return rate of debris (Phinney, 1989) is $\propto t^{-5/3}$ and the light curve is assumed to be powered by accretion and to follow the same decay.

This picture alone does not predict several properties of observed TDEs, mainly related to their puzzling optical emission (van Velzen et al., 2011; van Velzen, 2018; van Velzen et al., 2021). These properties include: i) low peak bolometric luminosities (Chornock et al., 2014) of $\sim 10^{44}$ erg/s, $\sim\,$1 per cent of the value expected from radiatively efficient accretion (Svirski et al., 2017); ii) low temperatures, more consistent with the photosphere of a B-type star than with that of an accretion disc at a few tens of gravitational radii ($R_{g}\equiv GM_{\mathrm{BH}}/c^{2}$) (Gezari et al., 2012; Miller, 2015), and consequently large emission radii, $\sim\,$$10$–$100$ au for a $10^{6}M_{\odot}$ black hole (Guillochon et al., 2014; Metzger & Stone, 2016); and iii) spectral line widths implying gas velocities of $\sim\,$$10^{4}$ km/s, much lower than expected from an accretion disc (Arcavi et al., 2014; Leloudas et al., 2019; Nicholl et al., 2019).

As a consequence, numerous authors have proposed alternative mechanisms for powering the TDE lightcurve, via either shocks from tidal stream collisions during disc formation (Lodato, 2012; Piran et al., 2015; Svirski et al., 2017; Ryu et al., 2023; Huang et al., 2023), or the reprocessing of photons through large scale optically thick layers, referred to as Eddington envelopes (Loeb & Ulmer, 1997), super-Eddington outflows (Strubbe & Quataert, 2009), quasi-static or cooling TDE envelopes (Roth et al., 2016; Coughlin & Begelman, 2014; Metzger, 2022) or mass-loaded outflows (Jiang et al., 2016; Metzger & Stone, 2016). Recent spectro-polarimetric observations suggest reprocessing in an outflowing, quasi-spherical envelope (Patra et al., 2022).

The wider problem is that few calculations exist that follow the debris from disruption to fallback for a parabolic orbit with the correct mass ratio. The challenge is to evolve a main-sequence star on a parabolic orbit around a SMBH from disruption and to follow the subsequent accretion of material (Metzger & Stone, 2016). The dynamic range involved when a $1M_{\odot}$ star on a parabolic orbit is tidally disrupted by a $10^{6}M_{\odot}$ SMBH is greater than four orders of magnitude: the tidal disruption radius is 50 times the gravitational radius, where general relativistic effects are important, while the apoapsis of even the most bound material is another factor of 200 further away. This challenge led previous studies to consider a variety of simplifications (Stone et al., 2019): i) reducing the mass ratio between the star and the black hole by considering intermediate mass black holes (Ramirez-Ruiz & Rosswog, 2009; Guillochon et al., 2014); ii) using a Newtonian gravitational potential (Evans & Kochanek, 1989; Rosswog et al., 2008; Lodato et al., 2009; Guillochon et al., 2009; Golightly et al., 2019), pseudo-Newtonian (Hayasaki et al., 2013; Bonnerot et al., 2016) or post-Newtonian approximations (Ayal et al., 2000; Hayasaki et al., 2016); iii) simulating only the first passage of the star (Evans & Kochanek, 1989; Laguna et al., 1993; Khokhlov et al., 1993; Frolov et al., 1994; Diener et al., 1997; Kobayashi et al., 2004; Guillochon et al., 2009; Guillochon & Ramirez-Ruiz, 2013; Tejeda et al., 2017; Gafton & Rosswog, 2019; Goicovic et al., 2019); and iv) assuming stars initially on bound, highly eccentric orbits instead of parabolic orbits (Sadowski et al., 2016; Hayasaki et al., 2013, 2016; Bonnerot et al., 2016; Liptai et al., 2019; Hu et al., 2024).

These studies have, nevertheless, provided useful insights into the details of the tidal disruption process. In particular, it has been shown that the distribution of orbital energies of the debris following the initial disruption is roughly consistent with $dM/dE=$const, consistent with the analytic prediction of a $\propto t^{-5/3}$ mass fallback rate, although the details can depend on many factors such as stellar spin, stellar composition, penetration factor and black hole spin (Lodato et al., 2009; Kesden, 2012; Guillochon & Ramirez-Ruiz, 2013; Golightly et al., 2019; Sacchi & Lodato, 2019). The importance of general relativistic effects in circularising debris has also been demonstrated. The self-intersection of the debris stream, which efficiently dissipates large amounts of orbital energy, is made possible by relativistic apsidal precession (Hayasaki et al., 2016; Bonnerot et al., 2016; Liptai et al., 2019; Calderón et al., 2024). But until recently debris circularisation has only been shown for stars on bound orbits, with correspondingly small apoapsis distances and often deep penetration factors (we define the penetration factor as $\beta\equiv R_{\mathrm{t}}/R_{\mathrm{p}}$, where $R_{\mathrm{t}}=R_{*}(M_{\mathrm{BH}}/M_{*})^{1/3}$ is the tidal radius and $R_{p}$ is the pericenter distance).

Recent works have shown that circularisation and initiation of accretion is possible in the parabolic case, by a combination of energy dissipation in the ‘nozzle shock’ that occurs on second pericenter passage (Steinberg & Stone 2024; but see Bonnerot & Lu 2022 and Appendix E for convergence studies of the nozzle shock) and/or relativistic precession leading to stream collisions (Andalman et al., 2022).

In this paper, we present a set of simulations that self-consistently evolve a one solar mass polytropic star on a parabolic orbit around a $10^{6}$ solar mass black hole from the star’s disruption to circularization of the returning debris and then accretion. We follow the debris evolution for one year post-disruption, enabling us to approximately compute synthetic light curves which appear to match the key features of observations.

We modelled the disruption of stars with mass $M_{*}=1M_{\odot}$ and radius $R_{*}=1R_{\odot}$ by a supermassive black hole with mass $M_{\mathrm{BH}}=10^{6}M_{\odot}$ using the general relativistic implementation of the smoothed particle hydrodynamics code Phantom (Price et al., 2018) as described in Liptai & Price (2019).

We used the Kerr metric (Kerr, 1963) in Boyer-Lindquist coordinates, with the black hole spin vector along the $z$-axis, and unless otherwise specified assumed an adiabatic equation of state for the gas with $\gamma=5/3$. We deleted particles that fall within a radius of $R_{\mathrm{acc}}=5R_{g}$ (within the innermost stable circular orbit of a non-spinning black hole) from the simulation to avoid small timesteps close to the horizon.

We placed the star on a parabolic orbit initially at a distance of $r_{0}=10R_{\mathrm{t}}$ with velocity $v_{0}=\sqrt{2GM_{\mathrm{BH}}/r_{0}}$, where $R_{\mathrm{t}}$ is the tidal radius. The parabolic trajectory was oriented such that the star travels counter-clockwise in the $x$-$y$ plane, and that the Newtonian pericentre $R_{\mathrm{p}}$ was located on the negative $y$ axis. For inclined trajectories, we rotated the initial position and velocity vectors by an angle $\theta$ about the $y$-axis.

We modelled the star as a polytropic sphere with polytropic index $n=3/2$ using $4,188,898$ equal mass particles at our highest resolution and approximately $512$K particles in our medium resolution calculations (real stars are expected to be more concentrated, which may alter the fallback rate, delaying the peak with respect to a polytropic model; e.g. Law-Smith et al. 2019), and include the perturbation in the metric due to Newtonian self-gravity in our simulations to hold the star together during the pre-disruption phase.

For our main simulations we employed an adiabatic approximation, implying that energy is trapped or advected rather than radiated. Appendix A reports a set of isentropic calculations where we assume shock heating to be immediately radiated away. The adiabatic simulations capture all accretion feedback down to $5R_{\rm g}$ inside which particles are assumed to accrete without radiating.

Our point of closest approach is typically smaller than the Newtonian pericentre $R_{\mathrm{p}}$ due to relativity, and thus the ‘true’ penetration factor is larger (by up to $\sim\,$1 in our $\beta=5$ calculation, where the true $\beta$ is closer to 6).

Figure 1 shows snapshots of the long-term debris evolution in the $\beta=1$, high-resolution, adiabatic simulation with a non-rotating black hole. The star takes approximately 7 hours to reach pericentre at which point it is tidally disrupted and forms a long, thin debris stream. Approximately half of the stream is bound to the black hole while the other half is unbound (visible as the thin, dense upward stream in the figure). The head of the stream (i.e., the most bound material) takes $\sim\,$26 days to return to pericentre. After the stream passes through pericentre for a second time, the small amount of apsidal precession leads to stream self-intersection at the ‘outer shock’ (most evident in inset panels shown at 60 and 90 days; see Appendix B for more detail on the outer shock), which cause material to eventually fall towards the black hole, driving the formation of kinetic outflows which spread almost isotropically, eventually encasing the black hole (main panel). A fraction of the material joins the returning debris and returns to pericentre, however most of it goes into a low density mechanical outflow. Following the fallback for one year post disruption, we find that the bubble continues to expand, with only a small amount of material forming an eccentric disc (lower panels of Figure 1). The large-scale debris structure (main panel) is similar to the ‘Eddington envelope’ first proposed by Loeb & Ulmer (1997) except that ours is outflowing rather than static, as predicted by subsequent authors (Strubbe & Quataert, 2009; Jiang et al., 2016; Metzger & Stone, 2016).

Figure 2 shows the density (top), temperature (middle) and radial velocity (bottom) distributions of particles a year after disruption, with temperature calculated from Equation (3). The density of the outflowing envelope is $\sim\,$$10^{-16}$ g cm^-3, approximately six orders of magnitude lower than the density of the stream, which is $\sim\,$$10^{-10}$ g cm^-3. The typical outflow velocity is $\sim\,$10⁴ km s^-1. After one year there is approximately 67% ($0.67$ M_⊙) of the stellar mass in the incoming or unbound stream, 31.2% ($0.31$ M_⊙) in the envelope, 0.8% ($0.008$ M_⊙) in the eccentric disc, while 1% ($0.01$ M_⊙) has been accreted by the black hole.

The top panel of Figure 3 shows the energy dissipation rate in our simulation, $L_{\mathrm{shock}}$. This was computed by computing the time derivative of the thermal energy caused by irreversible heating terms due to viscosity and shocks. We find that it is dominated by regions close the accretion radius $R_{\mathrm{acc}}$. The energy dissipation rate rises rapidly from $\sim\,$26 days, when accretion begins, until $\sim\,$50 days after disruption when it reaches a peak of $\sim\,$$10^{44}$ erg s^-1 (roughly the Eddington limit). This is followed by a power law decay that does approximately match the expected $t^{-5/3}$ energy dissipation from the mass fallback rate (Rees, 1988; Phinney, 1989), but with $\dot{M}_{\rm acc}$ lower than the mass fallback rate by approximately two orders of magnitude. We compare this to the luminosity derived from our measured mass accretion rate $\dot{M}$, given by $L_{\mathrm{acc}}=\epsilon\dot{M}c^{2}$. Here, $\epsilon=0.2$ is the efficiency based on our inner accretion radius. However, since our $R_{\rm acc}$ is inside the last stable orbit, matter crossing this radius is mostly advected into the black hole rather than radiated, so converting it into a luminosity simply measures the amount of energy that would have been dissipated to reach this radius if all the matter were on a perfectly circular orbit. The resulting ‘luminosity’ $L_{\mathrm{acc}}$ follows the same trend as the actual energy dissipation rate above, but is $\sim$0.5–1 orders of magnitude greater throughout the simulation. This suggests some material is accreted almost radially, with little to no orbital energy dissipated before accretion (Svirski et al., 2017).

The lower panel of Figure 3 shows the approximate time taken for photons to diffuse from near the black hole out to the photosphere radius, assuming the material remains static. We computed this via

by summing the contributions from each SPH particle intersecting a line-of-sight along each coordinate axis with opacity $\kappa=0.4$ cm²/g. The typical photon diffusion time varies between $\sim\,$10–200 days, peaking approximately 80 days after disruption and decaying roughly as a power law. This indicates that, although most energy is dissipated in regions close to the black hole (see Figure 9), where the reservoir of gravitational binding energy is largest, this material is optically thick. The energy is transported outward by mechanical outflows before some of it is radiated through a photosphere.

Figure 4 shows the temperature at the ‘photosphere’ of the outflowing envelope, produced by ray-tracing through our simulation as described in Sec. 2.1. An expanding emission region of 10–100 au is consistent with those inferred from TDE observations (van Velzen et al., 2021). Photospheric temperatures are in the range expected to produce blackbody emission at UV/optical wavelength. Figure 2 shows that temperatures would remain in the 10⁴–10⁵ K range even accounting for scattering from deeper layers (see Appendix C).

Figure 5 compares the corresponding time evolution of our measured $L_{\mathrm{bb}}$, $R_{\mathrm{photo}}$ and $T_{\mathrm{bb}}$ (inferred bolometric luminosity, blackbody radius and temperature of the blackbody fit to optical wavelengths, respectively) from six different lines of sight (solid lines) to the evolution of observed TDEs inferred from black body fits taken from van Velzen et al. (2021) and Neustadt et al. (2020) (dots). The equivalent plot for the more penetrating $\beta=5$ encounter is shown on the right.

The typical luminosity (top panel) is around $10^{43}$–$\,10^{44}$ erg s^-1 ($10^{9.4}$–$\,10^{10.4}$ $L_{\odot}$), however we caution that the exact value is sensitive to our photospheric temperature (since $L\propto T^{4}$). Nonetheless, our approximate luminosities are consistent with the observed values (see figure 9 of Holoien et al. 2019), as well as with the analytically computed light curves of Lodato & Rossi (2011) (see left panel of their Figure 7). The bolometric luminosity of the envelope in our models appears almost constant, and close to the Eddington value, key features of the analytical model of Loeb & Ulmer (1997).

The photospheric radius (centre panels), when computed by fitting to synthetic spectra in the optical band (Appendix C), is relatively insensitive to the line-of-sight, giving typical radii of $\sim\,$10–100 au ($10^{14}$–$10^{15}$ cm), in the same range as the observations. The photospheric radius reaches a peak after roughly 150–200 days in the $\beta=1$ calculation, and $50$ days in the $\beta=5$ calculation, the latter more closely matching the typical $20-50$ day rise time seen in most of the observations. On the other hand some TDE candidates such as ASASSN-18jd (Neustadt et al., 2020) (see legend) do show a more or less flat luminosity evolution over a 1 yr timescale, similar to our $\beta=1$ lightcurves. The transients AT2018hco and AT2019iih show the closest match to the photospheric radius evolution in our $\beta=5$ calculation, particularly for $t\gtrsim 50$ days.

We find photospheric temperatures (lower panels) in the range of $1$–$4\times 10^{4}K$, with the same range as the observations, but with slightly lower temperatures in our models ($T\sim 10^{4.1}$K) compared to the observations ($T\sim 10^{4.2}$–$10^{4.6}$K) for $t\gtrsim 50$ days. This may be a result of our underestimate of the photospheric temperature due to scattering of photons from deeper layers (Appendix C). The observed trend of temperatures slowly rising as a function of time for $t\gtrsim 50$ days is reproduced in our $\beta=5$ calculation (right panels of Figure 5). Our temperature evolution is similar to that obtained by Mockler et al. (2019) who fitted models of a ‘dynamic reprocessing layer’ to observed TDEs (see their Figure 8), finding a 10–20 day sharp rise in temperature followed by a dip and subsequent flat or slowly rising evolution.

A notable feature in both the $\beta=1$ and $\beta=5$ cases are the almost flat lightcurves predicted at late times along most lines of sight. This is also observed (e.g. Holoien et al., 2019; van Velzen et al., 2021). In our simulations this corresponds to the stalling of the photospheric radii in an expanding gas cloud; however, as discussed below, models without cooling may not accurately describe light curves after $\sim 1$ year.

Our simulations address the main challenge of tidal disruption events: to simulate the disruption of a $1M_{\odot}$ star, on a parabolic orbit, by a $10^{6}M_{\odot}$ black hole in general relativity, and to model the subsequent accretion flow from first principles.

Our main new insight from doing so is that, as in Hu et al. (2024), we have been able to self-consistently form the hypothesised ‘Eddington envelope’ for the first time, more or less as predicted by Loeb & Ulmer (1997) and Metzger & Stone (2016). With the disrupted star moving on an orbit with a mean energy of zero and low angular momentum, we found dissipation of energy and angular momentum to occur through relativistic precession and ensuing stream-stream collisions (see Appendix A), as found by previous authors for stars on bound trajectories (Hayasaki et al., 2016; Bonnerot et al., 2016; Stone et al., 2019; Liptai et al., 2019) and similar to the process originally envisioned by Rees (1988).

Once accretion commences in the optically thick environment with inefficient cooling (we assumed adiabatic gas), the gravitational energy of accreting material launches mechanical outflows, yielding a large, asymmetric, optically thick outflowing envelope (Fig. 4).

We found a similar outflow in recent calculations of initially eccentric encounters (Hu et al., 2024), but with faster rise of the lightcurve due to the shorter orbital periods. The main difference in the parabolic case is that only $\sim 1\%$ of the star is accreted, compared to more than half the star for the eccentric encounters (Hu et al., 2024). The lower efficiency of accretion occurs because of the smaller amount of bound material reaching the inner regions (Appendix A).

The envelope structure is similar to other rapidly expanding, optically thick, quasi-spherical outflows powered by a central energy source, such as supernovae and kilonovae. In this case, the expanding envelope is powered by gravitational energy, dissipated in the radiatively inefficient accretion flow onto the black hole (see Figure 9). Our black hole mass accretion rate, shown in Figure 3, is 1–1.5 orders of magnitude lower than the predicted mass fallback rate. This broad conclusion on the formation of an Eddington envelope is unchanged by whether or not the black hole is spinning, or by whether the star approaches closer to the black hole (see Figure 11), or the numerical resolution (see Appendix E).

That our outflowing envelope, or ‘Black Hole Sun’, is optically thick is demonstrated by our estimated photon diffusion timescales of $\sim\,$10–100 days shown in the bottom panel of Figure 3. This means that the debris produced by the disruption cannot cool efficiently, leading to the formation of a photosphere. Synthetic spectral energy distributions (Appendix C) show this leads to ratios of optical to X-ray emission of order $10$–$10^{3}$, similar to those observed (van Velzen et al., 2021).

Our numerical experiments have several serious limitations in how well they represent real tidal disruptions. The main ones are that we employed an adiabatic approximation, meaning that none of the energy produced is actually radiated in our simulations, that we assume $\gamma=5/3$ and used a polytrope not a real star.

Our adiabatic assumption is likely a good approximation near the peak of the light curve since the photon diffusion timescale is long, although it does not account for the change to the radiation-pressure dominated regime after shock heating, which implies that the adiabatic index should change to $4/3$. On the other hand, the adiabatic assumption likely breaks down by around 1 year after disruption. At t $\gtrsim$ 1 yr, our estimated luminosity (Figure 5) becomes comparable to the energy input rate from shock heating (Figure 3), while the photon diffusion timescale becomes comparable to the time for the kinetic outflows to reach the photosphere, allowing material to radiatively cool.

Therefore, at t $\gtrsim$ 1 yr, with relatively efficient diffusion, the observable luminosity should transition to tracking the shock heating rate, which in turn roughly tracks the mass return rate with $\propto t^{-5/3}$ scaling (top panel of Figure 3). Efficient cooling after 1 year may also lead to thin disc formation (see Appendix A), rather than the eccentric thick disc formed at $t\gtrsim 240$ days in Figure 1. Cooling may also cause our photospheres to retreat at later times, better matching the photospheric radius and temperature observations.

Despite this limitation, our photospheric radius estimates yield a luminosity evolution with the same order of magnitude as optical/UV light curves from several recently observed TDEs (see Figure 5). We do not expect that the dynamics of our simulations would be significantly affected by allowing the photosphere to radiate at early times, mainly because the energy release via radiation is small — most of the energy in our Eddington envelope is released via mechanical outflows, with typical outflow velocities of $\sim\,$1–2$\times 10^{4}$ km s^-1 (see bottom panel of Figure 2). Such velocities match those observed in spectral lines (e.g. Leloudas et al., 2019).

We find that the optical light curve is sensitive to the observing angle, as suggested by Dai et al. (2018a), but due to different photospheric temperatures along different lines of sight (see Figure 5) rather than the presence of a disc. Mildly supersonic outflows indicate that temperatures are not in equilibrium across the photosphere.

The photospheric radii plateau at $\sim 20$ to $\sim 100$ au, depending on the line of sight and the amount of material participating in the outflow (Figure 5). Applying a single-zone model with uniform density $\rho=\frac{4\pi}{3}\frac{M}{R^{3}}$ to a homologously expanding envelope, this stalling is estimated to occur when the photospheric radius is $\frac{2}{3}\sqrt{\frac{\kappa M}{4\pi}}$, where the mass $M$ in the envelope depends on the penetration factor $\beta$. This estimate yields a photospheric radius plateau of $\sim 100$ au for our choice of $\kappa$ and $M=0.1M_{\odot}$. While our predicted radius evolution might change if we allowed radiation to escape (especially at these late times), a significant late-time flattening has been observed in optical/UV for black holes of this mass (van Velzen et al., 2019).

Many TDE candidates show soft X-ray emission. While in general the ratio of X-ray to optical luminosities in TDE is low (Auchettl et al., 2017), there are cases where optical and X-ray emission are both present. ASASSN-14li has similar X-ray and optical luminosities (Holoien et al., 2016), though with X-rays lagging by 32 days (Pasham et al., 2017), which is difficult to explain via reprocessing. Our model can help to explain some of these TDE features. At late times, as the material cools and the photosphere retreats to reveal regions close to the black hole, we expect the formation of a disc or a thick torus similar to that described in Appendix A, which will be a significant central X-ray source.

Indeed, we found in our $\beta=5$ calculation (Appendix C) that retreat of the photosphere produces a rising X-ray luminosity similar to some observed TDEs Gezari et al. (2017) and to the class of ‘veiled’ TDEs (Auchettl et al., 2017). Our Eddington envelope naturally obscures any central X-ray source, with absorption strongly varying for different lines of sight with typical $N_{H}\sim 10^{24}$–$10^{26}$ cm^-2. This suggests an explanation for the dichotomy between X-ray and optical TDE in terms of different viewing angles to the source, in a sort of ‘unified’ model for TDE emission (Dai et al., 2018b).

In summary, we simulated the tidal disruption of 1M_⊙ polytropic stars around a $10^{6}$ solar mass black hole, using general relativistic hydrodynamics in the Kerr metric, following the simulations for one year post disruption. Our main finding, similar to our recent finding for real stars on eccentric orbits (Hu et al., 2024) is the production of an optically thick ‘Eddington envelope’ hypothesised by previous authors (e.g. Loeb & Ulmer, 1997; Metzger & Stone, 2016). Our calculations predict:

1. 1.

   outflow velocities of $\sim\,$$10^{4}$ km s^-1;

2. 2.

   large photosphere radii of $\sim\,$$10$–$100$ au;

3. 3.

   a relatively low mass accretion rate ($\sim\,$$10^{-2}M_{\odot}/$yr) mediated by the mechanical outflow of material;

4. 4.

   peak luminosities of $\sim\,$$10^{43}$–$10^{44}$ erg/s.

We find that our $\beta=5$ encounters provide a better match to TDE observations than $\beta=1$, and that black hole spin does not significantly change the outcome. We thus confirm Eddington envelopes as the likely solution to the puzzle of TDE optical/UV emission.