Dynamical Unification of Tidal Disruption Events

About a hundred tidal disruption events (TDEs) have been observed and they exhibit a wide range of emission properties both at peak and over their lifetimes. Some TDEs peak predominantly at X-ray energies while others radiate chiefly at UV and optical wavelengths. While the peak luminosities across TDEs show distinct properties, the evolutionary behavior can also vary between TDEs with similar peak emission properties. At late time, some optical TDEs rebrighten in X-rays, while others maintain strong UV/optical emission components. In this Letter, we conduct three-dimensional general relativistic radiation magnetohydrodynamics simulations of TDE accretion disks at varying accretion rates ranging from a few to a few tens of the Eddington accretion rate. We make use of Monte Carlo radiative transfer simulations to calculate the reprocessed spectra at various inclinations and at different evolutionary stages. We confirm the unified model proposed by Dai et al. (2018), which predicts that the observed emission largely depends on the viewing angle of the observer with respect to the disk orientation (X-ray strong when viewed face-on and UV/optically strong when viewed disk-on). What is more, we find that disks with higher accretion rates have elevated wind and disk densities, which increases the reprocessing of the high-energy radiation and thus generally augments the optical-to-X-ray flux ratio along a particular viewing angle. This implies that at later times, as the accretion level declines, we expect the funnel to effectively open up and allow more X-rays to leak out along intermediate viewing angles. Such dynamical model for TDEs can provide a natural explanation for the diversity in the emission properties observed in TDEs at peak and along their temporal evolution.


