© T. Wevers et al. 2022

1. Introduction

Double-peaked emission lines, usually seen in H Balmer and He II optical transitions, are observed in a small fraction (∼3%) of active galactic nuclei (AGNs; Chen et al. 1989; Chen & Halpern 1989; Eracleous & Halpern 1994; Strateva et al. 2003). Although a number of possible explanations for their origin exist in the literature, including binary supermassive black holes (SMBHs; Begelman et al. 1980; Gaskell et al. 1983), bipolar outflows (Norman & Miley 1984; Zheng et al. 1990), or highly anisotropic continuum sources (Goad & Wanders 1996), the leading explanation is that they originate in the outer parts of an inclined accretion disk (several thousand gravitational radii, where the gravitational radius R_g = $\frac{G{M}_{\mathrm{BH}}}{c^2}$; see e.g. Eracleous & Halpern 2003 for a detailed discussion).

The original literature models (e.g. Chen et al. 1989; Chen & Halpern 1989) envisaged a circular accretion disk and were successful in reproducing ∼40% of the known samples (Eracleous & Halpern 1994; Strateva et al. 2003). The observed sample morphology of double-peaked AGNs is diverse, including both stronger blue than red peaks and vice versa (see e.g. Eracleous & Halpern 2003 for an overview); the latter in particular cannot be explained with a circular accretion disk model alone. Some show line profile variability, observed on timescales from a few days (the ‘reverberation’ timescale; Schimoia et al. 2015) up to months and years (the ‘dynamical’ timescale; e.g. Gezari et al. 2007; Schimoia et al. 2017). Motivated by this diverse behaviour, more sophisticated and, more importantly, non-axisymmetric accretion disk models were developed (Eracleous et al. 1995; Strateva et al. 2003). Such elliptical models were able to reproduce the majority (but again, not all) of the sources where the circular models failed.

Two main hypotheses were put forward by Eracleous et al. (1995) for the formation of such eccentric accretion disks around SMBHs: binary SMBHs, where the eccentricity of the disk is pumped by the tidal torques of the secondary; and tidal disruption events (TDEs), in which the formation of a highly eccentric structure is a natural expectation in the absence of a mechanism to efficiently and rapidly remove angular momentum from the stellar debris. The vastly different timescales involved provide a mechanism to discriminate between these hypotheses through spectroscopic monitoring. In particular, binary SMBH disks are expected to evolve in thousands of years, whereas in TDEs evolution can be expected on timescales of weeks to months.

There are some previous claims in the literature for the presence of disks with significant eccentricity following TDEs. Cao et al. (2018) model the prototypical TDE ASASSN–14li, finding a large disk (r_out ∼ 1700 R_g) and an eccentricity e = 0.97. However, the profiles in this event are single-peaked and, moreover, can also be modelled as an optically thick, spherically symmetric outflow, where the line evolution is explained through electron scattering depth variations (Roth & Kasen 2018). Given the absence of significant asymmetries and double-peaked profiles, the evidence for an accretion disk origin of the emission lines is unclear. Using the same model, Liu et al. (2017) modelled the TDE PTF–09djl (see also Arcavi et al. 2014); for this source the data are sparse and noisy, but the Hα line does appear strongly asymmetric. It can be fit with a compact, highly elliptical (e = 0.96) accretion disk. Unfortunately, no similar line profile is found in other lines (e.g. H Balmer lines, He I, or He II). The absence of He II is explained through the inferred high inclination angle (an idea that is not compatible with the conclusions of this work), but the difference between the Hα and Hβ line profiles is more difficult to explain.

With increasing TDE detection rates and spectroscopic follow-up datasets, clearer evidence has emerged to associate optical emission line profiles directly with an accretion disk. Wevers et al. (2019a) and Cannizzaro et al. (2021) reported narrow Fe II emission lines likely associated with a disk chromosphere, while Holoien et al. (2019a, PS18kh), as well as Short et al. (2020) and Hung et al. (2020, AT 2018hyz), reported on flat-topped or double-peaked Hα (and other H Balmer) emission line profiles that are very likely disk-related (although see Hung et al. 2019 for an outflow scenario to explain the line profiles in PS18kh). In both PS18kh and AT 2018hyz, the inferred eccentricities are low (∼0.1–0.2) and uniform, and the disk inclinations are low to moderate (20–60 degrees). This may appear somewhat surprising given the current lack of understanding of the detailed dynamics of the post-disruption debris; in particular from a theoretical point of view, it is unclear how the stellar debris can shed its (expected) large amount of energy in such a short timescale to form a quasi-circular disk (e.g. Krolik et al. 2020). In the absence of such a mechanism, the naive expectation is for the debris to form a highly elliptical structure. Using hydrodynamical simulations, Shiokawa et al. (2015) found that the returning debris is unlikely to settle into a compact, circular disk, but instead forms an extended eccentric accretion flow. Piran et al. (2015) elaborated upon these results by showing that this is consistent with the relatively small amount of energy released in stream self-intersection shocks, and furthermore that such an elliptical disk scenario can provide a natural explanation of the observed properties (e.g. luminosity, temperature, and line widths; see also Krolik et al. 2016; Svirski et al. 2017; Ryu et al. 2020a; Zanazzi & Ogilvie 2020).

In this work we present the analysis of photometric and spectroscopic data of a new TDE, AT2020zso. We describe the observations and their data reduction in Sect. 2. Our analysis methods and results are presented in Sect. 3, and we discuss these results and their implications for accretion disk formation in TDEs, as well as AGNs, in Sect. 4. We summarise our conclusions in Sect. 5. Figures of the full posterior distributions for all model fitting results are provided in the appendix, along with a table containing all the photometry used. We assume a flat Λ cold dark matter cosmology with H₀ = 67.11 km s⁻¹Mpc⁻¹, Ω_m = 0.32, and Ω_Λ = 0.68 (Planck Collaboration XVI 2014) throughout the article.

2. Observations and data reduction

AT 2020zso was first reported as a transient by the Zwicky Transient Facility (ZTF20acqoiyt, Forster et al. 2020), and also detected by ATLAS (ATLASbfok; Smith et al. 2020) and Gaia (Gaia20fqa, Hodgkin et al. 2021). The host galaxy (SDSS J222217.13-071558.9) is an elliptical galaxy located at a redshift of z = 0.0563. A classification spectrum (Gromadzki et al. 2020) was obtained as part of the extended Public ESO Survey for Transient Objects (ePESSTO+; Smartt et al. 2015), and further spectroscopic follow-up was triggered within ePESSTO+. A detailed observing log is presented in Table 1. All phases are reported with respect to the phase of peak light (measured from the bolometric light curve) at MJD 59 184. The optical spectroscopy will be made publicly available through WISErep.

2.1. Spectroscopy

2.1.1. New Technology Telescope/EFOSC2

Low-resolution optical spectra were taken with the ESO Faint Object Spectrograph and Camera (EFOSC2) spectrograph mounted on the New Technology Telescope (NTT) at La Silla Observatory, Chile as part of the ePESSTO+ collaboration. We used the Gr#11, Gr#13, and Gr#16 grisms and a 1 or 1.5 arcsec slit width. The data reduction is performed using a dedicated pipeline (Smartt et al. 2015), which includes standard tasks such as bias-subtraction, flat-fielding, and a wavelength calibration based on arc frames and a comparison to sky emission lines. Cosmic rays are removed using the lacos routine (van Dokkum et al. 2012). To minimise host galaxy contamination, the source extraction was performed using an extraction aperture of 1 arcsec. For some epochs (in particular, those with slid widths > 1 arcsec) the seeing was > 1 arcsec, which may lead to different galaxy light contamination in these spectra. The Gr#16 observation is dominated above 7000 Å by second order contamination and is not used in our analysis. The flux calibration and extinction correction are performed using standard star observations.

2.1.2. Very Large Telescope/X-shooter

Shortly after the first EFOSC2 spectrum we triggered target-of-opportunity (ToO) observations with X-shooter, mounted on the Very Large Telescope (VLT) Unit 3 (Melipal) at Paranal Observatory, Chile. A total of four spectra were obtained using slit widths of 1.0, 0.9 and 0.9 arcsec for the UVB, visible (VIS), and near-infrared (NIR) arms, yielding a spectral resolution of R = 5400, R = 8900 and R = 5600, respectively. The data were taken in on-slit nodding mode. To increase the signal-to-noise ratio (S/N) of the UVB and VIS arms, we reduce these data using the X-shooter pipeline with recipes designed for stare mode observations. The NIR arm is reduced with both the stare and the nodding mode X-shooter pipeline recipes. The latter method provides a slightly better sky subtraction. Regardless of the method used, no transient emission features are found in the NIR spectra; they are shown in Fig. A.1 for completeness. For uniformity with the EFOSC2 spectra and to minimise host galaxy contamination, we adopt an extraction box with side 1 arcsec. Only for the last epoch (in which no TDE signal appears to be present) we use a 2 arcsec extraction box to boost the galaxy signal and determine the host galaxy properties. We use the molecfit software (Smette et al. 2015) to calculate atmospheric profiles and subtract telluric absorption bands in the VIS arm, which contaminate a small region redward of the Hα rest wavelength. The deep absorption band around 7300 Å is not well corrected, but does not contain any important emission lines. Figure A.1 shows the flux-calibrated X-shooter spectra.

2.1.3. Las Cumbres Observatory/FLOYDS

Two spectra were obtained with the low-resolution FLOYDS spectrograph mounted on the Las Cumbres Observatory 2m Faulkes Telescope North in Haleakala, Hawaii. The spectra were reduced using the floydsspec custom pipeline, which performs flux and wavelength calibration, cosmic-ray removal, and spectrum extraction¹.

2.1.4. Nordic Optical Telescope/ALFOSC

One epoch of spectroscopy was obtained using a ToO program on the Nordic Optical Telescope (NOT) in La Palma, Spain. This spectrum was taken with the Alhambra Faint Object Spectrograph and Camera (ALFOSC) spectrograph in combination with Grism 4 and a 1 arcsec slit. This observation was reduced using custom scripts based on the pypeit Python package (Prochaska et al. 2020b,a).

Following the standard data reduction recipes, we normalise all spectra to the continuum by fitting low order spline functions to the spectra, excluding known host galaxy and transient emission and absorption lines such as the He IIλ4686, He Iλ5876, and Hα regions.

2.2. Photometry

We retrieve the public Zwicky Transient Facility (ZTF) photometry via the ZTF forced-photometry service (Masci et al. 2019). The multi-band light curves are shown in Fig. 1.

Fig. 1. Host-subtracted light curves of AT 2020zso as observed by Swift, ZTF, ATLAS, and LCO. Vertical lines indicate epochs of spectroscopic observations. As a result of colour evolution, the peak in the optical bands occurs significantly later than in the UV bands. This is consistent with the observed cooling of the blackbody temperature over time.

