Accretion Disc Winds in Tidal Disruption Events:
Ultraviolet Spectral Lines as Orientation Indicators

All calculations for this work were completed using Python¹¹1Python is a collaborative open-source project. The source code is available at https://github.com/agnwinds/python., a state-of-the-art Monte Carlo ionization and radiative transfer code for moving media that uses the Sobolev approximation (e.g. Sobolev, 1957; Rybicki & Hummer, 1978). Originally described by Long & Knigge (2002), Python’s evolving structure and capabilities have been summarised several times in the literature (Higginbottom et al., 2013, 2014; Matthews et al., 2015, 2016). We therefore only provide a brief overview.

Python operates in two stages. The first determines the ionization state, level populations and temperature structure of the modelled outflow. This is done by iteratively flying several populations of Monte Carlo energy quanta (“photon packets”) through a spatially discretised grid. A number of different coordinate systems are supported. The most commonly used are a 2.5D logarithmic polar or logarithmic cylindrical grid. The outflow is discretised over $n_{1}\times n_{2}$ grid cells in these coordinates. Photon packets are generated over a wide wavelength range, sampling the spectral energy distribution (SED) of each radiation source in the model. As these packets interact with the plasma, Monte Carlo estimators that describe the radiation field in each cell are updated. The heating effect of the photon packets is also logged and used to move the temperature towards thermal equilibrium, where heating and cooling are balanced in each cell²²2Python includes all significant radiative heating and cooling processes, along with adiabatic cooling due to the expansion of the outflow. Compton heating and inverse Compton cooling are also included, together with Auger ionization and dielectronic recombination (see Higginbottom et al., 2013).. After each iteration, the revised electron temperature and radiation field estimators are used to update the ionization state of the model. This process is then repeated until the model converges, with heating and cooling balanced throughout the wind.

The second stage of the code produces synthetic spectra. Additional photon packet populations (typically over a narrow wavelength range to give high signal-to-noise) are flown through the converged model to generate spectra for a selection of sight lines.

Throughout this work, we examine two distinct wind geometries. Both are based on a kinematic biconical outflow prescription originally designed by Shlosman & Vitello (1993) and illustrated in Figure 2. In this prescription, wind streamlines emerge from disc radii between $r_{\text{min}}$ and $r_{\text{max}}$ at angles relative to the disc normal given by

where,

Here, $r_{0}$ is the launch radius of the streamline, $\theta_{\text{min}}$ and $\theta_{\text{max}}$ are the minimum and maximum opening angles of the wind, and $\gamma$ controls the concentration of streamlines towards either of the two boundaries $r_{\text{min}}$ and $r_{\text{max}}$.

The launch velocity, $v_{0}$, and terminal velocity, $v_{\infty}$, of a streamline are set to the local sound speed and a multiple of the escape velocity at the streamline footpoint, respectively. The poloidal velocity, $v_{l}$, at a distance $l$ along a streamline, is defined by

where $R_{v}$ is the acceleration length scale, and $\alpha$ controls how rapidly the wind accelerates. The rotational velocity, $v_{\phi}$, is Keplerian at the footpoint of a streamline and assumed to conserve specific angular momentum as it rises above the disc. At a cylindrical radius $r$, $v_{\phi}$ is then given by

where $v_{k}$ is the Keplerian velocity at $r_{0}$. The density at any point in the wind is obtained by requiring mass continuity, giving

In this equation, $\dot{m}^{\prime}$ is the mass-loss rate per unit surface area of the accretion disc,

where $\lambda$ is a parameter that controls where mass is lost from the disc. A value of $\lambda=-2$ results in mass being lost roughly uniformly across the accretion disc surface.

We use two sets of model parameters to create (i) a wide-angle and (ii) an equatorial disc wind, inspired by previous efforts to model the outflows from accreting WDs and QSOs, respectively. The difference between these two geometries is illustrated in Figure 3.