INTRODUCTION
The tidal disruption of stars by massive black holes (MBHs) offers a unique probe of MBH demographics (Mockler et al. 2019;Gezari 2021), host galaxy properties (French et al. 2020), stellar dynamics (Stone et al. 2020;Pfister et al. 2020), as well as black hole accretion and jet physics (Dai et al. 2021). When a star with mass m and radius r approaches a MBH with mass M BH , the star is disrupted within the tidal radius R t ≈ (M BH /m ) 1/3 r , where the MBH's tidal force exceeds the stellar self-gravity. About half stellar debris orbits back to the vicinity of the MBH following a characteristic pattern which first quickly increases to a peak and then declines with time following a mass fallback rateṀ fb ≈ t −5/3 (Rees 1988;Evans & Kochanek 1989;Guillochon & Ramirez-Ruiz 2013). When M BH ∼ 10 6 M ,Ṁ fb can exceed the Eddington accretion rate by two orders of magnitude at peak and stays super-Eddington for years after peak. Here we define the Eddington accretion rate of a black hole asṀ Edd = L Edd /(η NT c 2 ), where L edd = 4πGM BH c/κ is the Eddington luminosity for an opacity κ, and η NT is the nominal accretion efficiency for the Novikov-Thorne thin disk solution (Novikov & Thorne 1973).
With the recent launches of all-sky transient surveys including ZTF and eROSITA, the number of detected TDE candidates has reached around one hundred (e.g., Gezari 2021;Sazonov et al. 2021;Hammerstein et al. 2022). TDEs have been categorized into two classes based on their main emission type at peak: X-ray TDEs (Auchettl et al. 2017;Saxton et al. 2021) and optical TDEs (van Velzen et al. 2020). Most X-ray-selected TDEs emit soft X-rays with effective temperatures at 10 5 − 10 6 K, while only three of them emit beamed, hard X-rays which are associated with relativistic jets (e.g., Bloom et al. 2011;Burrows et al. 2011;Cenko et al. 2012;De Colle et al. 2012). The optical TDEs have lower effective temperatures at few×10 4 K, and are further characterized by their spectroscopic features (producing strong and broad H, or He, or Bowen fluorescence emission lines) (Leloudas et al. 2019;Charalampopoulos et al. 2022). Interestingly, as the luminosity of a TDE typically declines by around one order of magnitude over tens to hundreds of days after peak, its effective temperature usually undergoes a peculiar non-evolution. Recently, a large number of TDEs have been followed up for longer than a year and they show different late-time behaviors. A few initially optically-strong TDEs rebrighten in X-rays (e.g., Gezari et al. 2017;Holoien et al. 2018;Wevers et al. 2019;Hinkle et al. 2021;Liu et al. 2022), while many others maintain strong UV/optical emissions over years .
Many theoretical models have been proposed to explain these TDE emission properties (Roth et al. 2020). While the X-ray emissions have been predicted by classical accretion disk models (Ulmer 1999), UV/optical emissions are argued to be produced from either the reprocessing of X-ray photons in an extended envelope or outflows (Loeb & Ulmer 1997;Strubbe & Quataert 2009;Lodato & Rossi 2011;Coughlin & Begelman 2014;Guillochon et al. 2014;Roth et al. 2016) or the shocks powered by debris stream self-intersection (Piran et al. 2015;Bonnerot et al. 2021). The late-time rebrightening of X-rays in TDEs can be accounted by either a change in the disk morphology as the accretion rate declines from super-Eddington to sub-Eddington, or the delayed onset of accretion Hayasaki & Jonker 2021). The latter model, however, is disfavored by recent simulations which show that a large fraction of the debris already settles into a disk with moderate eccentricity within dynamical times for at least a subset of TDE parameters (Bonnerot et al. 2021;Andalman et al. 2022).
In search for a unified model that can explain both the optical and X-ray TDEs, Dai et al. (2018) have carried out the first three-dimensional (3D) general relativistic radiation magnetohydrodynamics (GRRMHD) simulation tailored for TDE super-Eddignton accretion flow. The simulated disk, around a black hole with M BH = 5 × 10 6 M and spin parameter a = 0.8, has an accretion rate ofṀ acc ∼ 10Ṁ Edd , representing a typical accretion level in TDEs. The spectra of the disk have been obtained by post-processing the simulated disk with a novel radiative transfer code. It is found that the observed emission type largely depends on the viewing angle of the observer with respect to the disk orientation. When viewed face-on, X-ray emissions can escape from the optically-thin funnel surrounded by winds. When viewed edge-on, X-ray emissions are heavily reprocessed in the geometrically and optically thick wind and disk, so only UV/optical emissions can be observed.
While this study gives a good first-order description of TDE disks and spectra, recent simulations also show that the properties of a super-Eddington disk depends on the accretion rate, the black hole mass and spin, as well as the magnetic flux threading the disk (Jiang et al. 2014(Jiang et al. , 2019McKinney et al. 2015;Sadowski & Narayan 2016). Although the black hole spin and disk magnetic flux might only mildly affect the main structure of the disk, they determine whether a relativistic jet can be launched (Blandford & Znajek 1977;Curd & Narayan 2019). At first glance, the black hole mass is not expected to affect the disk spectra significantly, since most TDE host MBHs have masses in a narrow range ≈ 10 5 − 10 7 M . However, the peak fallback rate of the TDE debris sensitively depends on the black hole mass with the relationṀ fb,peak /Ṁ Edd ∝ M −3/2 BH . Therefore, TDEs from smaller MBHs should in general have much higher accretion rates at peak than those from larger MBHs (Ramirez-Ruiz & Rosswog 2009). Also, further variance of the accretion rates in different TDEs asṀ fb,peak can be introduced by the difference in the masses and ages of the disrupted stars .
In this Letter, we investigate how the TDE disk structure and the accompanied emission are influenced by the accretion rate at super-Eddington rates. We conduct three 3D GRRMHD simulations of super-Eddington disks with similar structures but varying accretion rates, and post-process the simulated disks to obtain their spectra at different inclination angles. These simulations are useful for understanding the diversity of the emissions observed from different TDEs, as well as the evolution of single TDEs as their fallback and accretion rate decline after peak. In Section 2 we introduce the setup of the GRRMHD and radiative transfer simulations. In Section 3 we give the main results and compare with TDE key observables. In Section 4 we draw a summary and discuss the implications and future work.

Disk simulation setup
We carry out three dimensional simulations of super-Eddington disks using a fully 3D GRRMHD code HARMRAD with M1 closure (McKinney et al. 2014). In all simulations the MBH has M BH = 10 6 M and a = 0.8. The radial grid has 128 cells spanning from R in ≈ 1.2R g to R out = 8500R g (R g = GM BH /c 2 is the gravitational radius of the MBH), with cell size increasing exponentially until R break ≈ 500R g and then increasing hyper-exponentially. The θ-grid has 64 cells spanning from 0 to π with finer resolution in the jet and disk regions. The φ-grid has 32 grids spanning uniformly from 0 to 2π with periodic boundary conditions. The gas is assumed to have solar chemical abundances (mass fractions of H, He, and "metals", respectively, X = 0.7,Y = 0.28, Z = 0.02). Frequency-mean absorption and emission opacities are used as in McKinney et al. (2015), except that the Chianti opacity is removed as it is unimportant for the disk temperature of our TDE models. Thermal Comptonization is also included.
We tailor the disk initial conditions to be consistent with realistic TDE scenarios following the setup described in Dai et al. (2018). The accretion disk is initialized with Keplerian velocity profile with a rest-mass density which is Gaussian in angle with a height-to-radius ratio of H/R ≈ 0.3, so that the initial density profile is given by where H is the scale height, ρ 0 is the initial reference density. The disks are checked to have total masses and specific angular momentum consistent to the first order with the conditions in TDEs. Adjusting ρ 0 leads to different accretion rates after the disk reaches the quasi-equilibrium state. AsṀ fb,peak can reach ≈ 100Ṁ Edd in a TDE around a 10 6 M black hole, and a large fraction of the debris mass is expected to escape in outflows, we set the aimed accretion rates to be a few to a few tens ofṀ Edd (Table 1). A large-scale poloidal magnetic field is initially seeded. As adopted in McKinney et al. (2015), for r smaller than a breaking radius R b = 500R g , the φ component of the vector potential is given by A φ = MAX(r 10 40 − 0.02, 0)(sin θ) 5 . For r ≥ R b , the field transitions to a split monopolar field, which is given by A φ = MAX(R b 10 40 − 0.02, 0)(sin θ) 1+4(R b /r) . The field is normalised so that the initial ratio of gas+radiation pressure to magnetic pressure β ≡ (p gas + p rad )/p b has a minimum value of 10.

Radiative transfer setup
We use the Monte Carlo radiative transfer code SEDONA (Kasen et al. 2006) to post-process the simulated disks and calculate the escaping radiation. More specifically, we calculate the bound electron level populations under non-local thermal equilibrium (LTE) as in Roth et al. (2016) and the Comptonization of electrons as in Roth & Kasen (2018). We track the free-free interactions, as well as the bound-bound and bound-free interactions with their opacities obtained from. the atomic database CMFGEN. The gas is assumed to consist of only H, He, and O with solar abundances. We focus on calculating the SED, and leave the line modeling to a future study which requires higher-solution simulations.
For each accretion rate, we calculate the spectra for four inclination bins of θ bin = 10, 30, 50 and 70 • . In each bin, we take an average over θ bin ± 5 • for the already time-and-φ−averaged profile of the simulated disk. Since the simulated jet density is likely numerically boosted, for the bins at θ bin = 10 and 30 • , we also reduce the jet density by 100 times before taking the average. We note that the first-order behavior of the spectra does not depend on this arbitrary choice of density re-scaling inside the jet. We set the source to be a blackbody spectrum with a single temperature of 10 6 K, which is consistent with the inner disk temperature. The source photons are injected from an inner boundary set at the boundary between the inflow and outflow, which is typically at a few R g except for the bin at i = 70 • which is partly in the disk inflow region. For the latter, we place the inner boundary at R = 20R g instead, and set the gas velocity to be always zero. The photons are then propagated in 3D, by assuming that the gas density, temperature, and radial velocity profiles are spherically symmetric. The photons propagate outwards until they leave the computational domain set at R = 4000R g or have become a part of the kinetic/thermal pool.
Based on the Monte Carlo calculations, we iterate the gas temperatures, gas ionization state, bound electron states, and radiative transfer solution under the assumption of radiative equilibrium until a steady solution has been reached (after approximately 20 iterations). Since GRRMHD simulations show that the luminosity of super-Eddington disks are not always capped by L Edd , we also tweak the source photon luminosity and obtain the spectra under two limits -the escaped bolometric luminosity is either L bol = L Edd ≈ 10 44 ergs −1 or L bol = 0.1Ṁ acc c 2 ≈ (Ṁ acc /Ṁ Edd ) L Edd .

Properties of the accretion flow
In all three GRRMHD simulations, it takes an initial time period of t = few × 1, 000R g /c for the accretion flow to get established. As the disk evolves, magnetic fluxes accumulate near the MBH and their strength further grow via the magnetorotational instability (MRI) (Balbus & Hawley 1998). Here we focus on the final stage of the simulations when the accretion flow has reached steady state and the wind has established equilibrium at most inclination angles. The black hole parameters and some basic quantities of the accretion disks, averaged over the period of t = 15, 000 − 20, 000R g /c, are listed in Table 1. The disk profiles used for post-processing are also first φ−averaged and then time-averaged over the same period. More numerical details of the simulated disks are given in Appendix A. The averaged accretion rates onto the event horizon of the MBH in these three simulations areṀ acc ≈ 7, 12 and 24Ṁ Edd . To first order, these super-Eddington accretion flows have similar structures as those described in Dai et al. (2018). As shown in the upper panels of Fig. 1, geometrically and optically thick disks are formed with half-angular thickness H/R ≈ 0.20 − 0.25, with, as expected, higher gas density towards the equatorial plane. The rotation profile of the resultant thick disks is sub-Keplerian. Wide-angle winds are launched by the large radiation and magnetic pressure, which are also optically thick at most angles except close to the polar axis. There are roughly two components of the wind: an ultra-fast outflow (UFO) with speeds faster than 0.1c within 45 • from the pole, and a slower but denser wind outside 45 • inclination. At larger accretion rates, the disk and wind densities increase, and the winds become slower as a result of mass-loading. However, the wind geometry stays fairly similar.
The emission properties of the accretion flow should be examined near the photosphere. The optical depth for an opacity κ is calculated using τ (r) = ρκdl along the radial direction, r, from the outer boundary R out = 8500R g towards the event horizon. (Here relativistic effects are included so that: dl = − f γ dr, with f γ ≈ u t (1 − v/c) and v/c ≈ 1 − 1/(u t ) 2 ). The electron scattering photosphere is defined by τ es (r) = 1 with κ = κ es ≈ 0.2(1 + X) cm 2 g −1 . The effective photosphere is then defined by τ eff (r) = 1 with κ = κ eff = √ 3κ ff (κ ff + κ es ), where we have only considered free-free opacity κ ff ≈ 3.82 × 10 22 Z(1 + X)(1 − Z)ρT −3.5 g in the scattering-dominated gas, where T g is the gas temperature. Both photospheres are plotted over the disk profiles in the lower panels of Fig. 1. As expected, we see the presence of an optically thin "funnel" near the rotation axis where the wind density is correspondingly lower. As the accretion rate increases, the wind becomes more opaque, which reduces the angular size of the funnel.
The effective photospheres reside mostly within r ∼ 5000R g except along the equitorial direction for the disk with the largesṫ M acc . Therefore, we select to evaluate various physical quantities at r = 5500R g . This gives averaged wind mass rates ofṀ w = 1.4, 4.5 and 14Ṁ Edd , and averaged bolometric luminosities at L RAD = 5.4, 3.3 and 8.1 L Edd , respectively, for the three simulations. The radiation temperature of the accretion flow is also plotted, which varies from ≈ 10 6 K near the black hole to ≈ 10 4 K in the outer part of the disk and wind. The radiation flux varies with the inclination as shown in Dai et al. (2018) and leaks out through the region of least resistance, which is the funnel.
Magnetic fluxes are dragged by the accretion flow and accumulate near the MBH. This brings the inner regions to swiftly become a magnetically arrested accretion disk (MAD) (Narayan et al. 2003). Relativistic jets are launched magnetically through the Blandford-Znajek mechanism (Blandford & Znajek 1977) in all simulations. In this Letter we focus on calculating the emission properties from the disk, and leave the analysis of the jets for future study. Figure 1. The 2D vertical profiles of time-and-φ-averaged gas rest-mass density ρ0 (upper panels, zoomed into the inner regions) and radiation temperature TRAD (lower panels, whole range of the simulation box), for the three runs with different accretion rates (from left to right: Macc = 7ṀEdd, 12ṀEdd and 24ṀEdd) in the quasi-equilibrium state. In the upper panels, we show the contours of constant lab-frame radial velocity (vr ≡ u r /u t ) with white lines, and mark the jet regions where the electromagnetic energy is larger than the rest-mass energy of the gas with dark blue color. In the lower panels, the black lines indicate the electron-scattering photosphere with τes = 1, and the red lines indicate the effective photosphere with τeff = 1. Larger accretion rates induce larger disk/wind density and higher gas/radiation temperature, while the gas distribution and velocity structure remain rather robust against the variance in accretion rates. The sizes of photospheres generally increase as accretion rate increases.

Spectra from post-processing
In this Section we investigate how the emission reprocessing depends on a few key parameters, namely, the viewing angle, the accretion rate, and the luminosity. The dependence on viewing angle has been previously studied by Dai et al. (2018). They show that for a fixed accretion rate, there is a clear trend of the spectral energy distribution (SED) moving towards larger wavelengths with increasing inclination angle. At low inclinations, the gas density is lower, so X-rays produced from the inner disk easily escape. At relatively large inclinations, on the other hand, the optically thick wind and outer disk serve as an effective reprocessing envelope. More specifically, in the fast wind region, the photons lose energy as they go through multiple scatterings in the expanding outflow before escaping. In the disk and the slow wind region, the reprocessing is mainly caused by the absorption of X-ray photons and the thermal Comptonization of electrons.
We re-examine the emission viewing-angle dependence in the three new simulated disks. Fig. 2a, 2b and 2c show how the escaped spectrum varies with inclination, whenṀ acc = 7, 12, and 24Ṁ Edd respectively. At any of these accretion rates, the SED still evolves from X-ray to optically dominated as the inclination increases. However, the exact angle at which this transition occurs depends on the accretion rate. One can clearly observe that at low accretion rates ≈ few ×Ṁ Edd , the escaped emission is dominated by X-ray emission at most inclinations unless the inclination angle is substantial (e.g., i 70 • forṀ acc = 7Ṁ Edd ), while at very high accretion rates ≈ few × 10Ṁ Edd , the emanating radiation is mainly optically dominated except near the polar region (e.g., i 10 • forṀ acc = 24Ṁ Edd ).
We expect a TDE to be observed along a fix viewing angle during its entire evolution, unless the disk experiences some axial precession. Therefore, we show the spectral evolution withṀ acc at fixed inclination angles in Fig. 2e-2h. As the accretion rate declines after peak, the amount of reprocessing material is reduced, so the SED universally shifts towards lower wavelengths. Figure 2. The simulated escaping spectra of the accretion disk at different accretion rates (Ṁacc = 7, 12, and 24ṀEdd) and inclinations (i = 10 • , 30 • , 50 • , 70 • ). The bolometric luminosity of the spectra Lbol = LEdd for all spectra. The purple shaded region indicates the X-ray band with an energy above 0.3 keV. The orange shaded region corresponds to Swift UVOT band at 1700-6500 Å. Panel (a)-(c) show the spectral evolution with inclination angles at fixed accretion rate. All spectra change from X-ray strong to UV/optical strong as the inclination goes from the polar direction to the disk direction. Panel (d) is the same as panel (b), but only includes the spectrum at i = 10 • with a blackbody spectrum fitting its X-ray continuum component, and the spectrum at 70 • with another blackbody spectrum fitting its UV/optical continuum component. Panel (e)-(h) show the spectral evolution with accretion rates at fixed inclinations. Three types of evolution can happen as accretion rate decreases: X-ray strong all the time (small inclination), optical/UV strong at early time and X-ray rebrightening at late time (intermediate inclination), optical/UV strong all the time (large inclination).
However, the exact behavior of the spectrum depends sensitively on the inclination angle. At very small inclinations (i = 10 • ), the TDE stays X-ray dominated throughout its evolution. At intermediate inclinations (i = 30 • and 50 • ), the TDE can be optically strong at early times (ifṀ acc can be sufficiently large), and then rebrightens at X-ray energies at late time when the accretion rate diminishes. At very large inclinations (i = 70 • ), the TDE stays UV/optically strong throughout its entire evolution.
While we have assumed that the bolometric luminosities of the escaped radiation always close to L Edd in the analysis above, GRRMHD simulations show that the true escaped luminosity from super-Eddington disks can exceed the Eddington limit, with more flux leaking out through the funnel (e.g., Sadowski & Narayan 2016). Therefore, we also calculated the spectra when the escaped radiation has a bolometric luminosity L bol = 10%Ṁ acc c 2 , which are shown in Fig. B.1. When the reprocessing is driven by adiabatic expansion, the SED shape stays unchanged while the magnitude of the flux scales with the luminosity. When the reprocessing mechanism is driven by absorption and Comptonization, increasing the luminosity makes the gas more ionized and reduces bound-free and free-free absorption, which shifts the spectral energy peak to slightly higher energies. However, the spectral shape is rather insensitive to the luminosity. In general, the setting of the luminosity within the explored range does not alter how the escaped radiation depends on the viewing angle and accretion rate.

Comparison with observations: blackbody luminosity, effective temperature, and photosphere radius
In this section, we compare our model predictions to the observed properties of TDEs. We start from the TDE catalog in Gezari (2021) which also lists the observed blackbody luminosities and temperatures, and then only include TDEs which have their masses estimated from the M − σ relation as in Wong et al. (2022). This gives us 7 X-ray-selected TDEs (Table C.1) and 16 optically-selected TDEs (Table C.2). We plot their observed blackbody luminosities, effective temperatures and photosphere radii as functions of M BH in Fig. 3a -3c. The observed L BB varies between 10 −3 − 10 L Edd , and can exceed L Edd when M BH 10 6 M . Interestingly, L BB commonly has a dependence on M BH following the fallback rate trend:Ṁ fb /Ṁ Edd ∝ M −3/2 BH . The observed T BB clearly depends on whether the TDE is optically-selected (≈ few × 10 4 K) or X-ray-selected (≈ 10 5 − 10 6 K). The observed R BB , calculated from L BB and T BB using the Stefan-Boltzmann law, also has a bimodal distribution. The optically-selected TDEs have R BB reaching 10 2 − 10 4 R g and always exceeding both the circularization radius R circ = 2R t and the stream self-intersection radius R int . The X-ray-selected TDEs, however, sometimes have R BB even within the event horizon, which might be caused by absorption in circumnuclear medium.
As the observed properties of TDEs are usually inferred from their spectra in monochromatic bands, we also fit the simulated spectra in either the X-ray or UV/optical band with blackbody spectra. For example, in Fig. 2d, we show the simulated escaped spectra forṀ acc = 12Ṁ Edd at i = 10 • and 70 • and the blackbody radiation fits to the spectra in the X-ray band (0.3-10 keV) or the part of the spectrum in the UV/optical band (the 170-650 nm Swift UVOT band). It can be seen that the X-ray and UV/optical emission can be individually well approximated by blackbody radiation. Similarly, for each spectrum (atṀ acc = 7, 12, or 24Ṁ Edd and inclination i = 10 • , 30 • , 50 • , or 70 • ), we obtain two blackbody radiation fits in the X-ray band and the UV/optical band. The escaped bolometric luminosity is assumed to be varying between L Edd and 0.1%Ṁ acc c 2 . The luminosities, temperatures and radii of the best-fit blackbody radiation spectra are are listed in Table D.1 and plotted in Fig. 3d -3f (UV/optical fit) and Fig. 3g -3i (X-ray fit). We further categorize whether a modeled TDE spectrum is a X-ray strong or optically strong by comparing the luminosity in the X-ray band (0.3-10 keV) and the blackbody luminosity inferred from the UV/optical band. We compare the model predictions to the observed properties as below: 1. Luminosity: The modeled L BB mostly lies between 0.01-few L Edd . The simulated spectrum usually peaks in EUV and has a broader shape than a single-temperature blackbody spectrum. Therefore, the inferred blackbody luminosity (L O,BB or L X,BB ) is always smaller than the bolometric luminosity L bol of the escaped radiation. L O,BB /L bol for optically-strong TDEs and L X,BB /L bol for X-ray-strong TDEs is typically few×(1 − 10)%. This naturally explains the missing energy problem in TDE ; Lu & Kumar 2018; Mockler & Ramirez-Ruiz 2021) -the majority of the energy is emitted in the EUV band which is difficult to be detected.
2. Temperature: Our modeling reproduces the bimodal distribution of TDE temperatures, i.e., optically-strong TDEs have temperature T O,BB ≈ few × 10 4 K, and X-ray-strong TDEs have temperature T X,BB 10 6 K. Furthermore, our modeling shows that the inferred temperatures of optical TDEs are not highly sensitive to accretion rates or observer inclination, which explains why TDEs have relatively constant T BB throughout the evolution (van Velzen et al. 2020).
3. Photosphere radius: The optically-strong TDEs have far larger photosphere radii (R O,BB ≈ 10 2 − 10 4 R g ) than X-raystrong TDEs (R X,BB ≈ few × R g ). A comparison between the observed R O,BB = 10 3 − 10 4 R g and the modeled R O,BB suggests that optically-selected TDEs are either commonly observed from large inclinations, or have L bol > L Edd . The distinction between the observed and modeled R X,BB is possibly accounted by the absorption of X-rays in the circumnuclear medium. TDEs have RBB larger than the circularization radius (red curve) or the stream self-intersection radius (green curve) (both calculated using a 0.1M star). X-ray TDEs can sometimes have RBB smaller than the black hole Schwarzschild radius (black line). Middle row (d)-(f) and lower row (g)-(i): The inferred quantities based on the blackbody radiation spectrum fitting the simulated spectra in the UV/Optical band or X-ray band respectively vs. inclination angle i. Different symbols are used to mark different accretion rates: 7ṀEdd (blue circle), 12ṀEdd (green triangle) and 24ṀEdd (red square). Vertical lines connect the values calculated with an escaped luminosity of Lbol = LEdd (smaller symbol size) and Lbol = 10%Ṁaccc 2 (larger symbol size) to indicate possible ranges. In (d)-(f), TDEs with LO,BB < LX,0.3−10keV are masked with low opacity to indicate they are less likely to be selected optically. Similarly, in (g)-(i), TDEs with LO,BB > LX,0.3−10keV are masked with low opacity to indicate they are less likely to be selected by X-ray instruments. The blackbody luminosity, temperature and radius inferred from our modeling, to the first order, reproduce the observed ones.
We note that our predictions for the X-ray quantities, in particularly in the i = 10 • bin, are sensitive to the setting of the radiative transfer calculations. Here we always inject a blackbody spectrum with a constant T = 10 6 K, so our predicted T X,BB at small inclinations also fall into a very narrow range. However, the temperatures at the center of accretion disks generally increase with increasingṀ acc and decreasing M BH , which will induce more variance to the observed X-ray temperatures. Also, for the setting with L = L Edd , the predicted L X,BB decreases asṀ acc increases, as a result of a constant bolometric luminosity and more reprocessing from X-ray to UV/optical emissions at higher densities. However, simulations show that the radiation fluxes leaking out through the funnel are not Eddington-limited (McKinney et al. 2015;Dai et al. 2018). Therefore, we expect that X-ray luminosities should scale positively with accretion rates, as illustrated by the L = 0.1Ṁ acc c 2 case.