Following the spectroscopic classification, Swift follow-up observations were triggered. The Swift Ultra-violet Optical Telescope (UVOT) photometry is measured using the uvotsource task in HEAsoft package v6.29 using a 5 arcsec aperture. Because no X-ray source was detected in the first observations, we derive an upper limit to the X-ray flux using the online X-ray Telescope (XRT) tool². Combining all observations, we find an upper limit of 1.56 × 10⁻³ cts s⁻¹, which translates into a flux of 4.5 × 10⁻¹⁴ erg cm⁻² s⁻¹, assuming a thermal (blackbody) spectral model with a temperature of kT = 75 eV³, typical for the soft X-ray emission in TDEs. This translates into a luminosity upper limit of 3.8 × 10⁴¹ erg s⁻¹ in the 0.3–10 keV band, uncorrected for foreground Milky Way extinction. Assuming instead an AGN like power-law spectral model (with power-law index Γ = 1.7), this translates into an upper limit of 5.3 × 10⁴¹ erg s⁻¹ in the 0.3–10 keV band, and 4.3 × 10⁴¹ erg s⁻¹ in the 3–20 keV band.

2.2.1. Las Cumbres Observatory

Las Cumbres Observatory BVgri-band data were obtained using the Sinistro cameras on Las Cumbres 1m telescopes. Point-spread-function fitting was performed on host-subtracted images using the lcogtsnpipe pipeline (Valenti et al. 2016), which uses HOTPANTS (Becker 2015) for the subtraction, with template images obtained also at Las Cumbres after the event faded. BV-band photometry was calibrated to the Vega system using the AAVSO Photometric All-Sky Survey, then converted to the AB system using the corrections from Blanton & Roweis (2007), while gri-band photometry was calibrated to the AB system using the Sloan Digital Sky Survey (Smith et al. 2002).

2.2.2. NIR photometry

Two epochs of NIR photometry were taken on 2021 May 22 (MJD 59 356) and 2021 July 15 (MJD 59 416), at phases +172 and +232 days after peak light. The first epoch, comprising observations in the H (three series of six dithered 20 s exposures, for a total of 1440 s on source) and K_s (also 1440 s exposure time) bands, was taken using the Son of Isaac instrument mounted on the NTT in La Silla, Chile. The reduction and combination of dithered images were carried out with the PESSTO pipeline. The second epoch of observations, including J, H and K_s band observations, was taken with the NOTCam instrument mounted on the NOT in La Palma via the NUTS2 programme. The NOTCam data were reduced using a version of the NOTCam Quicklook v2.5 reduction package⁴ with a few functional modifications (e.g. to increase the FOV of the reduced image).

In order to check for infrared variability, we performed aperture photometry on the NIR images. We measured the brightness of the central regions of the host galaxy with a 2 arcsec aperture and calibrated the resulting magnitude against the magnitudes of field stars taken from the Two Micron All-Sky Survey (2MASS) catalogue. The measurements were consistent within the measurement uncertainties, which were typically ∼0.05 mag. We do not find any significant brightening in either epoch. Similarly, no brightening is observed in the Near-Earth Object Wide-field Infrared Survey Explorer (NEOwise) mid-infrared light curves, that is, there is no evidence for an ongoing infrared echo.

2.3. Radio observation

A weak radio detection by the Very Large Array was reported on 2020 December 21 (phase +20 days, Alexander et al. 2021). A reported flux density of 22 ± 7 micro-jansky at 15 GHz corresponds to a monochromatic luminosity of 2.73 ± 0.87 × 10³⁷ erg s⁻¹.

2.4. Gaia astrometry

AT 2020zso was detected by the Gaia Photometric Science Alerts (GSAs; Hodgkin et al. 2021) as Gaia20fqa at coordinates (RA, Dec) = (22:22:17.130, –07:15:59.08). This allows an accurate evaluation of the positional offset with respect to the host galaxy nucleus, which has positional coordinates listed in the Early Gaia Data Release 3 (Gaia Collaboration 2021). An offset of 42 milli-arcseconds (mas) is measured, which corresponds to 46 parsec (pc) at the host redshift. The per-transit accuracy of GSAs is 55 mas (Wevers et al. 2019a; Hodgkin et al. 2021), leading to an offset of 46 ± 60 pc, consistent with a location in the nucleus of the galaxy.

3. Analysis and results

3.1. Host galaxy

3.1.1. Spectral energy distribution

We compiled the host galaxy (SDSS J222217.13-071558.9) spectral energy distribution (SED) using archival observations in the UV through infrared bands (see Table 2 and Fig. 2). In the NIR we used 2MASS (Cutri et al. 2003) J, H, and K_s magnitudes, and we used the Pan-STARRS DR1 magnitudes in g, r, i, y, and z optical bands. Finally, for the UV we performed aperture photometry in the Galaxy Evolution Explorer (GALEX; Bianchi et al. 2011) NUV and FUV images with the gPhoton package (Million et al. 2016). We performed forced aperture photometry in all available bands with two distinct apertures. For one set, we used 5 arcsec apertures; this was used to subtract the host contribution from the observed light curves (in particular, the Swift photometry is performed with a 5″ aperture). For the other set, we used an elliptical aperture with major and minor axes of 13 and 9 arcsec, optimised to include the entire host galaxy flux; this was used to model the host galaxy and derive its properties.

Fig. 2. Host galaxy SED and best-fit template used to synthesise host galaxy magnitudes in the Swift bands. The data are tabulated in Table 2.

We modelled the SED using the flexible stellar population synthesis (Conroy et al. 2009) module. We used the Prospector (Johnson et al. 2021) software to run a Markov chain Monte Carlo (MCMC) sampler (Foreman-Mackey et al. 2013). We assume an exponentially decaying star formation history (SFH), and a flat prior on the five free model parameters: stellar mass (M_⋆), stellar metallicity (Z), V-band extinction (A_V, assuming the extinction law from Calzetti et al. 2000), the stellar population age (t), and the e-folding time of the exponential decay of the SFH (τ_sfh).

Using the median and 1-σ confidence intervals of the posteriors of the fit to the 9″ by 13″ aperture photometry, we derive a host stellar mass of $log(M_{\star }/M_{\odot })=10.{11}_{-0.02}^{+0.02}$, a metallicity of $log(Z/Z_{\odot })=-0.{34}_{-0.07}^{+0.04}$, $A_V=0.{23}_{-0.04}^{+0.03}$ mag, $t=1.{90}_{-0.94}^{+1.56}$ Gyr, and $log({\tau }_{\mathrm{sfh}})=2.{42}_{-0.32}^{+0.35}$ Gyr. The extinction is roughly consistent with the Galactic foreground extinction of E(B − V) = 0.06 (Schlafly & Finkbeiner 2011). The estimated mass combined with the extinction-corrected rest-frame colour u − r = 1.89 ± 0.02 mag places the host galaxy near the ‘green valley’ region (Schawinski et al. 2014) of the mass colour diagram, in which the TDE host galaxies seems to be over represented compared to the general galaxy population (Law-Smith et al. 2017; van Velzen et al. 2021; Hammerstein et al. 2021).

To estimate the host galaxy fluxes in the UVOT bands, we similarly model the host galaxy SED but using the 5″ aperture data. The host contribution is then subtracted from the measured photometry, which is also corrected for foreground Galactic extinction. The uncertainty on the host galaxy model is propagated into our measurement of the host-subtracted TDE flux (see Table 2).

3.1.2. Black hole mass

Using the late time X-shooter spectrum, in which no broad emission lines are apparent, we measure the velocity dispersion of the host galaxy following the method of Wevers et al. (2017), using the penalised pixel fitting routine (Cappellari 2017). We find a velocity dispersion of σ = 60 ± 1 km s⁻¹. Using the M − σ relation from McConnell & Ma (2013) this translates to a black hole mass of log(M_BH) = 5.3 ± 0.4, or alternatively log(M_BH) = 6.2 ± 0.3 M_⊙ using the Kormendy & Ho (2013) relation, indicating a low mass black hole similar to many other UV/optical discovered TDEs (Wevers et al. 2019b).

3.1.3. Emission lines

From the narrow host galaxy emission lines, which are resolved in the X-shooter spectra, we measure a redshift of z = 0.0563, which corresponds to a luminosity distance of 263 Mpc.

We identify a plethora of narrow emission lines originating in the host galaxy, including in order of increasing wavelength: the [O II] λλ3726, 3729 doublet, [Ne III] λ3869, He IIλ4686, Hβ, the [O III] λλ4959, 5007 doublet, He Iλ5876, [O I] λ6300 line, the [N II] λλ6548, 6584 doublet, Hα, the [S II] λλ6717, 6731 doublet, a (very) weak [Ar III] λ7136 line, the [S III] λλ9069, 9532 doublet lines, and Pa αλ1.875 μm. High-ionisation-potential lines such as He II, [Ar III] and [S III] indicate that a hard photo-ionising continuum source is present, while there is no sign of a broad component to any of these lines at late times. We measure a full width at half maximum (FWHM) from the [O III] line of 159 ± 3 km s⁻¹, while for the narrow Hα, Hβ, N II, and S II lines we measure an average of FWHM = 127 ± 5 km s⁻¹. Closer inspection shows that some of these narrow lines are asymmetric/double-peaked, with a velocity separation of ∼70−80 km s⁻¹. Figure 3 shows some of the prominent narrow emission line profiles. We measure the asymmetry of the [O III] λ5007 line, the strongest narrow emission line, by using the non-parametric measurement of Liu et al. (2013), and find a very small asymmetry A = 0.047 (other lines yield similar values). This suggests that the narrow line region is rotation dominated, but probably not kinematically disturbed (Blecha et al. 2013; Nevin et al. 2016), which makes it very unlikely that the system hosts a dual AGN or a wide-separation SMBH binary. We measure an [O III] line luminosity from the late time spectrum of L_[OIII] = 1.17 ± 0.05 × 10⁴⁰ erg s⁻¹. Assuming a correlation between L_[OIII] and the 3–20 keV X-ray luminosity observed in AGNs (Heckman et al. 2005), we expect an AGN X-ray luminosity of L_X = 1.6 × 10⁴² erg s⁻¹. The upper limit for L_X derived from Swift observations is 4.3 × 10⁴¹ erg s⁻¹ in the 3–20 keV band, which is marginally inconsistent given the large scatter (a factor of ≈3) in the correlation. Reconciling these two values (again assuming a power-law spectrum with index Γ = 1.7) would require an absorbing column of at least ∼1–1.5 × 10²² cm⁻².

Fig. 3. Insets of the narrow emission lines observed in AT2020zso.

We also use the narrow host galaxy line ratios to put the source on a Baldwin-Philips-Terlevic (BPT) diagram (Baldwin et al. 1981). We find that the source falls well into the AGN/Seyfert part of the diagram, in line with the presence of high-ionisation narrow emission lines (Fig. 4). We conclude that this galaxy hosts a Seyfert AGN.