We assume a black hole mass of $3\times 10^{7}\text{M}_{\odot}$ and model an accretion disc using parameters inspired by the characteristic (blackbody) luminosity inferred for iPTF15af (Blagorodnova et al., 2019) to estimate a reasonable accretion rate. More specifically, we assume that the blackbody luminosity of iPTF15af, $L_{\text{BB}}\sim 10^{43}\text{ergs s}^{-1}$, roughly corresponds to the luminosity of the accretion disc, $L_{\text{acc}}=\eta\dot{\text{M}}_{\text{disc}}c^{2}$, for an accretion efficiency of $\eta=0.1$. We thus calculate an estimate of the accretion rate to be roughly $\dot{\text{M}}_{\text{disc}}\simeq 2\times 10^{-2}\text{M}_{\odot}\text{yr}^{-1}$, corresponding to an Eddington accretion limit of $0.02~\dot{\text{M}}_{\text{Edd}}$ or $0.01~L_{\text{Edd}}$. The outer radius of the disc for both models is set to be $R_{\text{disc}}=10^{15}$ cm to roughly match the continuum level of iPTF15af. As shown later in Section 3.1, BALs are produced by both models, but for different viewing angles. We therefore adjusted the accretion rates of each model so that the predicted model flux for a BAL-forming sight line roughly matches the flux level of iPTF15af at $d=350$ Mpc. For the equatorial model, we found that the previously calculated accretion rate already matched well with observations. However, the wide-angle model required a slightly lower mass accretion rate of $\dot{\text{M}}_{\text{disc}}=10^{-2}\text{M}_{\odot}\text{yr}^{-1}$. This difference is mainly due to foreshortening: BALs are observed for more face-on inclinations in the wide-angle model.

For both models, we initially assumed a mass-loss rate of $\dot{\text{M}}_{\text{wind}}=0.1\,\dot{\text{M}}_{\text{disc}}$; based on previous experience of modelling disc winds in accreting WD and QSO systems. We conducted an initial coarse parameter search, focusing on parameters which can impact line formation, i.e. the velocity law and mass-loss rate of the wind, consequently finding that a larger mass-loss rate of $\dot{\text{M}}_{\text{wind}}=\dot{\text{M}}_{\text{disc}}$ resulted in stronger line formation for the wide-angle model. Our final parameter choices for the two wind models are presented in Table 1. We conducted a finer parameter search around these values in Table 1, finding that the ionization state and emergent spectra are fairly insensitive to moderate changes in these values.

In both models, the disc wind emanates from the entire surface of the accretion disc. We choose values of $\gamma=1$ for uniform streamline spacing and $\lambda=0$ corresponding to uniform mass loss across the disc surface. We use a logarithmic cylindrical grid with a greater concentration of cells at smaller radius, where most of the line formation is expected to occur. Several resolution tests were conducted. For grids with a spatial resolution lower than $50\times 50$ cells, we found that the line-forming region was inadequately spatially resolved. For resolutions higher than $75\times 75$ cells, we found only small changes to the emergent spectra. We use $100\times 100$ cells for our final calculations, ensuring we adequately resolve the line-forming regions. We adopt a maximum wind radius of $R_{\text{wind}}=5\times 10^{17}$ cm to ensure our computational domain is large enough to not to affect the results. We ran several models to test this, varying the maximum wind radius by factors of a few, finding little effect on the emergent spectra.

We now describe a clumpy wide-angle disc wind with solar abundances and clumping factor $f_{v}=0.1$ in the context of the formation of BALs and BELs. Based on our limited exploration of parameter space (see Sections 3.2, 3.3, 3.4 below), we use this as our fiducial disc wind model.

In Figure 4, we present a selection of physical parameters for the wide-angle wind on a log-log distance scaling. In Figure 5, we show the synthetic UV spectra produced by this model compared to the spectra of a BAL TDE (iPTF15af) and a BEL TDE (ASASSN14li). For our BAL model, we adopt an inclination of $i=60^{\circ}$, i.e. a sight line which looks into the wind cone. To represent BEL models, we show spectra for both face-on (i = $10^{\circ}$) and edge-on ($i=75^{\circ}$) sight lines, both of which lie outside the wind cone.

The spectra produced by both wide-angle and equatorial wind models are shown in Figure 6 (for smooth vs clumpy outflows) and Figure 7 (for solar vs CNO-processed abundances). The dependence of the spectra on clumpiness and abundances will be discussed in Sections 3.3 and 3.4, respectively. Here, we briefly comment first on the basic differences between the spectra produced by the two types of outflow geometry. The implications of these differences are discussed in more detail in Section 4.1.