Temporal evolution of TDEs
We show in Fig. 4 the evolution of the modeled luminosity, temperature, radius, as well as the ratio of optical to X-ray luminosity, as functions of the accretion rate. In order to connect the snapshots at an specific accretion rate to the temporal evolution of TDEs, we assumeṀ fb (t) =Ṁ acc (t) +Ṁ wind (t). This assumption is valid only if the fallback timescale dominates over other timescales, such as the disk viscous timescale and the photon diffusion/advection timescales. The exact conversion fromṀ fb to t depends on the mass of the MBH, the properties of the disrupted star, and the impact parameter . Focusing on the post-peak evolution, the three disk simulations correspond to 45.9, 102.4 and 174 days after the peak of the flare, assuming a solar-type star is fully disrupted by a 10 6 M black hole. Calculations of various timescales are given in Appendix E.
As the accretion rate drops from 24Ṁ Edd to 7Ṁ Edd , the optical luminosity also decreases, and the UV/optical light curve roughly follows the canonical t −5/3 decay. As discussed in the previous section, the fitted blackbody temperature stays rather constant. Interestingly, at large inclinations T O,BB slightly decreases with decliningṀ acc , while at small to intermediate inclinations T O,BB shows the opposite trend. This can provide an explanation to the different observed evolution of T O,BB (van Velzen et al. 2020). As a result, the photosphere radius shrinks asṀ acc decreases, with a faster evolution at smaller inclinations.
The ratio between UV/optical luminosity (L O,BB ) and X-ray luminosity (L X,0.3−2keV ) also decreases as the accretion level drops and the amount of obscuring material is reduced. The fastest L O,BB /L X,0.3−2keV evolution is observed at intermediate inclinations, and especially if the TDE has a high accretion rate at peak. In such cases, we expect to see strong X-ray re-brightening of initially optically-strong TDEs, as reported in Gezari et al. (2017). The X-ray luminosity reaches the same level as the UV/optical luminosity at t 100 days after peak, and possibly even later if the disk formation or viscous timescale is long. At very small inclinations, the event is always X-ray strong. At very large inclinations, the event should stay optically strong for a long period, although it is theoretically predicted that the disk should eventually become geometrically thin and emit mostly in X-rays/UV when the accretion level drops to around Eddington (Shakura & Sunyaev 1973).