Fig. 4. Narrow emission line ratios as measured from the late-time X-shooter spectrum. The results place the host galaxy firmly into the AGN (Seyfert) part of the diagram. Dividing lines are taken from Kewley et al. (2001) and (2006).

3.2. Light curve evolution

The host-subtracted light curves are shown in Fig. 1. The best sampling is achieved in the ZTF bands, particularly at very early times. There appear to be three distinct phases in both the g- and r-band light curves: a very steep initial rise, followed by a break to a slower increase in brightness and finally a turnover to a decline in brightness. To characterise the light curve behaviour at early times, we fit a power-law model to the two parts of the rising ZTF light curve (before and after the break) independently, of the form

$$\begin{array}{c}\hfill L=a+b\times {(t-t_0)}^{\alpha }.\end{array}$$(1)

We find that the early rising part of the light curve is consistent with L ∝ t² evolution (α = 1.9 ± 0.4). After the break (which happens around phase = –12 days), the slope flattens to α = 1.55 ± 0.25. We remark that the emission lines contribute < 10% of the total light, and therefore do not significantly influence the light curve evolution nor the inferred parameters.

3.2.1. Parameter inference with MOSFit and TDEMASS

We fit the light curve using a TDE model in the MOSFit package (Guillochon et al. 2018; Mockler et al. 2019), employing the same free parameters and priors as in Mockler et al. (2019). We used the DYNESTY dynamic nested sampling algorithm with default stopping criteria to explore the parameter space and sample the model posteriors (see Speagle 2020, for details). The MOSFit TDE model only includes the fallback luminosity as an energy source, while at later times ≳100 days (e.g. Mummery & Balbus 2020a) there may be a significant contribution from an accretion disk, leading to a flattening of the light curve. We therefore exclude data points more than 100 days after peak luminosity from the fit, while noting that a fit including these data does not significantly alter the black hole mass and stellar mass estimates. Fits were run on the University of Birmingham BlueBEAR cluster.

The results are a poor fit (Fig. 5); there are short-term variations that are not encapsulated by the model itself, so these are not expected to be reproduced, and the temperature variation is more rapid than the model can accommodate. These results should therefore be interpreted with some caution.

Fig. 5. MOSFit model light curves overlaid on the data. The fits are not great, likely due to the relatively rapid temperature evolution and peculiar light curve behaviour (including a break in the rising part), which the model cannot accommodate.

A black hole mass of log₁₀(M_BH) = 6.0 ± 0.3 is inferred from the light curve, which is consistent with the estimate from the stellar velocity dispersion. Furthermore, we obtain estimates of the disrupted stellar mass (at the lower allowed limit of ∼0.1 M_⊙) and the impact parameter $\beta \equiv \frac{R_{\mathrm{p}}}{R_{\mathrm{t}}}=0.63\pm 0.05$ (where R_p is the orbital pericentre radius, and R_t the tidal radius), although there are large systematic uncertainties of 0.66 dex or 0.36 in a linear scale for a value of 0.1 for the stellar mass, and 0.35 for the impact parameter (see Mockler et al. 2019 for a detailed discussion of the systematic uncertainties produced by MOSFit). The results of these fits are reported in Table 3, and full posterior distributions for the fits can be found in the Appendix. These values are very similar to those found by Gomez et al. (2020) for the other double-peaked TDE AT 2018hyz, and indicate that AT 2020zso may likewise be the result of a partial, rather than a full, disruption (as inferred from the fact that β < 0.9; e.g. Guillochon & Ramirez-Ruiz 2013). Simulations by Ryu et al. (2020b) suggest that the surviving stellar remnant may have lost ∼40 per cent of its original mass (but keeping in mind the large systematic uncertainties, this could range from a few up to > 60 per cent). Finally, MOSFit suggests a disruption date of 8 ± 2 days before the first data point, at MJD 59 157 ±2 days.

Alternatively, Ryu et al. (2020a) presented a framework to infer the black hole and disrupted stellar masses on the basis of eccentric accretion disk dynamics. Using the peak bolometric luminosity of 7 ± 1 × 10⁴³ erg s⁻¹ and a peak colour temperature of 25 000 ± 5000 K, a black hole mass of $M_{\mathrm{BH}}=1.7_{-0.9}^{+2.0}\times {10}^6$ M_⊙ and stellar mass of $M_{\star }=0.{92}_{-0.13}^{+0.21}$ M_⊙ are inferred; the former is consistent with alternative estimates from galaxy scaling relations, while the latter differs significantly from the MOSFit estimate. Given that the TDEMASS framework was explicitly developed on the basis of eccentric accretion disk dynamics, which appear to be particularly suited for application to AT2020zso, we give preference to these inferences.

3.2.2. Blackbody modelling

The properties of UV/optical TDE flare emission can be empirically be well described using evolving blackbody models. This is somewhat surprising, given that there is growing evidence for a viewing angle dependence of the observational consequence of stellar destruction (Dai et al. 2018; Leloudas et al. 2019), implying that asymmetry is present in the ensuing structure. As a result, it is unclear to what degree the results of blackbody fitting the TDE SED can be physically interpreted.

With this caveat in mind, for each epoch of Swift observations, we model the SED using a blackbody curve, although the physical interpretation of this model is unlikely to be straightforward, as the powering source of this emission remains unclear. We include all the host-subtracted Swift photometry, and linearly interpolate the ZTF g- and r-band measurements to these epochs to provide coverage from 2000–7000 Å. We do not extrapolate the ZTF measurements beyond their latest observing epochs. We use a maximum likelihood approach to fit a blackbody model to each epoch, assuming a flat prior for all parameters; the resulting fits are illustrated in Fig. 6. Uncertainties are assessed by sampling from the posterior distributions of the parameters directly. Assuming isotropic emission also yields the characteristic blackbody radius. The temperature, radius and bolometric blackbody luminosity are shown in Fig. 7; the values based only on ZTF (without temperature fit, but with a bolometric correction) data are shown as green triangles.

Fig. 6. Blackbody fitting of the TDE SED over time. The uncertainties become very large for the last three epochs as two data points do not provide much constraining power for the blackbody model.

Fig. 7. Blackbody parameter evolution of AT2020zso over time. Strong cooling is observed during the first part of the light curve.

The temperature decreases by ∼10 000 K during the first part of the light curve. Several sources in the van Velzen et al. (2021) sample show similar behaviour. Such cooling is typically seen in TDEs after peak, likely as a result of an expanding photosphere, although the effect is particularly strong for AT 2020zso (a similar effect was seen in AT2019qiz; Nicholl et al. 2020). The MOSFit results show a similar temperature evolution, although the temperature changes somewhat slower (this is intrinsic to the TDE model). There is some indication of an increasing trend at later times (again similar to typical TDE behaviour), but the uncertainties are large.

The photosphere radius evolution follows a linear expansion profile initially. We measure an expansion velocity of v_exp = 2900 ± 300 km s⁻¹ before peak light (Fig. 8, top panel). Afterwards, the radius reaches a plateau before moving back inwards to scales ≈5 × 10¹³ cm. This behaviour is very similar to that observed in AT2019qiz (Nicholl et al. 2020) and AT2019ahk (Holoien et al. 2019b). Based on the expansion velocity before maximum, we estimate that the first observations were taken approximately 15 days after expansion began, suggesting a disruption date around MJD 59 149, and a rise time of approximately 35 days from disruption to peak (compared with an explosion date of MJD 59 157 from MOSFit). This value is insensitive to the assumption about the temperature evolution before the first Swift observations.

Fig. 8. Blackbody radius evolution over time. Top panel: evolution of the blackbody radius, overlaid with the best-fit linear expansion model (orange line). At late times the radius reaches approximately 5 × 1013 cm (but with large error bars), which corresponds to several hundred Rg for MBH = 106 M⊙. Bottom panel: same evolution, but now converted to units of gravitational radius for two different black hole masses. The orange and red stars mark the inferred disk outer and inner radii, respectively (see Sect. 3.6). At late times, the blackbody radii overlap with the estimated accretion disk size (highlighted with the shaded area, assuming that the disk size has not changed significantly from the last epoch that can be modelled).

We note that an expanding photosphere does not necessarily require a physical outflow (i.e. outward fluid motion) to be present. Alternatives to explain the photosphere expansion include the accumulation of matter around the peak of the mass fallback rate, which extends the photosphere to larger radii; it could be the result of time-dependent photon diffusion due to changing density and/or optical depth in the debris; or due to the orbital motion of heated matter, which at the inferred radius of ∼5 × 10¹⁴ cm is comparable (∼5000 km s⁻¹) to the measured growth rate.

3.3. Transient emission features

In addition to the narrow host galaxy lines, broad evolving emission lines are present in the spectra, which are shown in Figs. 9 and 10. These lines are typically (quasi)-Gaussian in TDEs, with velocity shifts up to ∼15 000 km s⁻¹. However, in AT 2020zso there are several emission features whose identification would be contrived, or completely unclear, when taking this approach. For example, a broad feature centred on 4500 Å appears around –4 d. While this could in principle be broad Fe II emission often seen in AGNs, the profile appears smooth with a broad blue wing that would be atypical. Similarly, broad features centred on 4250 Å, 5050 Å, 6080 Å, and 6820 Å are present in several of the spectra. No similar features have been readily identified in other TDEs to date, and no immediately obvious line identifications are available for these wavelengths.

Fig. 9. Spectral evolution of AT 2020zso from early to late times. All spectra are continuum normalised. The different emission features in the 4600 Å blend are indicated with vertical lines: solid black lines indicate H Balmer lines, dashed black lines N III, dotted black lines Fe II, red dashed lines He II, and red dotted lines He I. Black spectra were taken with NTT/EFOSC2, blue with VLTX-shooter, orange with LCO/FLOYDS, and magenta with NOT/ALFOSC. The X-shooter spectra are resampled at lower spectral resolution for clarity. The gap around 5200 Å is due to the low S/N of the data at the edges of the UVB and VIS arms, and the gap around 7300 Å is due to the poor correction of molecular bands.

Fig. 10. Same as Fig. 9, but now in velocity space. Left panel: the He II region, and the right panel shows the region around Hα.

Noting that these broad emission features appear to be roughly symmetrical around rest wavelengths of ∼4700 Å and 6560 Å, we instead explored the idea that these are multi-peaked structures constituting a single feature, that is, emission lines of He II/Hβ and Hα. This is motivated by previous studies that have identified double-peaked emission lines, attributed to accretion disk structures, in TDEs (Arcavi et al. 2014; Short et al. 2020; Hung et al. 2020).