In principle, both wide-angle and equatorial winds are able to produce BALs and BELs. In both geometries, BALs are observed preferentially for sight lines looking into the wind cone. However, wide-angle flows produce BELs for both low and high inclinations, whereas our equatorial model only produces BELs at low inclinations. Perhaps more importantly, in terms of line-to-continuum contrast, the equatorial model only produces very weak emission lines. This is partly due to the anisotropic radiation pattern produced by the disc in our model. As a result, the continuum level is much higher for the face-on orientations that can produce BELs in the equatorial wind model. The line emission produced by the wind is always more isotropic than the disc continuum, so it is harder to achieve high line-to-continuum contrasts in equatorial models. As discussed further in Section 4, it is actually not clear if real TDE discs emit this anisotropically. For a more isotropic continuum source, the emission line-to-continuum contrast will be considerably more uniform across all orientations. However, absorption features would still be formed preferentially for sight lines looking into the wind cone.

Figure 6 illustrates the effect of clumping on the emergent spectra. More specifically, we compare here the spectra produced by smooth winds ($f_{v}=1$) to those produced by clumpy winds with a clump to inter-clump density ratio of 10 (i.e. $f_{v}=0.1$).

The main effect of clumping is to strengthen the UV resonance transitions. This happens because clumping increases the density of the line-forming gas and hence lowers its ionization state. In the wind models presented here, which tend to be somewhat overionized, this increases the abundance of the ionic species responsible for the UV lines. For the “in-wind” (BAL) spectra, this manifests as deeper and broader absorption troughs. The BALs are also more highly blue-shifted, as the abundance of the relevant ions in the high velocity regions of the wind is increased. Emission lines are additionally strengthened in clumpy winds, because collisional excitation scales as the square of the density in the line-forming region.

The presence of clumps is particularly important for the formation of Si iv. As noted before, Si v tends to be the dominant ionization stage of Silicon in our models, but clumpy winds contain enough Si iv to produce BAL features for sight lines looking into the outflow. This feature is not generally observed in our smooth wind models.

Some UV spectra of TDEs present strong Nitrogen emission lines, but weak or undetectable Carbon lines; e.g. iPTF15af and ASASSN14li. This has been attributed to an enhanced Nitrogen to Carbon ratio in the line-forming region, as expected from TDEs associated with stars that have undergone CNO-processing. (Kochanek, 2016; Yang et al., 2017). We have therefore also carried out some simulations with abundances representing CNO-processed material. In these, we changed the relative abundances of Helium, Carbon and Nitrogen to $(X/X_{\odot})_{{}^{2}He}\approx 2$, $(X/X_{\odot})_{{}^{6}C}\approx 0.5$ and $(X/X_{\odot})_{{}^{7}N}\approx 7$, where $X_{\odot}$ is the relevant solar abundance. These values are based on abundances calculated by Gallegos-Garcia et al. (2018) for a $\simeq 1.8M_{\odot}$ star near the terminal main sequence. In Figures 7 and 8, we explore the impact of abundance variations by showing synthetic spectra produced by smooth and clumpy winds with solar and CNO-processed abundances, respectively.

For smooth equatorial winds, the change to CNO-processed abundances results in a slight suppression of the C iv line. The temperature and ionization structure of the wind is roughly the same for both abundance patterns, suggesting that this suppression is due to the change in the abundance of Carbon atoms, rather than a change in the ionization state. This suppression is less obvious in the wide-angle wind model, probably because C iv is more dominant, resulting in the line being more optically thick here. In any case, the effect of the CNO-processed abundances on the N v BAL is much greater. In both wind geometries, this feature becomes both deeper and broader as the density of Nitrogen atoms increases in the wind.

The effect of CNO-processed abundances on clumpy winds is similar, but stronger. Figure 8 shows a comparison between the two abundances for a clumping factor $f_{v}=0.1$. For the equatorial wind, we again see that the C iv absorption has weakened, whereas it remains unchanged for the wide-angle wind. As before, the strongest effect is on the N v BAL, which is once again deeper and broader with CNO-processed abundances.