SUMMARY AND FUTURE WORK
Inspired by the unified model for TDEs proposed by in Dai et al. (2018), we carry out three addtional 3D GRRMHD simulations of TDE super-Eddington accretion flow at different Eddington ratios, and conduct radiative transfer calculations to obtain the emanating spectra. Based on the results, we further propose a dynamical unified model which can explain the diversity and evolution of TDEs: • The viewing angle of the observer with respect to the orientation of the disk is the most important parameter in determining whether we observe either an X-ray or an optical bright TDE. At small inclinations, X-rays can escape from the funnel of the super-Eddington disk. At large inclinations, X-rays are mostly reprocessed into UV/optical radiation by the geometrically and optically thick wind and disk.
• The blackbody radiation fits of the TDE super-Eddington disk spectra produce effective temperature, blackbody luminosity and photosphere radius distributions consistent with the observed ones. Most radiative energy escapes in the EUV range, and only a few to a few tens of percentage of radiative energy can be detected, which provides a solution to the TDE missing energy problem.
• The observed diversity of the emission from different TDEs can be associated with the different Eddington ratios of their accretion rates,Ṁ acc = few × (1 − 10)Ṁ Edd , at the flare peak conditions. In general, higher accretion levels induce larger (fitted blackbody) luminosities and larger photosphere radii, but do not significantly change the fitted effective temperatures.
• The early-time evolution (t ≈ 100 days after peak) of optical TDEs can be explained by this reprocessing model. As the luminosity drops by about 0.5-1 order of magnitude, the fitted temperature slightly increases at small to intermediate inclinations or decreases at large inclinations. , radius (c), and the ratio between the UV/optical and X-ray luminosity (d). Different colors denote different inclination angles. The escaped radiation has luminosity Lbol = LEdd for all curves. The lower axis shows the accretion rate and the upper axis shows the corresponding time elapsed since the peak, assuming a solar-type star disrupted by a 10 6 M black hole. In (a)-(c), we do not include the evolution i = 10 • , where the event is always X-ray strong. In (a) the gray line shows the trend of t −5/3 to guide the eye. In (d) the X-ray luminosity includes only the flux in the 0.3-2 keV band for direct comparison with observations. Also at i = 70 • the X-ray luminosities at the two higher accretion rates are negligible.
• The evolution of the optical-to-X-ray flux ratio also depends sensitively on the viewing angle. At large inclinations, the TDE stays UV/optically strong for a very long time. At intermediate inclinations, we expect to see the fastest X-ray rebrightening, and L O /L X reaches unity at few hundred days. At small inclinations, the TDE is always X-ray strong. The exact evolution timescale also depends on the accretion rate at peak, which further depends on the black hole and stellar parameters.
In this study we focus on understanding the disk continuum emission. We will investigate the disk spectroscopic features in a future work. Also, while our GRRMHD simulations are conducted in 3D, we only do post-processing along different inclination bins by assuming the profiles are spherically symmetric, which means that photons have an overall radial motion. Implementing 2D or 3D radiative transfer calculations can allow us to study the reprocessing of photons along more realistic paths. Furthermore, the evolution of the disk and jet emissions can be studied in more details by conducting simulations covering a more complete parameter space for black hole accretion rate and spin.
The unified model we have proposed for TDEs also sheds light on other types of black hole accretion systems, especially other super-Eddington sources such as ultra-luminous X-ray sources, narrow line Seyfert 1 galaxies, changing-look AGN and quasars in the early universe. Works like this illustrate that novel simulations are crucial for making direct comparisons between the observed emissions and model predictions, which offers important insights to how black holes grow through cosmic time, use accretion as a source of energy to power outflows and radiation, and give feedback to their surrounding environment.