The velocity structure of these lines is shown in Fig. 11; it is encouraging that the profiles appear very similar in this representation. We identify the emission features with two main contributors, consistent with rest wavelengths of He IIλ4686 and Hα. Because Hα appears centred near rest velocity, we disfavour an identification as Hβ for the emission feature near 4700 Å; such an identification would require a large systematic blueshift (–11 000 km s⁻¹), whereas for He II the line would also be centred near rest velocity. Nevertheless, the broad feature around 5050 Å may identify as the red wing of a similar double-peaked velocity profile consistent with Hβ, albeit much weaker than He II. No other broad He II emission lines (e.g. at λλ3203,10 123 Å) are evident in the spectra. We also note that there are other emission features seen in TDEs around this region, most notably the Bowen N III 4640 line. We explore possible contamination in more detail in Sect. 3.5.

Fig. 11. VLT/X-shooter spectra at various epochs, plotted in velocity space for He II (black) and Hα (orange). To guide the eye, we show: the rest velocity as a solid line; the photospheric velocity with a vertical dashed line, which happens to coincide with the expected position for N IIIλ4640; and the position for Hβ with respect to He II as a black dotted line. For display purposes, narrow emission lines have been removed and the spectra are resampled to lower spectral resolution.

Given the strongly non-Gaussian line profiles, we measure the emission line equivalent widths (EWs) of the He II complex and Hα through direct integration. We mask out telluric absorption features in the spectra when present. The EW evolution and their ratio is shown in Fig. 12. The He II/Hα ratio decreases rapidly from ≈11 (phase –14 days), to 7 (phase –13 days), to stabilise around 1.5 (phases later than –10 days). For the earliest epochs, we also measure the EW ratio by using a Gaussian profile, which yields somewhat lower values (≈8 and 4 at –14 and –13 d, respectively) but shows a consistent, rapidly decreasing trend.

Fig. 12. Emission line measurements and their evolution over time. Top panel: equivalent width of the He II complex (black circles) and Hα (orange triangles), measured through direct integration from the spectra. Bottom: their ratio (He II/Hα) as a function of time.

The first four epochs are well described (reduced χ² < 1.3) by a broad, single Gaussian, and for these spectra we also attempt to measure the line velocities and FWHM of He II and Hα. The results indicate that both lines are likely at rest velocity. All measurements are consistent with 0 within 3σ, although we find large variations between different spectra that are unphysical given that they were taken only days apart, likely due to the relatively low S/N. The ‘best’ early spectrum available was taken with X-shooter (phase –13 days), from which we measure line velocities of 800 ± 400 km s⁻¹ and –170 ± 135 km s⁻¹ for Hα and He II, respectively. From the same spectrum, we measure Gaussian FWHM values of 36 000 ± 7000 (Hα) and 31 000 ± 1000 (He II) km s⁻¹ (for a reduced χ² of 1.2).

3.4. Emission line evolution

We plot the broad emission line profiles of He II and Hα in Fig. 11. For clarity, we focus on the highest spectral resolution (X-shooter) spectra. We start by noting that there appears to be a delay between the emergence of He II (which appears earlier) and Hα, which is very weak in the earliest epochs but strengthens ∼10 days after the first observation (as is also apparent from the EW evolution in Fig. 12). The line profile of He II is slightly asymmetric and peaks around –3000 km s⁻¹, which would be consistent with the photospheric outflow velocity measured through blackbody modelling.

Near peak light (phase –3 days) the spectra become distinctly non-Gaussian. The He II line profile is inconsistent with the presence of a very broad Gaussian (similar to the one observed at –13 days) – this component would extend well beyond 4300 Å, where no excess flux is observed. Both lines show a double-peaked structure, centred near rest wavelength; He II appears roughly symmetric whereas Hα shows a strong asymmetry, with a bright red peak and a broad blue shoulder (no clear blue peak is visible in the spectrum). The red peak of both lines occurs at similar velocities (∼8000 km s⁻¹), and the red wing is nearly identical in velocity structure, extending out to ≈26 000 km s⁻¹. On the blue side, the He II profile extends roughly 5000 km s⁻¹ further blueward (out to 18 000 km s⁻¹) compared to Hα, and is a factor of ∼2 brighter than Hα. This could indicate either that the He II emission originates from a region with a different velocity profile or potential contamination of other emission lines known to be present in this region (including N IIIλ4640 and Fe II lines), discussed in Sect. 3.5.

After peak light (phase +10 days), the profiles have significantly evolved. Hα shows a pronounced triple-peaked structure. In particular the red peak of the profile is remarkable, being significantly brighter and more narrowly peaked than its blue equivalent. Similar to the previous epoch, the red wings of Hα and He II have comparable velocity structures. However, the red peak has moved to higher velocities, particularly for Hα (∼11 500 km s⁻¹). He II remains broader and brighter in the blue wing, extending out to 21 000 km s⁻¹. The blue Hα peak is situated around –8000 km s⁻¹, whereas He II peaks closer to –10 000 km s⁻¹. The ‘central’ peak of Hα is more pronounced and near rest velocity.

We note that the photospheric radius has reached a plateau at this phase (+10 d), that is, the outflowing photosphere reaches its maximum radius before becoming optically thin and receding inwards. The origin of these central features could therefore be either from the accretion disk itself⁵ or (less likely) in this outflowing component near maximum radius, consistent with the low observed velocities.

It is worth noting that there is no pronounced central peak in He II, although some feature may be present. This is unlikely to be the blue wing of a double-peaked Hβ line profile, because (i) this would imply a velocity of –14 000 km s⁻¹, significantly larger than both the Hα and He II blue peak velocities, and (ii) for Hα the red peak is brighter than the blue peak, whereas such an identification would imply a stronger blue peak for Hβ. Alternatively, this feature could be consistent with Bowen N IIIλ4640.

This remarkable triple-peaked structure is also seen in He Iλ5876 emission lines (right panel of Fig. 11). An almost identical feature is also present near 4100 Å. This could be identified as Hδ, although this is inconsistent with the weakness of Hγ and Hβ. A more likely possibility is the N III Bowen line at 4100 Å (see Sect. 3.5). These profiles may also be present in the earlier X-shooter epochs, albeit very weak. The wings of these profiles extend from –10 000 km s⁻¹ to 12 000 km s⁻¹, so are significantly narrower than Hα and He II. This could imply that they originate in lower velocity regions of the disk, that is, farther out than Hα and He II. Their central peak also appears near rest velocity. We note that He II shows broader peaks compared to the profiles of the lower-ionisation lines. At the same time, the He II blue peak appears brighter than the red peak, which is opposite to the Hα profile. This is likely due to contamination of the Bowen N III line at 4640 Å, as quantified below.

3.5. De-blending the He II complex

As noted above, the line profiles of He II differ somewhat from those of Hα: in particular, the blue peak is either equal to or brighter than the red peak, whereas in Hα the red peak is always stronger than the blue peak.

The region around He II contains a number of emission lines that are observed in TDEs, including N IIIλ4640, Hβ and Fe II lines. To investigate this in more detail, we tried to fit the entire emission feature with a superposition of Gaussians for the aforementioned elements. We are not able to find consistencies in the derived parameters for the respective lines (e.g. line identification, velocity, or FWHM) if we include single Gaussian components in addition to He II.

Instead, we turn to the emission lines observed during the epoch at +10 days. In particular, the double-peaked feature centred on 4100 Å is unlikely to be Hδ, given the absence of both Hβ and Hγ. Instead, this line can be identified as a N III Bowen line. From a spectroscopic study of a sample of TDEs, it was found that the N III 4100 and 4640 Å lines have a roughly 1:1 flux ratio (Charalampopoulos et al. 2022) when they are present. Given that the 4100 Å line is readily identified in the X-shooter spectrum at +10 days, we attempt to subtract this line profile from the He II line, assuming a 1:1 flux ratio and identical velocity structure (so we assume that N III 4640 is identical to N III 4100). The red wing of the N III 4100 appears to be contaminated by another potential emission line. This feature is centred around 4250 Å, and it is not immediately clear to which element this line can be attributed. In order to not overfit the spectra, we normalised the entire UVB spectrum (in the range 3700–5250 Å) using a low order spline function, and this feature therefore remains in the normalised spectrum. The N III 4100 and He II features are shown in the top panel of Fig. 13. The contamination of N III manifests as the increasing flux trend at velocities above 10 000 km s⁻¹, and may lead to some systematic uncertainties in the fitting described later. To avoid biasing the results in the blue wing, we linearly interpolate over the Ca H+K region (shown by the dashed line in the top panel). The bottom panel of Fig. 13 shows the result of the subtraction (in blue) overlaid on the (unmodified) Hα line profile. The result of the red wing contamination leads to a somewhat narrower red wing profile in velocity space. Notwithstanding the simplifying assumptions for the various line profiles, the subtracted He II profile is remarkably similar to the Hα line profile. We conclude that AT2020zso belongs to the spectroscopic class of He+Bowen line TDEs, and that the He II line likely originates from the same physical region as the Hα and Bowen N III lines.

Fig. 13. Comparison of the emission line profiles. Top panel: He II (blue) and N III 4100 Å lines from the X-shooter spectrum at +10 days. Bottom panel: subtracted spectrum (He II – N III, blue) and Hα for comparison (in black). The line profiles are nearly identical, suggesting the presence of double-peaked Bowen N III lines, the first such line profiles seen in a TDE to date. The discrepancy in the red wing (compared to Hα) is likely related to the imperfect normalisation and subtraction of the N III profile.

We also performed the same subtraction procedure for for the X-shooter spectrum at phase –3 days and present the results in Fig. A.3 of the appendix. We use these subtracted He II line profiles for the accretion disk modelling. It is likely that the subtraction procedures introduce some systematic errors (e.g. the contamination in the red wing of N III), so we treat the fits to the Hα lines as our primary results. We also present the He II fitting results, but keeping in mind the caveats mentioned here.

3.6. Elliptical accretion disk fitting

Before we describe our model fitting results, we remind the reader that a general prediction of relativistic circular accretion disk models is that the blue peak is equal to, or brighter than, the red peak, due to Doppler boosting of emission from the region moving in our direction (e.g. Chen et al. 1989; Chen & Halpern 1989). This is no longer true for eccentric accretion disks: for non-negligible eccentricities, there exist orientation and inclination angle combinations such that the red peak can appear brighter than the blue peak. In other words, when interpreting a line profile as originating from an accretion disk, a double-peaked line profile with a dominant red peak is a strong indication of a non-axisymmetric (eccentric) configuration. There exist other mechanisms that can produce similar asymmetric line profiles, which will be discussed in Sect. 4.

Double-peaked emission lines have been observed in other TDEs, where they were interpreted as signatures of a spiral wave or an accretion disk. For the TDE with the most convincing double-peaked emission profiles so far (AT2018hyz; Short et al. 2020; Hung et al. 2020), an almost circular geometry was derived through model fitting. In our analysis we take a similar approach to Hung et al. (2020), and attempt to fit a general relativistic accretion disk model of Eracleous et al. (1995) to the Hα and He II lines. We do not model the N III and He I lines because of their limited S/N, and furthermore we do not present modelling results for the lines in low-resolution spectra because this leads to degeneracies and inconsistent results.