BALs have been now detected in three out of four UV spectra of TDEs. As pointed out by Hung et al. (2019), this suggests that the appearance of these outflows is less sensitive to viewing angle than that of the outflows in QSOs, where BALs are seen in only $\simeq 20$ per cent of sources (e.g. Knigge et al., 2008; Allen et al., 2010). If BALs are observed preferentially for sight lines to the central engine that pass through the outflow, the fraction of systems displaying them is a measure of the wind "covering factor", i.e. the solid angle subtended by the wind, as seen by the central engine. In practice, selection effects will complicate this picture. For example, if bright systems are over-represented in observational samples, fore-shortening, limb-darkening and attenuation by the wind itself will all lead to orientation-dependent incompleteness. Moreover, TDEs are evolving systems, so their outflows and associated BALs are bound to be transient features. Nevertheless, for TDEs observed at similar stages of their eruptions, as is the case for the objects shown in Figure 1, the incidence rate of BALs should provide a reasonable initial estimate of the outflow covering fraction.

Our finding that a wide-angle wind model seems to match the observations better than an equatorial one is in line with this picture. The covering factor, $f_{\Omega}=\int_{\theta_{min}}^{\theta_{max}}\sin(\theta)~\text{d}\theta$, of our wide-angle wind model is $f_{\Omega}=0.52$, whereas that of our QSO-inspired equatorial model is $f_{\Omega}=0.20$. A wide-angle wind model spanning an even wider range of opening angles is clearly worth exploring in future.

In the context of our particular wide-angle wind model, where $\theta_{min}>0^{\circ}$ and $\theta_{max}<90^{\circ}$, BELs are seen for both very low and very high inclinations. Do TDEs displaying BELs correspond to extremely face-on systems, extremely edge-on ones, or both? For our fiducial model, the narrower line widths seen for face-on orientations appear to be more in line with the very limited observational data (c.f. Figure 5). On the other hand, extremely face-on orientations are a priori unlikely. Only 6 per cent of randomly oriented systems would be viewed at inclinations $i<20^{\circ}$, whereas 42 per cent would be viewed at $i>65^{\circ}$. In any case, the detailed line shapes produced by our models should not be over-interpreted: the parameter space of even just wide-angle models is large, and, as discussed in Section 4.3 below, our calculations still have significant limitations that are likely to affect the line profiles. As the sample of TDEs with UV spectroscopy grows, it will be interesting to check for correlations between the incidence of BALs/BELs and other system properties that might depend on orientation (such as the X-ray to optical ratio; Dai et al., 2018).

Several lines of evidence suggest that the outflows of TDEs may be optically thick. For example, the likely presence of P v $\lambda 1118$³³3Note that this transition is currently not included in Python’s line list. in iPTF15af suggests a column density of $N_{H}>10^{23}$ cm^-2 (Blagorodnova et al., 2019). A similar constraint has been obtained for AT2018zr from the detection of BAL signatures in Balmer and metastable Helium lines (Hung et al., 2019). From a modelling perspective, Roth & Kasen (2018) have shown that the detailed line profile shapes seen in TDEs can be understood if these features are formed in an optically thick outflow. Finally, the reprocessing of radiation produced by the central engine in an optically thick, non-spherical wind is also central to the unification scenario presented by Dai et al. (2018), which is based on 3-D general relativistic magneto-hydrodynamic simulations of the super-Eddington accretion phase near the peak of TDE flares.

How optically thick is our fiducial clumpy wide-angle wind model? Figure 9 shows the continuum (free-free, bound-free and electron scattering) optical depth through this model, as a function of frequency, for several lines of sight. Starting from a point at the inner disc edge, each sight line here runs radially outward, making an angle $i$ with the normal to the disc plane. At low frequencies, electron scattering is the dominant source of opacity. The electron scattering optical depth is $\tau_{es}>1$ for all of the sight lines shown ($30^{\circ}\leq i\leq 85^{\circ}$). For $i=30^{\circ}$, the optical depth is moderate, $\tau_{es}\simeq 3$, but it increases monotonically with inclination, reaching $\tau_{es}\simeq 30$ for $i=85^{\circ}$. This increase continues even for sight lines that exit the wind cone, i.e. $i>\theta_{max}$, because the shorter path through the wind for these sight lines is more than compensated by the higher densities found at the base of the wind. On average, each photon in our fiducial wide-angle wind models undergoes $\simeq\!\!100$ electron scatters, with a few photons scattering upwards of 2000 times.