APPENDIX A. PROPERTIES OF SIMULATED DISKS
A few more time-averaged quantities of the three simulated disks are listed in Table A.1. Φ H is the normalised magnetic flux at horizon Φ H (Tchekhovskoy et al. 2011), with Φ H 30 − 40 being the condition for MAD disks. The net accretion efficiency η H evaluates how much rest-mass energy going into the black hole is converted to out-going energy near the event horizon. L jet is the total power of the relativistic jet at r = 5500r g , most of which is in the form of electromagnetic energy. Γ jet,max is the maximum Lorentz factor of the jet. L K is the thermal+kinetic luminosity of the wind calculated at r ∼ 5500r g . α eff is the effective α-parameter of the disk as defined in McKinney et al. (2012). R eq indicates the radius within which the disk inflow equilibrium has been established. In our simulations, the winds are launched mostly from the inner disk regions and have traveled beyond the photospheres at most inclination angles.

C. LIST OF OBSERVED TDES REPORTED IN LITERATURE
For plotting Fig. 3a-3c we use 16 optically-selected TDEs and 7 X-ray-selected TDEs. For completeness we list their names and relevant parameters as reported from pervious literature in Table C.1 and C.2. The photon transport time through optically thick medium is the shorter one between the diffusion timescale and advection timescale. The photon diffusion timescale is calculated as t diff = τ es R es /c, where R es is the size of the electron-scattering photosphere along a particular inclination, and τ es is the electron-scattering optical depth integrated radially from r = 0 to r = R es .  The advection timescales is calculated as t adv ≈ R es /v r , since the photons trapped by the optically-thick gas has a similar speed as the gas. Here the gas radial velocity v r is averaged over the radial path within R es and weighted by gas density. The values for theṀ acc = 12Ṁ Edd disk are given in Table E.1. One can see that for the inclinations considered in this work, photons are preferably advected out by the optically-thick wind. The photon transport time varies from ∼0.1 day to a few days depending on the inclination. The disk viscous timescale can be analytically calculated by t visc = t dyn α −1 (H/R) −2 , where t dyn is the orbital timescale of the disk, α is a free parameter between 0 and approximately 1, and H/R is the disk thickness (Shakura & Sunyaev 1973). For our disk parameters, we have t visc ≈ 5.44 days (R disk /8500r g ) (M BH /10 6 M )(α/1) −1 ((H/R)/0.3) −2 . The viscous time is therefore only a few days. We caution the readers that our simulated disks are MAD which typically have effective α 1. The viscous time is potentially longer if the disks do not have such large magnetic fluxes.
The disk formation/circularization timescale induces uncertainty into the evolution of TDEs. Recent simulations show that a large fraction of the debris materials can form a disk within dynamical timescale but the disk still pocesses some moderate eccentricity (Bonnerot et al. 2021). As the topic is out of the scope of this paper, we assume here that the disk forms quickly for the calculation of the emission evolution.
Stellar debris typically falls back on timescales of t mb , which is the orbital time of the most bound debris (Evans & Kochanek 1989;Guillochon & Ramirez-Ruiz 2013;Rossi et al. 2021 Lastly, assuming that the fallback timescale governs the evolution of TDEs, we calculate the time corresponding to each simulated disk with a further assumption that the fallback rateṀ fb equals the instataneous accretion rateṀ acc plus the wind mass rateṀ w . Here we focus on the post-peak evolution. By setting that the fallback rate peaks at t mb , the post-peak time t pp associated with a particularṀ fb can be calculated with: With the wind mass rate given in Table 1, the three simulated disks correspond toṀ fb =Ṁ acc +Ṁ w = (7+1.4, 12+4.5, 24+14)Ṁ Edd . The post-peak time obtained with M BH = 10 6 M and a few different m are shown in Table E.2.