This model has seven free parameters, including the emissivity power-law index q, the broadening parameter σ, the major axis orientation of the elliptical rings ϕ₀ (for 0 degrees, the nodal line is along our line of sight, and the apocentre is in our direction), the inclination angle i (where 0 degrees is face-on), the eccentricity e, and the inner and outer pericentre distances (that is the line emitting region of the disk is bounded by elliptical annuli of radii r₁ and r₂), and is described by the following specific intensity profile I_ν:

$$\begin{array}{c}\hfill I_{\nu }=\frac{1}{4\pi }\frac{{ϵ}_0{\xi }^{-q}}{\sqrt{2\pi }\sigma }\exp [-\frac{({\nu }^2-{\nu }_0^2)}{2{\sigma }^2}].\end{array}$$(2)

The line profile flux F ≡ ∫dν∫∫dΩI_ν, that is, the integral of I_ν over frequency, specific intensity and solid angle, is calculated by numerical integration and rescaled to fit the observations. The fit is performed by varying two normalisation constants A and B, A + B × F (where A ≈ 1 is the adjacent continuum level for each spectrum and B is the amplitude of the profile) to the data to account for small (a few per cent) differences in the normalisation level.

We focus our analysis on the X-shooter spectra as they have superior spectral resolution and S/N. The first epoch (at –13 days) is omitted, because at these early times the line emitting regions likely originate in an outflowing photosphere rather than in an accretion disk like structure (as is implied by the EW evolution and will be discussed later). We nevertheless attempted a fit for this epoch, but the results are highly degenerate and do not allow robust parameter inference. We fit each emission line profile (He II and Hα) separately, for each epoch, for a total of four fits. We use a nested MCMC sampling approach implemented in dynesty, with uniform priors for all parameters as summarised in Table 4. Because the computational time of this approach is proportional to the size of the parameter space, we first create a grid of ∼231 000 models and perform a least-squares minimisation to assess the best-fit models. These results are used to inform the prior ranges, which are nevertheless taken very conservatively to encompass most of the plausible parameter space. Following the results presented in Hung et al. (2020), we report on the results of a composite accretion disk + outflowing component model fit to the data, where the outflow is represented by a Gaussian component.

The best-fit model parameters are presented in Table 5, and they are overlaid on the data in Fig. 14. The posterior distributions of the parameters for each fit can be found in Appendix A.

Fig. 14. Best-fit accretion disk plus Gaussian models, overlaid on the emission line profiles of He II (left panels) and Hα (right panels). Black represents the data, blue the accretion disk model, and orange the accretion disk plus outflow model. The best-fit Hα models are overplotted in the He II profiles as dashed lines for comparison. The red wing of He II is poorly subtracted, leading to some differences with the Hα profiles, but overall the results are consistent. The model parameters can be found in Table 5.

For each individual epoch, we infer global accretion disk model parameters that are very similar for the Hα and the He II profiles. This may appear somewhat surprising, given that the He II profile is significantly contaminated by N III. This indicates that the subtraction (and the 1:1 flux ratio assumption for the N III lines) performs well.

We find that the disk must be highly inclined with respect to our line of sight, that is, we are seeing it nearly edge-on (i = 85 ± 5 degrees for Hα, and i = 88 ± 2 for He II). The inferred eccentricity is high and consistent between Hα and He II, e = 0.97 ± 0.01. This is very similar to the characteristic eccentricity of ballistic orbits travelled by the tidal debris, ∼0.98 [$\frac{(M_{\mathrm{BH}}/M_{\star })}{{10}^6}$]^1/3, expected for returning debris whose orbit has not been significantly altered by hydrodynamics. We furthermore find largely consistent orientation angles (ϕ between 210 and 260 degrees) and line broadening parameters (σ∼1000–1500 km s⁻¹) for all epochs/lines. There is some scatter in the inner and outer disk radii, but this may not be surprising given the limited S/N of the data. The inner disk radius is several 100 gravitational radii, while for the outer radius we find values between 1400–2200 R_g, with one outlier at 4100 R_g. Finally, we note that while we have added a Gaussian component to the line profile fitting, in most cases the amplitude and contribution of this component is small. Only in the last epoch is there a clear triple-peaked structure. This additional component is found to be consistent with being at rest velocity, with a FWHM of 1000–2500 km s⁻¹.

In summary, we find that a highly inclined, highly elliptical accretion disk model can reproduce the Hα and He II line profiles of both epochs, with general disk parameters that are largely consistent within their uncertainties. Given the contamination He II, the fitting results are reasonably similar to those inferred from Hα.

4. Discussion

4.1. Alternative origins of the broad double-peaked emission lines

As previously indicated, there exist multiple mechanisms or structures that can explain the presence of double-peaked broad emission lines. We now discuss each of these in more detail, and why we prefer the accretion disk model as an explanation.

4.1.1. Supermassive black hole binary

An SMBH binary would spend most of its time in the hard binary phase, where the separation is typically 0.1−1 pc (10¹⁷⁻¹⁸ cm, Eracleous et al. 1995). In this scenario, the tidal field of the secondary would drive a disk around the primary to become eccentric if the mass ratio is > 4 (Eracleous et al. 1995). However, the timescale for this eccentricity to evolve is thousands of years, clearly inconsistent with the observed line profile variations on timescales of weeks.

We conclude that an SMBH binary alone cannot provide an explanation for the observed properties of AT 2020zso, in particular the rapid evolution of the line profiles.

4.1.2. Turn-on or changing-look AGN

One scenario that could help explain the observed properties is that rather than a newly formed accretion disk, we are seeing either a turn-on AGN (without the need for a tidal disruption) or a dormant, eccentric accretion disk being reinvigorated with fresh material, so an AGN turning on where a fossil accretion disk is being resupplied by the debris of a star. In either case, the binary SMBH hypothesis is then required to explain the observed eccentricity. This scenario could be consistent with the hypothesis that the primary black hole was a dormant AGN that shut off in the recent past, as inferred from the narrow line region diagnostics. In other words, the AGNs must have shut off within a light travel time (typically a few hundred years) to the narrow line region. Such timing would be coincidental; the high inferred eccentricities also imply that the origin of the emission lines is very unlikely to be a pre-existing broad line region that is reactivated by the flare, as the values inferred from double-peaked AGN sources are typically much more modest (e.g. Eracleous et al. 1995; Strateva et al. 2003).

4.1.3. Outflows and spiral arms

Bipolar outflows can also result in double-peaked line profiles, but the brighter red-than-blue peaked profile cannot be reproduced through Doppler boosting of emission. Nevertheless, profiles similar to AT2020zso (in that they have a brighter red than blue wing) have been observed in some supernovae (Smith et al. 2015; Bose et al. 2019), although they are generally seen in Hα, not He II. It is unclear what would power the Bowen fluorescence lines in this scenario. This scenario would require an ad hoc adjustment of the relative brightness of the blue and red peaks to produce the observed variability. Furthermore, it would require the UV/optical blackbody photosphere to expand and recede independently of the outflow (as we see the former moving inwards after peak, which is incompatible with an outflow scenario powering the lines at those times). The quasi-Gaussian line profile, combined with the light curve evolution at early times suggests that any outflow present in AT 2020zso would likely have a near spherical geometry. While a wide-angle bipolar outflow can therefore not be excluded based on current data, an aspherical structure could be detectable in polarimetric observations.

Spiral structures have also been invoked to help explain the variability in double-peaked AGN sources (e.g. Storchi-Bergmann et al. 2003). This variability is typically associated with the rotation of the gas (or precession of the spiral structure) on timescales larger than several dynamical timescales. This dynamical timescale is roughly

$$\begin{array}{c}\hfill t_{\mathrm{dyn}}=2M_6{\xi }^{3/2}\mathrm{days},\end{array}$$(3)

where M₆ is the black hole mass in units of 10⁶ M_⊙ and ξ is the disk outer radius in units of 1000 R_g. This yields a value of ∼5–15 days for AT 2020zso. While this appears compatible with the observed timescales for variability, it remains unclear how such a spiral structure would form in the very brief period of time between disruption and peak light. Furthermore, this scenario generally invokes an axisymmetric disk configuration, and hence circularisation would have to be extremely rapid – not accounting for the formation timescale of the spiral structure itself. We therefore deem it unlikely that spiral arm patterns can provide a plausible explanation of the observed behaviour.

4.2. Light curve and blackbody evolution

The light curve and radius evolution are remarkably similar to other TDEs with pre-peak observations, including AT 2019qiz (Nicholl et al. 2020) and AT2019ahk (Holoien et al. 2019b): consistent with L ∝ t² and constant outflow velocities of a few thousand km s⁻¹. A quasi-spherical outflow with constant velocity and temperature will lead to the observed L ∝ t² behaviour. The fact that the early evolution (before the first Swift observations) is consistent with L ∝ t² suggests that the temperature was roughly constant during this phase.

The temperature cools significantly over the first 40 days, behaviour that is similar to AT 2018hyz and ASASSN–14ae (Gomez et al. 2020) as well as ASASSN–15lh (Leloudas et al. 2016). This may be related to a comparatively low amount of debris due to a partial tidal disruption, leading to shorter diffusion times (Short et al. 2020; Gomez et al. 2020) and therefore faster temperature evolution. This cooling phase may also help explain the transition from initially broad Gaussian line profiles dominated by He II to the appearance of Hα slightly later, and finally to the emergence of the double-peaked disk profiles.

The peak bolometric UV/optical luminosity reaches 7 × 10⁴³ erg s⁻¹ – this corresponds to roughly the Eddington limit of a ∼5 × 10⁵ M_⊙ black hole, or an Eddington ratio of 0.4 for a ∼10⁶ M_⊙ black hole. The blackbody radius (Fig. 8) reaches a maximum around 10¹⁵ cm, then rapidly decreases after peak light, and asymptotes to 5 × 10¹³ cm at late times. Assuming a 10⁶ (5 × 10⁵) M_⊙ black hole, the peak and late-time values correspond to approximately 3000 (6000) and 200 (400) gravitational radii, respectively. This latter value is similar to that inferred for the inner edge of the accretion disk at +10 days, indicating that the UV/optical emission at these epochs (+150 days) is consistent with being produced directly by the accretion disk. Here it is assumed that the inner disk radius does not significantly increase in size between the +10 days and +160 days epochs. We justify this assumption by noting that the accretion disk emission is dominated by the hottest, inner regions of the disk, and there are no obvious accretion related processes that would lead to an approximately order of magnitude increase in the inner disk radius. Similar behaviour has been observed in other TDEs, where the late-time UV emission (in this case meaning several years after peak light) is also found to be consistent with an accretion disk origin (van Velzen et al. 2019). In AT 2020zso, however, the accretion disk emission appears to dominate the UV bands already much earlier, similar to the TDE AT 2018fyk (Wevers et al. 2021), where rapid disk formation was inferred from other spectroscopic emission features (Wevers et al. 2019a).