At higher frequencies, the opacities are dominated by photo-ionization. In particular, the He ii edge at 54.4 eV (228 Å) is highly optically thick for essentially all inclinations. The Lyman edge at 13.6 eV (912 Å) and the He i edge at 24.6 eV (504 Å) only start to make a noticeable contribution to the opacity for high-inclination lines of sight that pass through the dense, less ionized wind base. Figure 9 also includes the angle-averaged input disc SED and the angle-averaged emergent flux of the model for comparison. This shows that the contribution of He ii ionizing photons, which will be almost completely reprocessed by the wind, is not dominant, but still significant in our fiducial model. The majority of the reprocessed emission in the model is emitted at optical wavelengths and at high inclination angles. The column densities along the sight lines shown in Figure 9 range from $N_{H}\simeq 5\times 10^{25}$ cm^-2 (at $i=30^{\circ}$) to $N_{H}\simeq 5\times 10^{26}$ cm^-2 (at $i=85^{\circ}$). These values are compatible with the lower limits suggested by Blagorodnova et al. (2019) and Hung et al. (2019).

As is the case for any computational model, our calculations have several limitations. First, Python is a time-independent code, so we implicitly assume that the flow is in a steady-state on time-scales short to the time-scales over which the TDE luminosity and SED evolve significantly. Perhaps most importantly, however, Python is currently a purely Newtonian code. In particular, this means we neglect general relativistic effects on the photon paths in our simulations. Whilst the full random walk of photons is simulated, photons travel in straight lines between interaction points throughout the computational domain. There is also no concept of an event horizon. As a simple test of this approximation, we repeated the simulation for our fiducial wide-angle wind, but this time removed all photons which intersected a sphere of radius $R_{\text{ISCO}}$, centred on the BH. We found that the ionization state and emergent spectra changed only marginally relative to the fiducial model. We also do not take into account the full detailed effects of special relativity, which will have an effect on the Doppler shifts experienced by photons and hence on the emergent line profiles. However, we do not expect this to affect the overall ionization state of the flow, nor the qualitative appearance of the emission and absorption line spectra. Nevertheless, the implementation of relativistic effects into Python is a high priority in our continuing development of the code and future work.

Another clear limitation is our description of the accretion disc as a standard steady state $\alpha$-disc (Shakura & Sunyaev, 1973). In reality, the inner discs in these systems are dominated by radiation pressure and are almost certainly vertically extended. However, the structure, evolution, and even the stability, of such discs are the subject of intense research in radiation-dominated regimes (e.g. Hirose et al., 2009; Jiang et al., 2013; Blaes, 2014; Shen & Matzner, 2014). Additionally, relativistic effects of accretion onto a Kerr black hole is also likely to play a crucial role in TDEs (Balbus & Mummery, 2018; Mummery & Balbus, 2019a), particularly at late times (Mummery & Balbus, 2019b). In the absence of a simple, physically realistic description of such discs, our reliance on the standard $\alpha$-disc picture is mainly a practical choice. From the perspective of modelling the outflow, the most important aspect of this approximation may be the angular distribution it imposes on the photons produced by the disc. Due to fore-shortening and limb-darkening, geometrically thin and optically thick discs generate a highly anisotropic radiation field. We adopt the Eddington approximation for limb-darkening, so that the specific intensity produced by the disc scales as $I\propto\cos{i}(1+\frac{3}{2}\cos{i})$ in our models. This degree of anisotropy is likely to be an overestimate for the vertically extended inner discs in rapidly accreting TDEs. Given that the outflow itself reprocesses, and hence isotropises the disc radiation field, we again do not expect this to change our results qualitatively. However, a more physically realistic description of the accretion discs in TDEs may be needed to enable meaningful quantitative modelling of observational data.

Finally, in Python, only H and He are self-consistently described with full, multi-level model atoms. By contrast, bound-bound transitions in metals are currently treated via a 2-level atom approximation that is reasonable for resonance lines, but not for transitions involving excited and/or meta-stable states. As a result, our simulations are not able to produce semi-forbidden transitions at the moment, such as N iii] or C iii], which are present in the data and could, in future, be useful outflow diagnostics. We plan to add an approximate way to handle some of these transitions in the near future.