Emission lines originating in an accretion disk may be collisionally excited rather than through photo-ionisation. This can have a profound effect on the Balmer decrement (Hα/Hβ line ratio), as this depends sensitively on the temperature (see e.g. Fig. 11 in Short et al. 2020). If the temperature in AT2020zso was lower than in AT2018hyz, this may help to explain the weakness or absence of both Hβ and Hγ lines. As noted in Short et al. (2020), the typical blackbody temperatures in TDEs (and also in AT2020zso) are much higher than those inferred from the Balmer decrement, so the accretion disk must be significantly cooler than the blackbody emission. An alternative explanation may be that the blackbody modelling, while empirically a good fit to the data, is not intrinsically related to the observed SED shape. In this case, the inferred blackbody temperatures do not represent physical temperatures of the emitting regions.

Finally, we highlight the peculiar early time light curve evolution, with evidence for a clear break at very early times in the ZTF g- and r-band light curves from a L ∝ t² behaviour to a slower evolution afterwards. Unfortunately, we do not have temperature information at the earliest times before the change in behaviour. If the outflow did not cool significantly initially, a homologously expanding outflow would be consistent with the L ∝ t² evolution. This could imply the presence of an additional source of energy injection to keep the material from cooling at the earliest times. We speculate that this may be provided by the initial debris self-intersection and/or disk formation processes. Once the bulk of this energy is radiated, the further evolution is dominated by cooling as the envelope expands. Verifying such a scenario will necessitate observational constraints on the temperature evolution shortly after disruption in future TDEs. Alternatively, non-spherical expansion may also result in differences from the canonical L ∝ t² evolution before the peak in TDEs.

4.3. Evolution in the context of the accretion disk model

The blackbody emission likely has two components: a reprocessing envelope and an accretion disk.

At phase −14 days, the expanding outflow is very likely reprocessing the X-rays produced at very small scales (whose presence is inferred from the presence of Bowen lines). This outflow provides the dominant contribution to the emission lines before peak, gradually decreasing as the material becomes optically thin. The spectrum is hence dominated by very broad Gaussian-like signatures, the hall-mark sign of TDEs. The accretion disk at this time is weak, either contained within the expanding photosphere or still assembling. For this reason, there are not yet any double-peaked signatures in the spectra.

This evolution is consistent with the evolution of the EW ratio of He II/Hα. Shortly after our observations begin, the Hα line emerges, that is, the EW ratio decreases rapidly. This apparent evolution from H-poor to H-rich is a natural consequence of an expanding reprocessing envelope (Roth & Kasen 2018), where Hα suffers from more self-absorption when the envelope is more compact (while He II photons can escape unimpeded). As the envelope expands and cools, Hα becomes less self-absorbed and its EW increases, while the He II emitting region, located closer to the central ionising source, does not change significantly. This process has been observed in several TDEs to explain the evolution of H-rich to H-poor as the outflowing photosphere contracts after peak light (Nicholl et al. 2019; Charalampopoulos et al. 2022), but here we show that this process is very likely also at work at very early times when the envelope first starts expanding.

Around peak light, when the envelope has expanded sufficiently such that it becomes optically thin, the reprocessing becomes much more inefficient. The blackbody (continuum) emission is now a superposition of both the reprocessing outflowing layer (weakening contribution) and the accretion disk (increasing contribution). The envelope is no longer optically thick, so the spectra are dominated by the accretion disk; they show broad double-peaked emission lines that are now visible due to the large contrast with the host galaxy at peak brightness. If the outflowing photosphere was still partially optically thick at –3 days, it may contribute to the spectrum as a low amplitude, broad Gaussian. We speculate that this could help explain the peculiar outer disk radius evolution. As shown in Fig. 8, this outer radius is very similar to the blackbody radius at that epoch. Given the limited S/N of the spectrum, the contribution of the outflow may be below the level that can be detected during the fitting. When the outflow reaches its maximum extent and becomes completely optically thin at ≈ + 5 days, there is no more contamination in the last X-shooter spectrum at +10 days, and the outer radius as inferred from the modelling reflects the true accretion disk outer edge.

Because the mass fall-back rate scales as a negative power-law with time (and this is what mainly powers the accretion disk emission), after peak the contrast with the host starts to decrease. The continuum emission of the accretion disk remains visible in the UV (even at phases +150 days; see Fig. 1) because of the higher contrast with the host (see the host SED in Fig. 2), whereas the optical (continuum as well as line) emission falls below the host level. As a result, the blackbody UV (continuum) emission remains visible, but the optical spectra no longer shows emission line signatures of the disk.

4.4. A rapidly formed, elliptical accretion disk

The Hα and He II emission line profiles display a prominent asymmetry, contrary to predictions from relativistic, circular accretion disks. Warped disks are also able to produce brighter red-than-blue profiles if the warp preferentially obscures the blue-shifted side of the disk. Our fitting results show that the line profiles can be well reproduced by eccentric, inclined relativistic disk models. We fitted all epochs and lines independently and find consistent values for the main disk parameters despite the significant line profile variability that is observed on approximately two-week timescales. We focus below on the results from Hα modelling, given the potential systematic uncertainties introduced by de-blending the He II region.

The line broadening parameter σ is not expected to strongly influence the line profiles (Eracleous et al. 1995); we find values ∼1000−2000 km s⁻¹ for both Hα and He II. Similarly, the inferred inclinations and orientation angles agree well for all epochs and lines (average i = 85 degrees and average ϕ₀ = 240 degrees). The inner and outer radii are largely consistent with theoretical expectations for TDE disks, predicted to be more compact in nature than, for example, AGN disks. Figure 8 shows that the peak of the inferred (expanding) blackbody radius coincides roughly with the outer extent of the disk at a similar epoch. Similarly, at late times the blackbody radius is of the same order as the disk inner and outer radii. This comparison is somewhat ambiguous, as we compare the radii of an elliptical structure (the accretion disk) with a spherical structure (implicitly assumed when calculating the blackbody radius). We stress that this is an order of magnitude comparison only.

It has been shown, on theoretical grounds, that the initial debris following the tidal disruption of a star is distributed on highly eccentric rings, and that this eccentricity may be long-lived (e.g. Syer & Clarke 1992, 1993; Zanazzi & Ogilvie 2020). Hydrodynamical models have further corroborated the picture where an eccentric, extended disk forms around the time when the mass return rate peaks (Shiokawa et al. 2015; Piran et al. 2015; Krolik et al. 2016). A complicating factor in the identification of such a structure is the presence of an optically thick wide-angle outflow at early times (e.g. Sądowski et al. 2016). While it is theoretically unclear if, and if so how quickly, the debris can shed its orbital energy and form a disk (Guillochon et al. 2014; Bonnerot & Lu 2020), observationally it is now well-established that an accretion disk can form on ∼month timescales (Short et al. 2020; Hung et al. 2020; Cannizzaro et al. 2021; Wevers et al. 2019a, 2021) in line with hydrodynamical simulations. In the case of AT 2018hyz, accretion disk modelling similar to that performed here was used to infer a quasi-circular structure (e ≈ 0.1) around 50 days after the peak of the light curve (Hung et al. 2020). Here, we establish the presence of an elliptical accretion disk around peak light, around one month after disruption, providing further evidence that in spite of the theoretical uncertainties and our lack of understanding of the post-disruption dynamics, an accretion disk can form very quickly.

The compact, elliptical nature of the accretion disk in AT2020zso is in contrast with the majority of the literature, in which it is often assumed that apsidal precession (or some other mechanism) will quickly remove orbital energy from the debris, leading to a nearly circular orbit on the scale of the tidal radius. Instead, our findings suggest a highly eccentric structure with a semi-major axis of ∼100 times the tidal radius, similar to that found by hydrodynamical simulations (e.g. Shiokawa et al. 2015).

AT 2020zso is the first TDE where it can be quantitatively confirmed that the initial debris maintains highly eccentric orbits for a significant amount of time. Due to the lack of observational data before the peak in other TDEs, it remains unclear how often this occurs, that is, if inefficient circularisation is common among TDEs. Krolik et al. (2020) argue that whether or not circularisation is efficient depends sensitively on the pericentre radius of the fatal orbit, with circular accretion disks being a rare occurrence.

Future observations of double-peaked emission lines covering pre- and post-peak phases, where disk signatures can be identified unambiguously throughout the evolution, are necessary to establish the eccentricity evolution of TDE disks in more detail.

4.5. Bowen lines and the TDE unification model

We have confirmed AT2020zso as a TDE with Bowen emission features – and the first with double-peaked Bowen lines. The excitation of these lines requires a strong soft X-ray or extreme ultra-violet (EUV) source, as they are powered through a recombination cascade including He II. However, Leloudas et al. (2019) found that the majority of Bowen-strong TDEs were not detected at X-ray wavelengths. In the TDE unification model of Dai et al. (2018), the properties of Bowen-strong TDEs can be explained if the inclination of the newly formed accretion disk is closer to edge-on than face-on. For Eddington ratios of L_Edd > 0.1, the accretion disk is likely slim rather than thin, leading to an optically thick barrier (potentially aided by an optically thick outflow) that results in strong suppression of X-ray photons for an outside observer.

The presence of double-peaked Bowen lines implies that they are formed very close to the accretion disk surface, most likely in the same region as the He II and Hα emitting regions. With a peak UV/optical Eddington ratio of ≈0.5 for AT2020zso, the disk likely has a slim geometry. Combined with the very high inclination angle (∼85 degrees) this may provide the dense gas that produces the Bowen lines through X-ray irradiation, while at the same time explaining the lack of observed X-ray emission by Swift. This provides the first direct confirmation (albeit for a single source) that the orientation of Bowen-strong TDEs is indeed near edge-on, and the unification model laid out by Dai et al. (2018) and Leloudas et al. (2019) is consistent with these results. The high inclination of the newly formed disk may also help explain the lack of intrinsic (as well as TDE) X-ray emission at early and late times – we derived a lower limit for the column density of 10²² cm⁻² to reconcile the observed X-ray upper limit with the expected AGN X-ray luminosity at late times. Such a column could be provided by a high inclination compact accretion disk.

4.6. Comparison to double-peaked TDEs and AGNs

The elliptical accretion disk model that we have employed has been extensively used in the literature for fitting AGN optical emission lines. Typically, it is applied to the (low-ionisation) double-peaked Balmer emission lines, from which parameters are extracted and analysed. We can therefore compare our own results, obtained from fitting Hα in particular, to the typical values inferred for AGN accretion disks. Figure 15 compares double-peaked TDE spectra at different phases with some double-peaked AGN spectra. Considering only the TDE spectra, it becomes apparent that different disk parameters have very different observational signatures. In particular the inclination of the system with respect to the line of sight can dramatically alter the line widths (AT 2018hyz has an inferred inclination angle of ∼50−60 degrees, whereas both AT 2020zso and PTF–09djl have inclinations > 80 degrees).

Fig. 15. Comparison between the known double-peaked TDE spectra at various epochs and double-peaked AGN sources. The grey horizontal line denotes the continuum level for clarity in the top spectrum. Data for NGC 1097 are from Schimoia et al. (2015), for NGC 4958 from Ricci & Steiner (2019), for Arp 102B from Couto et al. (2013), for AT2018hyz from Short et al. (2020), and for PTF–09djl are from Arcavi et al. (2014).

While the emissivity profile indices (between 2–3), line broadening parameters (1000–2000 km s⁻¹) and the outer radii (1500–10 000R_g) we find are typical of AGN samples (Eracleous et al. 1995; Strateva et al. 2003; Storchi-Bergmann et al. 2017), the inner radii we find are significantly smaller (typically > 500 R_g for AGNs). This is in line with theoretical expectations, which predict that the stellar debris will form an accretion disk with a size of about twice the fatal orbit pericentre (of the order of several tens of gravitational radii). It is also consistent with the much broader emission lines observed in AT 2020zso.

One notable feature of the transient emission is the presence of double-peaked high-ionisation lines, He II and N III λ4640. To our knowledge, these are the first high-ionisation line with an observed geometry similar to Hα in both TDEs and AGNs. In AGNs, double-peaked profiles are observed in low-ionisation lines such as the Balmer lines and sometimes Mg II, but no clear double-peaked profiles have been found yet in the (high-ionisation) UV resonance lines (Eracleous et al. 2009). This may be related to the (typically) much higher optical depth in the high-ionisation lines, which are thought to form in the densest parts of the accretion disk (wind). The likeness of He II to Hα suggests that the optical depth in both lines is similar, and hence they originate from largely the same physical region, as we infer from our accretion disk modelling.

Finally, we find marginal evidence for an additional Gaussian component near systemic velocity, with a width of ∼1000−2000 km s⁻¹, for the epoch at +10 days (Fig. 14). This component is very similar to the geometries found for Seyfert 1 galaxies with double-peaked emission lines (Ho et al. 1997; Schimoia et al. 2017; Storchi-Bergmann et al. 2017); it can be interpreted as originating in clouds outside of the accretion disk or outside of the disk plane (e.g. produced from a slow accretion disk wind), and hence they have lower velocity widths and are found near rest velocity. This component appears variable, as it does not appear prominently in the spectra around peak light whereas the triple peaked structure is very clear at +16 days. The difference of this component between the Hα and He II spectra is likely a result of a degeneracy with the additional Gaussian model, leading to a slightly different inferred orientation angle (this component prominently appears for ϕ > 230 degrees). In the interpretation as a disk wind, this variability could be intrinsic: if the wind is initially optically thin (and hence appears as a weak contribution to the overall flux) but over time intensifies (e.g. if the disk reaches a steady state) and becomes optically thick, it will become stronger over time. Alternatively, the wind component may remain steady, but because it is superposed onto the accretion disk contribution (which is observed to diminish over time) it appears variable, becoming more prominent as the accretion disk flux decreases. Our spectra are not of sufficient quality (in terms of S/N) to distinguish between these scenarios.

Higher S/N, medium-resolution spectroscopy (R ≳ 5000) of future double-peaked TDEs can help shed more light on the presence and evolution of this component.

4.7. Constraints on the precession and alignment timescales and black hole spin

On the one hand, the late-time optical spectrum taken with X-shooter (215 days after peak) shows that there is no evidence for persistent broad emission lines, neither in the optical nor in the NIR. On the other hand, the blackbody radii that we derive for the UV/optical emission component is 10¹⁴⁻¹⁵ cm (Fig. 7), well inside typical torus size scales (≳10¹⁶⁻¹⁸ cm, Suganuma et al. 2006; Hickox & Alexander 2018). This suggest that any obscuration must be comparatively low, certainly lower than values observed in Seyfert 2 AGNs (typically in excess of A_V > 5 mag; e.g. Burtscher et al. 2016; Schnorr-Müller et al. 2016). Hence, a type 2 (heavily obscured) AGN configuration is hard to reconcile with the observed photometric and spectral evolution. The implication is that any pre-existing structure, assuming it is aligned perpendicular to the black hole spin vector, has at most a moderate inclination (< 40−60 degrees) with respect to our line of sight.

From our modelling results, we have inferred a very high (∼80 degrees) inclination for the newly formed TDE disk. If it is indeed the case that the torus and any pre-existing accretion disk structures are located in a plane perpendicular to the black hole spin vector, this implies a significant misalignment between the black hole spin axis and the new disk, even after peak light. Such misaligned configurations lead to relativistic torques and so-called Lense-Thirring precession, which will tend to align the newly formed disk with the black hole spin vector as the black hole mass is orders of magnitude larger than the disrupted star.

Our spectroscopic observations can therefore be used to constrain a minimum decay timescale for the misalignment, which depends sensitively on the black hole spin and disk viscosity, as well as weakly on the black hole mass (Franchini et al. 2016; Zanazzi & Lai 2019). Rigid body precession occurs if the local precession period is longer than the sound crossing time, and simulations suggest that this is the case for a large part of parameter space typical for TDEs (Franchini et al. 2016; Zanazzi & Lai 2019). The time between the two observing epochs (∼15 days), for which we infer no significant changes in inclination, eccentricity or disk orientation angle, can therefore be interpreted as an absolute lower limit to the alignment timescale. Ignoring for the moment the (small) effect of black hole mass, with these (very) conservative estimates we can already rule out spin values in excess of a > 0.8 for the central black hole, regardless of disk viscosity (see e.g. Fig. 12 in Franchini et al. 2016).

Longer spectral series may in the future be used to constrain changes in disk inclination. This can, in turn, provide an alternative way to constrain the disk alignment timescales (typically constrained through the detection of X-ray variability such as quasi-periodic oscillations; e.g. Pasham et al. 2019). These alignment timescales then provide a new method (see e.g. Leloudas et al. 2016; Pasham et al. 2019; Mummery & Balbus 2020b) to constrain the black hole spin and/or disk viscosity.

5. Summary

AT 2020zso is a nuclear transient classified as a TDE based on its UV/optical photometric and spectroscopic properties and evolution. We summarise the main results reported in this study as follows:

• The host galaxy has an AGN, based on the emission line content and narrow line region diagnostics (BPT diagram). From the host galaxy velocity dispersion, as well as light curve modelling, we infer a black hole mass of 5−10 × 10⁵ M_⊙.

• Near-infrared photometric observations show no evidence for a light echo ∼180 days after the UV/optical peak.

• The spectra show transient, asymmetric double-peaked line profiles in He II and Hα. These lines significantly evolve from the –14 days before peak light to +14 days post-peak, for which there is spectral coverage, and similar line profiles also become apparent in the He Iλ5876 and N III Bowen lines (λλ4100, 4640). The He II and Bowen profiles are, to our knowledge, the first high-ionisation double-peaked lines observed in an accreting SMBH.

• Light curve modelling indicates that AT 2020zso may have been the result of a partial stellar disruption. This was also the case for AT 2018hyz, another TDE that showed prominent double-peaked Balmer emission lines. This commonality may suggest that the accretion disk was visible directly due to a low amount of obscuring material, compared to more typical full stellar disruptions.

• The low amount of debris resulting from a partial disruption can also help explain the relatively rapid observed blackbody cooling as a consequence of shorter diffusion times in the optically thick envelope. At the same time, this can also explain the rapid spectroscopic evolution: the spectra are initially He-dominated when the envelope is hottest and densest; as it expands and cools, Hα strengthens; and finally, when the debris becomes optically thin near peak light, the double-peaked disk profiles appear.

• Modelling the emission line profiles with an elliptical accretion disk model, we find that the system is highly elliptical and highly inclined (nearly edge-on). From independent fits of He II and Hα at two epochs (−3 and +10 days with respect to peak light), we infer consistent disk parameters such as the inclination and orientation angles, emissivity profiles, line broadening parameters, and inner and outer radii. The consistency between independent fits strengthens the conclusion that the eccentricity in the accretion disk is long-lived.

• The high inferred inclination (a nearly edge-on orientation), combined with the presence of Bowen fluorescence lines and the lack of observed X-ray emission, is consistent with the unification picture of TDEs, where the inclination angle largely determines the observational appearance across wavelengths. To our knowledge, this is the first direct confirmation of this theoretical picture.

• The presence of double-peaked emission lines originating in an accretion disk before peak light confirms that an accretion disk can form very quickly and efficiently (∼1 month after disruption), in contrast to theoretical predictions and simulations. This indicates that our current knowledge of the post-disruption debris is far from complete.

• Around 150 days after the peak, the size of the blackbody radius, as inferred from the light curve, is consistent with the size of the accretion disk inferred from spectroscopic modelling (assuming it is similar to the values inferred at +14 days). This suggests that the UV emission may already be dominated by the accretion disk early on in the evolution.

• We use, for the first time, the lack of change in inclination as inferred from the spectroscopic signatures to constrain the alignment timescale of the newly formed disk with the black hole spin vector. This provides a novel way to probe disk precession in TDEs through spectroscopic monitoring. We find that high black hole spin values (a > 0.8) can be ruled out for the inferred black hole mass.

This work provides a strong link between TDEs and the elliptical accretion disks that are often inferred to explain the asymmetric double-peaked profiles in AGNs. However, the timescales for which these emission lines are visible for both types of sources are vastly different (hundreds of days for TDEs versus decades for AGNs), likely as a result of the rapidly decreasing mass accretion rate following a TDE. Hence, it remains unclear at present if TDEs can also be held responsible for the long-lived asymmetric disk structures observed in some AGNs – a steady influx of material into the elliptical disk is required for the emission lines to remain visible for prolonged periods of time.

Future observations of TDEs with double-peaked line profiles will help shed more light on the post-debris dynamics, including the efficiency of disk formation and subsequent circularisation, as well as on the connection between TDEs and double-peaked AGN sources.

Acknowledgments

We are grateful to T. Hung for sharing a Python implementation of the relativistic accretion disk model, and we thank T. Ricci, J. Schimoia, and G. Couto for providing spectra of NGC 4958, NGC 1097 and Arp102B. We also thank J. Krolik, T. Piran and T. Ryu for insightful comments, and the anonymous referee for comments and suggestions that improved the paper. This work is based on observations collected at the European Southern Observatory under ESO programmes 106.216C, 106.2169.001 and 106.2169.002 (PI: Inserra). The spectra will be made publicly available through WISErep. We thank the Swift team for scheduling the requested ToO observations. The Swift data are publicly available from the Swift science archive. This work is partly based on the NUTS2 programme carried out at the NOT. NUTS2 is funded in part by the Instrument Centre for Danish Astrophysics (IDA). Based on observations made with the Nordic Optical Telescope, owned in collaboration by the University of Turku and Aarhus University, and operated jointly by Aarhus University, the University of Turku and the University of Oslo, representing Denmark, Finland and Norway, the University of Iceland and Stockholm University at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofisica de Canarias. This work makes use of observations from the Las Cumbres Observatory global telescope network. Based on data products created from observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 179.A-2010 and made use of data from the VISTA Hemisphere survey (McMahon et al. 2013). T.M.B. was funded by the CONICYT PFCHA/DOCTORADOBECAS CHILE/2017-72180113. S.Y. is funded through the GREAT research environment grant 2016-06012. P.C. is supported by a research grant (19054) from VILLUM FONDEN. T.-W.C. acknowledges the EU Funding under Marie Skłodowska-Curie grant H2020-MSCA-IF-2018-842471. The LCO team was supported by National Science Foundation (NSF) grants AST-1313484, AST-1911225, and AST-1911151, as well as by National Aeronautics and Space Administration (NASA) grant 80NSSC19kf1639. M.N. acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 948381). N.I. is partially supported by Polish NCN DAINA grant No. 2017/27/L/ST9/03221. S.S. acknowledges support from the G.R.E.A.T research environment, funded by Vetenskapsrådet, the Swedish Research Council, project number 2016-06012. IA is a CIFAR Azrieli Global Scholar in the Gravity and the Extreme Universe Program and acknowledges support from that program, from the ERC under the European Union’s Horizon 2020 research and innovation program (grant agreement number 852097), from the Israel Science Foundation (grant number 2752/19), from the United States – Israel Binational Science Foundation (BSF), and from the Israeli Council for Higher Education Alon Fellowship.

References

Appendix A: Figures

Fig. A.1. Flux-calibrated spectra of AT2020zso (top panel: visible wavelengths, bottom panel: NIR arm of X-shooter). The X-shooter spectra have been resampled to a 0.5 Å dispersion and offset for clarity.

Fig. A.2. Illustration of the appearance of a large third middle bump with increasing eccentricity in the accretion disk model (all other disk parameters are similar to the ones inferred from the fitting). The dashed vertical line marks the rest wavelength, in this case Hα.

Fig. A.3. Comparison of the emission line profiles. Top panel: He II and N III profiles as observed at –3 days with X-shooter. The Ca H+K doublet is interpolated to avoid biasing the subtraction. Bottom panel: Subtracted profile (He II – N III) in blue, and Hα (black) for comparison.

Fig. A.4. Same as Fig. 10 but for the He I line at 5876Å.

Fig. A.5. Full posterior distributions of the model parameter values for the light curve fitting with MOSFit.

Fig. A.6. Full posterior distributions of the model parameter values for the Hα emission line (epoch –3 days). The inclination and orientation angle are sampled in units of radians. Note that to estimate parameter values and uncertainties, these samples need to be combined with their associated importance weights (and hence the values and uncertainties can differ from those reported in Table 5).

Fig. A.7. Full posterior distributions of the model parameter values for the Hα emission line (epoch +10 days). The inclination and orientation angle are sampled in units of radians. Note that to estimate parameter values and uncertainties, these samples need to be combined with their associated importance weights (and hence the values and uncertainties can differ from those reported in Table 5).

Fig. A.8. Full posterior distributions of the model parameter values for the (N III-subtracted) He II emission line (epoch –3 days). The inclination and orientation angle are sampled in units of radians. Note that to estimate parameter values and uncertainties, these samples need to be combined with their associated importance weights (and hence the values and uncertainties can differ from those reported in Table 5).

Fig. A.9. Full posterior distributions of the model parameter values for the (N III-subtracted) He II emission line (epoch +10 days). The inclination and orientation angle are sampled in units of radians. Note that to estimate parameter values and uncertainties, these samples need to be combined with their associated importance weights (and hence the values and uncertainties can differ from those reported in Table 5).

Appendix B: Photometry

All Tables

All Figures

	Fig. 1. Host-subtracted light curves of AT 2020zso as observed by Swift, ZTF, ATLAS, and LCO. Vertical lines indicate epochs of spectroscopic observations. As a result of colour evolution, the peak in the optical bands occurs significantly later than in the UV bands. This is consistent with the observed cooling of the blackbody temperature over time.
In the text

	Fig. 2. Host galaxy SED and best-fit template used to synthesise host galaxy magnitudes in the Swift bands. The data are tabulated in Table 2.
In the text

	Fig. 3. Insets of the narrow emission lines observed in AT2020zso.
In the text

	Fig. 4. Narrow emission line ratios as measured from the late-time X-shooter spectrum. The results place the host galaxy firmly into the AGN (Seyfert) part of the diagram. Dividing lines are taken from Kewley et al. (2001) and (2006).
In the text

	Fig. 5. MOSFit model light curves overlaid on the data. The fits are not great, likely due to the relatively rapid temperature evolution and peculiar light curve behaviour (including a break in the rising part), which the model cannot accommodate.
In the text

	Fig. 6. Blackbody fitting of the TDE SED over time. The uncertainties become very large for the last three epochs as two data points do not provide much constraining power for the blackbody model.
In the text

	Fig. 7. Blackbody parameter evolution of AT2020zso over time. Strong cooling is observed during the first part of the light curve.
In the text

Fig. 8. Blackbody radius evolution over time. Top panel: evolution of the blackbody radius, overlaid with the best-fit linear expansion model (orange line). At late times the radius reaches approximately 5 × 1013 cm (but with large error bars), which corresponds to several hundred Rg for MBH = 106 M⊙. Bottom panel: same evolution, but now converted to units of gravitational radius for two different black hole masses. The orange and red stars mark the inferred disk outer and inner radii, respectively (see Sect. 3.6). At late times, the blackbody radii overlap with the estimated accretion disk size (highlighted with the shaded area, assuming that the disk size has not changed significantly from the last epoch that can be modelled).

In the text

Fig. 9. Spectral evolution of AT 2020zso from early to late times. All spectra are continuum normalised. The different emission features in the 4600 Å blend are indicated with vertical lines: solid black lines indicate H Balmer lines, dashed black lines N III, dotted black lines Fe II, red dashed lines He II, and red dotted lines He I. Black spectra were taken with NTT/EFOSC2, blue with VLTX-shooter, orange with LCO/FLOYDS, and magenta with NOT/ALFOSC. The X-shooter spectra are resampled at lower spectral resolution for clarity. The gap around 5200 Å is due to the low S/N of the data at the edges of the UVB and VIS arms, and the gap around 7300 Å is due to the poor correction of molecular bands.

In the text

	Fig. 10. Same as Fig. 9, but now in velocity space. Left panel: the He II region, and the right panel shows the region around Hα.
In the text

Fig. 11. VLT/X-shooter spectra at various epochs, plotted in velocity space for He II (black) and Hα (orange). To guide the eye, we show: the rest velocity as a solid line; the photospheric velocity with a vertical dashed line, which happens to coincide with the expected position for N IIIλ4640; and the position for Hβ with respect to He II as a black dotted line. For display purposes, narrow emission lines have been removed and the spectra are resampled to lower spectral resolution.

In the text

	Fig. 12. Emission line measurements and their evolution over time. Top panel: equivalent width of the He II complex (black circles) and Hα (orange triangles), measured through direct integration from the spectra. Bottom: their ratio (He II/Hα) as a function of time.
In the text

Fig. 13. Comparison of the emission line profiles. Top panel: He II (blue) and N III 4100 Å lines from the X-shooter spectrum at +10 days. Bottom panel: subtracted spectrum (He II – N III, blue) and Hα for comparison (in black). The line profiles are nearly identical, suggesting the presence of double-peaked Bowen N III lines, the first such line profiles seen in a TDE to date. The discrepancy in the red wing (compared to Hα) is likely related to the imperfect normalisation and subtraction of the N III profile.

In the text

Fig. 14. Best-fit accretion disk plus Gaussian models, overlaid on the emission line profiles of He II (left panels) and Hα (right panels). Black represents the data, blue the accretion disk model, and orange the accretion disk plus outflow model. The best-fit Hα models are overplotted in the He II profiles as dashed lines for comparison. The red wing of He II is poorly subtracted, leading to some differences with the Hα profiles, but overall the results are consistent. The model parameters can be found in Table 5.

In the text

Fig. 15. Comparison between the known double-peaked TDE spectra at various epochs and double-peaked AGN sources. The grey horizontal line denotes the continuum level for clarity in the top spectrum. Data for NGC 1097 are from Schimoia et al. (2015), for NGC 4958 from Ricci & Steiner (2019), for Arp 102B from Couto et al. (2013), for AT2018hyz from Short et al. (2020), and for PTF–09djl are from Arcavi et al. (2014).

In the text

	Fig. A.1. Flux-calibrated spectra of AT2020zso (top panel: visible wavelengths, bottom panel: NIR arm of X-shooter). The X-shooter spectra have been resampled to a 0.5 Å dispersion and offset for clarity.
In the text

	Fig. A.2. Illustration of the appearance of a large third middle bump with increasing eccentricity in the accretion disk model (all other disk parameters are similar to the ones inferred from the fitting). The dashed vertical line marks the rest wavelength, in this case Hα.
In the text

	Fig. A.3. Comparison of the emission line profiles. Top panel: He II and N III profiles as observed at –3 days with X-shooter. The Ca H+K doublet is interpolated to avoid biasing the subtraction. Bottom panel: Subtracted profile (He II – N III) in blue, and Hα (black) for comparison.
In the text

	Fig. A.4. Same as Fig. 10 but for the He I line at 5876Å.
In the text

	Fig. A.5. Full posterior distributions of the model parameter values for the light curve fitting with MOSFit.
In the text

	Fig. A.6. Full posterior distributions of the model parameter values for the Hα emission line (epoch –3 days). The inclination and orientation angle are sampled in units of radians. Note that to estimate parameter values and uncertainties, these samples need to be combined with their associated importance weights (and hence the values and uncertainties can differ from those reported in Table 5).
In the text

	Fig. A.7. Full posterior distributions of the model parameter values for the Hα emission line (epoch +10 days). The inclination and orientation angle are sampled in units of radians. Note that to estimate parameter values and uncertainties, these samples need to be combined with their associated importance weights (and hence the values and uncertainties can differ from those reported in Table 5).
In the text

Fig. A.8. Full posterior distributions of the model parameter values for the (N III-subtracted) He II emission line (epoch –3 days). The inclination and orientation angle are sampled in units of radians. Note that to estimate parameter values and uncertainties, these samples need to be combined with their associated importance weights (and hence the values and uncertainties can differ from those reported in Table 5).

In the text

Fig. A.9. Full posterior distributions of the model parameter values for the (N III-subtracted) He II emission line (epoch +10 days). The inclination and orientation angle are sampled in units of radians. Note that to estimate parameter values and uncertainties, these samples need to be combined with their associated importance weights (and hence the values and uncertainties can differ from those reported in Table 5).

In the text