© The Authors 2024

1. Introduction

An active galactic nucleus (AGN) is powered by an accreting supermassive black hole (SMBH) at the center of a galaxy. A common characteristic exhibited by persistently accreting AGNs is stochastic variability, which is present across all wavebands. Such variability is aperiodic, with fluctuations spanning over several orders of magnitude in temporal frequency (e.g., Markowitz & Edelson 2004). However, more extreme flux and spectral variability in X-ray and/or optical fluxes has been observed in both Seyferts and quasars over timescales of months to years (e.g., Gilli et al. 2000; Graham et al. 2020). These drastic changes in the continuum flux are sometimes accompanied by significant changes in the optical spectra, specifically the fluxes of broad components of Balmer emission lines. Most cases are generally thought to be associated with strong variations in the global accretion rate onto the supermassive black hole (e.g., Trippe et al. 2008; Shappee et al. 2014; Denney et al. 2014; LaMassa et al. 2015; Runnoe et al. 2016). The community is still assessing the full range of physical mechanisms underlying the diversity of observed transient behaviors. Mechanisms proposed to explain the transient phenomena generally focus on structural changes in the accretion flow associated with significant changes in the accretion rate. Some models posit cooling or heating fronts, as in Stern et al. (2018), Noda & Done (2018), and Ross et al. (2018). In other cases, impacts from streams associated with tidal disruption events (TDEs) may be relevant (e.g., Ricci et al. 2020). A few open questions associated with these events are how the various structural components of the accretion inflow and outflow respond – in terms of geometry and emission characteristics – to such major changes in accretion rate and if certain components are even present or absent in some cases.

As a review of the standard emission components in Seyfert AGN, the optical/UV emission is thought to originate from a geometrically thin, optically thick accretion disk that is approximated as a multicolor blackbody spectrum (Shakura & Sunyaev 1973; Novikov & Thorne 1973). Due to absorption along the line of sight by the Galaxy as well as by the AGN host galaxy, the peak of the thermal UV/EUV emission from this component is hard to detect. The X-ray spectrum of Seyfert galaxies typically has a cutoff power-law continuum with reflection features, low energy absorption, and often a soft X-ray excess. The primary X-ray power-law continuum is believed to be produced due to the Comptonization of the optical-UV soft seed photons emitted thermally from the accretion disk (“cold phase”) by energetic electrons in a hot compact region called the corona. The corona typically has an electron temperature k_BT_e ∼ 100 keV (“hot phase”) and is optically thin (Haardt & Maraschi 1991; Haardt et al. 1994; Zdziarski et al. 1995, 1996). A part of this primary emission can be reprocessed by interacting with material close to the central black hole in the accretion disk (George & Fabian 1991; Matt et al. 1991) or with the distant dusty torus (e.g., Murphy & Yaqoob 2009), generating a reflection spectrum. The reflection components are composed of a Compton hump peaking at ∼20–30 keV and an Fe Kα fluorescent emission line around 6–7 keV, depending on the ionization state of iron. Additionally, the parsec-scale dusty torus absorbs the emission from the accretion disk and reemits it in the near- and mid-infrared at wavelengths longer than ∼1 μm. Depending on whether the torus intersects the line of sight and depending on the morphology and viewing effects, the torus can obscure X-rays and give rise to silicon emission or absorption features in the mid-IR (e.g., Hao et al. 2007, and references therein).

In the absence of torus absorption of the X-ray continuum, soft X-ray excess emission is typically observed below 1–2 keV (Arnaud et al. 1985; Singh et al. 1985). The soft excess has been observed in more than 50% of Seyfert 1 galaxies (Halpern 1984; Turner & Pounds 1989). The origin of the soft excess is still not completely understood. A blackbody is generally excluded, as its characteristic temperature has been found to be independent of the black hole mass (Gierliński & Done 2004; Bianchi et al. 2009). Preferred explanations include (i) a “warm Comptonization” region, wherein the seed disk photons are up-scattered in an optically thick (τ ∼ 10 − 40) and warm (k_BT_e ≲ 1 keV) plasma (e.g., Porquet et al. 2004; Mehdipour et al. 2011; Petrucci et al. 2018), and (ii) blurred ionized reflection, wherein the soft excess is produced due to relativistic blurred reflection of the primary X-ray continuum in an ionized disk (e.g., Crummy et al. 2006).

The X-ray and optical emission from AGNs is highly variable, not only in flux but also in broadband spectral shape. Multiwavelength monitoring, including searches for inter-band lags or leads, can yield insight into the nature, structure, and physical properties of various accretion flow components. For example, observations of optical/UV photons lagging behind the X-rays are consistent with irradiation of the disk by the X-ray emission from the compact corona and subsequent thermal reprocessing, where the lags are consistent with light travel times from corona to disk (e.g., Uttley et al. 2003; McHardy et al. 2014; Edelson et al. 2015; Fausnaugh et al. 2016). An alternate interpretation of these lags focuses on the broad-line region (BLR) emitting reprocessed continuum radiation (Netzer 2022). On the other hand, observations of X-rays lagging the optical/UV emission are sometimes attributed to inwardly propagating local mass accretion rate fluctuations (e.g., Shemmer et al. 2003; Marshall et al. 2008).

The extended ROentgen Survey with an Imaging Telescope Array (eROSITA; Predehl et al. 2021) is the soft X-ray instrument on the Russian Spectrum-Roentgen-Gamma (SRG; Sunyaev et al. 2021) mission. It was developed by a German consortium led by the Max Planck Institute for extraterrestrial physics (MPE). The SRG was launched successfully on 13 July 2019 from the Baikonur cosmodrome, Kazakhstan, and it now orbits the Earth-Sun L2 position. The primary goals of eROSITA are to map the entire sky in X-rays and to detect 50 000–100 000 clusters of galaxies up to a redshift of z ∼ 1.3 in order to track the cosmological evolution of large-scale structures. Additionally, eROSITA yields a sample of roughly 10⁶ AGN. Notably, eROSITA has already found several new highly variable extragalactic nuclear sources, including TDE candidates, changing-look AGNs, and other accretion-driven flaring sources (e.g., Sazonov et al. 2021; Malyali et al. 2021; Scepi et al. 2021; Homan et al. 2023a; Liu et al. 2023).

In this paper, we discuss a flare in the Seyfert galaxy LCRS B040659.9−385922 (eRASSt_J040846.1–385132; henceforth J0408−38) detected by eROSITA that shows significant X-ray flux variations between eROSITA’s successive all-sky scans. We performed a multiwavelength follow-up campaign, making use of space-based telescopes such as the X-ray Multi-Mirror Mission (XMM-Newton) and Neutron star Interior Composition Explorer (NICER) to cover the UV to X-ray wavebands as well as public optical survey data. We determined that the source was flaring in the optical-UV and X-ray continua simultaneously. As discussed below, we find that the spectral shape of the soft excess varied significantly between the X-ray observations taken during the flaring epoch and after the flaring of the source. X-ray spectral fits and broadband spectral energy distribution (SED) fits indicate that the soft excess can be explained using a thermal Comptonization model. In addition, our optical spectral observations of the source reveal that the He IIλ4686 emission line fades in intensity, roughly tracking the optical/UV/X-ray continuum.

The paper is structured as follows. In Sect. 2 we present our multiwavelength follow-up data sets and data reduction techniques. In Sect. 3 we present our X-ray spectral analysis. In Sect. 4 we present the results from modeling the broadband optical-UV to X-ray SED. In Sect. 5 we present the analysis from our optical spectroscopic follow-up observations. We discuss our results in Sect. 6 and summarize our conclusions in Sect. 7.

All magnitudes quoted in this paper are Vega magnitudes. All uncertainties for model parameter values determined from X-ray and broadband SED fits are 90% confidence levels.

2. Source detection, follow-up observations, and data reduction

2.1. eROSITA source detection and data reduction

The position of J0408−38 was scanned by eROSITA during each of the five eROSITA All-Sky Surveys (eRASS; Merloni et al. 2024). Data products were extracted using the eROSITA data analysis software eSASS version eSASSuser 211214 (Brunner et al. 2022). The eRASS tasks evtool and srctool were employed to handle the event data and extract the source products, respectively. Data processing version c020 was used. We combined data from all seven Telescope Modules. Source spectra were extracted using circular regions whose radii were varied in size, depending on the brightness of the source: 61″, 94″, 92″, 134″, and 78″, respectively, for eRASS1–5. Backgrounds were extracted from annular regions, with nearby point sources excised; background inner and outer annuli were 136″–759″, 196″–1168″, 190″–1143″, 266″–1663″, and 162″–969″, respectively.

The source was detected by eROSITA during its ignition phase: there was an increase in the X-ray flux as determined by comparison of successive eRASS scans. The flux increased by a factor of 4.5 over a period of six months from $f_{0.2-5.0}=1.{16}_{-0.18}^{+0.15}\times {10}^{-12}$ erg cm⁻² s⁻¹ in eRASS3 (January 2021) to $f_{0.2-5.0}=5.{19}_{-0.28}^{+0.32}\times {10}^{-12}$ erg cm⁻² s⁻¹ in eRASS4 (August 2021). Considering an earlier flux measurement of $f_{0.2-5.0}=0.{29}_{-0.06}^{+0.10}\times {10}^{-12}$ erg cm⁻² s⁻¹ during eRASS 1 (February 2020), the flux increase was a factor of 18 over 18 months. Then, the flux dropped back toward the lower flux level, decreasing a factor of almost 3 by eRASS5 (January 2022). The fluxes from eRASS along with the other X-ray flux values from our follow-up program are listed in Table 1 and plotted in Fig. 1. The X-ray position of the object is at coordinates (RA, Dec) = (62.19204°, −38.85875°) with a positional uncertainty of ∼1.0″ (includes both statistical and systematic uncertainties); the source is thus named eRASSt_J040846.1−385132. The closest counterpart in the optical and IR bands is the Seyfert galaxy WISEA J040846.10−385131.0 = LCRS B040659.9−385922, with optical coordinates of (RA, Dec) = (62.19208°, −38.85861°), at a redshift of z = 0.0574 (Shectman et al. 1996). The optical to X-ray offset is about $0\stackrel{″}{.}36$. We obtained X-ray follow-up observations using XMM-Newton and NICER during the flaring state and during the declining flux state of the source, as described below.

Fig. 1. Multiwavelength light curves of J0408−38. In the top panel, we plot the 0.2–5 keV X-ray flux from the XMM-Newton observations, the five eRASS scans, and the two NICER observations. In the middle panel, we plot the optical light curves from the ATLAS photometric monitoring, showing data taken in the cyan (c, covering 420–650 nm) and orange (o, 560–820 nm) filters, as well as data from our B-band photometric monitoring at CTIO, shown in purple. In the bottom panel, we plot the MIR light curve from WISE’s W1 and W2 bands. The source was flaring around the same time in multiple wavelengths. The epochs of optical spectroscopy are denoted by red vertical dotted lines.

2.2. XMM-Newton data reduction

The first XMM-Newton observation (ObsID 0862770501) was taken during J0408−28’s flaring state, on 21 August 2021 (XMM1), with a 26.40 ks exposure. By March 2022, the X-ray flux of the source had decreased significantly, back to roughly the same flux levels before the flaring one year earlier. XMM-Newton thus observed the source again on 10 March 2022 (ObsID 0903990901), with a 25.65 ks exposure (XMM2). In XMM1, the EPIC pn and MOS camera were operated in Small Window and Partial Window imaging modes respectively. Both cameras were operated in the Large Window imaging mode during XMM2. The data were processed using XMM-Newton science analysis system (SAS) version 19.1.0, applying the latest calibration files and using HEASOFT (v.6.28). We followed the standard procedure for processing the EPIC pn and MOS data using the epproc and emproc tasks within SAS, respectively. Then we removed the time bins having high background by creating a Good Time Interval (GTI) file using the task tabgtigen. Final good exposure times after screening for XMM1 were 16.8 and 24.2 ks for pn and for each MOS, respectively. Final exposure times for XMM2 were 17.1 and 19.7 ks for pn and each MOS, respectively. The source and background spectra were extracted using a circular region with radius 35″ for both observations. The data do not show any significant pile-up as indicated by the SAS task epatplot. The spectra were rebinned with grppha with a minimum of 20 counts per bin to ensure the use of χ² statistics. For the pn, we extracted single (pattern 0) and double (patterns 1–4) spectra separately (hereafter pn0 and pn14). Redistribution matrices and ancillary response files were generated with the SAS tasks rmfgen and arfgen, respectively.

In both observations, XMM-Newton’s Optical Monitor imaged J0408−38 using the UV M2 filter (effective wavelength 2310 Å). Data were taken in both “image” and “fast” modes and reduced using the XMM_SAS routines omichain and omfchain, respectively. These routines apply all necessary calibrations (including flat-fielding and correcting exposures for dead time), and they perform point-source aperture photometry. We obtained AB magnitudes of 17.464 ± 0.008 for XMM1 and 18.906 ± 0.022 for XMM2, but these values are not yet corrected for Galactic extinction nor for systematic uncertainties. We corrected for Galactic extinction using E(B − V) = 0.0081 at the sky position of J0408−38 (Schlegel et al. 1998), and use R ≡ A_V/E(B − V) = 3.1, which yields A_V = 0.025 mag. We use the Galactic extinction curve of Cardelli et al. (1989) to estimate the extinction near M2 as A_M2 = 2.1 × A_V = 0.053 mag. The corrected AB magnitudes are thus 17.411 ± 0.008 for XMM1 and 18.853 ± 0.022 for XMM2. The corresponding extinction-corrected flux densities are (22.125 ± 0.159)×10⁻¹⁶ erg cm⁻² s⁻¹ Å⁻¹ (0.394 ± 0.003 mJy) for XMM1, and (5.861 ± 0.114)×10⁻¹⁶ erg cm⁻² s⁻¹ Å⁻¹ (0.104 ± 0.002 mJy) for XMM2.

2.3. NICER data reduction

NICER observed J0408−38 in two campaigns, as summarized in Table 1. The first campaign consisted of ten observations (ObsIDs 4595040101–4595040110), starting 24 August 2021 (three days after the first XMM-Newton observation) and lasting until 21 September 2021, and totaling 14.87 ks of exposure. The second campaign covered four observations (ObsIDs 4595040201–4595040204), spanning 25 to 30 November 2021, and totaling 21.84 ks in exposure. Extraction of source spectra, background spectra, response files, and ancillary response files all followed standard procedures using the NICER data reduction pipeline. Background spectra were based on the 3C50 background estimator. Counts from the two noisy detectors, 14 and 34, were not excluded, although we verified that their inclusion or exclusion did not impact results. To screen out periods of extreme optical loading impacting energies above 0.3 keV, we rejected time intervals when the detector undershoot rate – detector resets caused when incoming optical photons trigger a release of the accumulated charge (Remillard et al. 2022) – exceeded 150 ct s⁻¹. We filtered against periods of high-energy particle rates by screening out time intervals when the overshoot rate exceeded 1.5 ct s⁻¹. We also applied a background hard band count rate cut of 0.5 ct s⁻¹. Good exposure times after screening were 7.42 and 21.58 ks for the August–September and November observations, respectively. In all spectral fits described in Appendix B, we ignore the data below 0.4 keV due to the effects of optical loading affecting the softest energy bins, and we can obtain good model constraints fitting up to 5 keV. Spectra are grouped to a minimum of 20 counts per bin.

2.4. Optical and IR photometric data

As part of this study, we use the publicly available photometric data from the Asteroid Terrestrial-impact Last Alert System (ATLAS) survey, which were obtained by running forced photometry on the reduced images (Tonry et al. 2018; Smith et al. 2020). The observations from ATLAS are done using the orange (o, 560–820 nm) and cyan (c, 420–650 nm) filters. The resulting photometrically calibrated fluxes are reported in Fig. 1.

We obtained infrared photometry taken using the WISE telescope’s NEOWISE-R project, using the NASA IRSA archives. The WISE bands W1 and W2 are centered on 3.368 and 4.618 microns, respectively. We use the fluxes generated by the automated forced photometry pipeline. We binned the data to a six-month timescale, averaging both fluxes and uncertainties. As shown in Fig. 1, the optical and IR light curves show qualitatively similar variability trends to the X-ray data. The object exhibits a sudden increase and then decrease in flux at all continuum wavelengths probed, with both the increase and decrease lasting approximately six months.

We also obtained 14 B-band observations of J0408−38 at the 0.4-m PROMPT6 telescope at Cerro Tololo Inter-American Observatory (CTIO), operated as part of Skynet, during MJD 59439–59694. All images were reduced following standard procedures, which included bias correction and flat-fielding. Aperture photometry followed standard procedures using circular extraction regions for the source and source-centered annular regions for the background. The resulting Vega B-band magnitudes were converted into mJy are plotted in Fig. 1; they are not corrected for Galactic reddening.

2.5. Optical spectroscopic observations

We obtained a total of five optical spectra of the counterpart of J0408−38 as part of our multiwavelength follow-up campaign. A single spectrum was taken on 23 August 2021 (MJD 59449), during the flaring phase, using the FORS2 instrument (Appenzeller et al. 1998), installed on the UT1 telescope of the ESO 8.2 m Very Large Telescope in Chile. The target was observed with a 1.3-arc-second slit, using three grism+filter configurations: G14OOV (spectral resolution R = 2100; 1000 s exposure), G300V + GG435+81 (R = 440, 600 s exposure), and G300I + OG590+32 (R = 660, 400 s exposure). We combined these exposures to produce a single high-quality spectrum.

Three spectra were taken with the SpUpNIC spectrograph (Crause et al. 2019) at the 1.9 m telescope, also located at the SAAO. One observation occurred on 11 September 2021 (MJD 59468), during the flaring phase. Two additional observations were made as the flare faded, on 11 April 2022 (MJD 59680) and 6 September 2022 (MJD 59828). We made use of grating 7 at an angle of 16.5°, resulting in wavelength coverage of 3500–8500 Å and an wavelength resolution of R = 500. Each observation had a total integration time of approximately 2400 s. All spectra were reduced using standard bias corrections, flat-fielding, and wavelength calibration using arc-lamp spectra taken during the same night. Spectrophotometric calibrations were performed using a standard star taken on the night in the case of the VLT and SAAO 1.9 m observations, and a standard star taken within the past year in the case of the SALT spectrum.

Finally, a further spectrum was taken during the declining phase of the flare on 6 February 2022 (MJD 59616) using the Robert Stobie Spectrograph (RSS; Burgh et al. 2003; Kobulnicky et al. 2003) mounted on the 10 m Southern African Large Telescope (SALT; Buckley et al. 2006), situated at the South African Astronomical Observatory (SAAO) in the Northern Cape province of South Africa. The observation used the pg0900 grating. We used two camera setups. One spectrum was taken using a grating angle of 15.875° (camera angle 31.75°), targeting the Hα to Hγ region of the spectrum, and achieving ∼850 < R < ∼ 1250. The integrated exposure time (excluding overheads) was 600 s. The second setup used a grating angle of 12.875° (camera angle 25.75°), targeting the Hβ region and the blue continuum toward roughly 3500 Å (rest frame), and had an integrated exposure time (excluding overheads) of 540 s. This spectrum had ∼700 < R < ∼ 1000. We used the 2 × 2 spatial binning, the faint gain, and the slow readout mode for both setups.

2.6. Continuum variability overview

Here, we briefly discuss the continuum light curve overview, presented in Fig. 1. Section 3 covers the fits to the XMM-Newton data, which have the highest signal/noise. Based on the results of these fits, we additionally conducted spectral fits for the eROSITA and NICER data to extract the 0.2–5.0 keV fluxes, as illustrated in Fig. 1. Comprehensive details of these spectral fits can be found in Appendices A and B.

The X-ray and optical continua display roughly consistent behavior: both show mild increases from early 2020 through early 2021. Both bands experience data gaps, for example, due to Sunblock for ATLAS, from March 2021 through August 2021. In August 2021, both bands are sampled at the highest fluxes we measure (including eRASS4); unfortunately we cannot be certain precisely as to when the true intrinsic peak in any monitored band – X-ray, UV, optical, or IR – occurred. Then, both optical and X-ray, along with the UV M2 band, all decrease through early 2022, although we caution that X-ray and UV M2 lack good time sampling here. Both ATLAS bands are relatively monotonic in their decreases.

The X-ray and optical bands are consistent with zero lag (as per visual inspection). Due to the data gaps in both light curves, we cannot establish if any non-zero inter-band lag up to a timescale of very roughly six months exists.

3. X-ray data analysis

The spectral analysis for the X-ray observations are performed using the XSpec (v.12.11.1) package (Arnaud 1996); we use the χ² statistics for the model fitting. The uncertainties on the best-fit parameters are estimated at the 90% confidence level (Δχ² = 2.71 for one free parameter) as derived from Monte Carlo Markov Chains (MCMC) using the Goodman–Weare algorithm. We model Galactic absorption with the total hydrogen column density fixed at N_H = 1.27 × 10²⁰ cm⁻² (Willingale et al. 2013), using the photoelectric cross sections of Verner et al. (1996) and solar abundances from Wilms et al. (2000).

3.1. XMM-Newton observations

We perform combined spectral fitting to the pn0, pn14, MOS1, and MOS2 data for both observations. We use a redshift value z = 0.057. For both observations, we use the energy range of 0.25–10 keV for pn0 and 0.2–10 keV for MOS1 and MOS2. Since the double-events spectrum has calibration issues below 0.4 keV when using the small window mode, we use the energy range of 0.4–10 keV for pn14. In all models, we apply a constant term to account for cross-instrumental systematics. We left its value for pn0 frozen at unity; values for other instruments were always extremely close to unity.

The X-ray primary continuum of an AGN can be roughly described by a simple power-law model. Hence, we first fit the spectra with a simple Galactic-absorbed redshifted power-law model (TBABS × ZPOWERLW) to check for the emission and absorption components in the residual of the fit. This does not give a good fit for either observation, as illustrated by the data-to-model ratios for XMM1 and XMM2 respectively, and suggesting the presence of extra components in the spectra. A soft excess is evident from the positive residuals toward the lower X-ray energies (≲1 keV). One of the best-fitting phenomenological models for both observations was the double power-law model (TBABS × [ZPOWERLW + ZPOWERLW]). In XMM1, the photon index of the power-law fit to the soft X-ray band (SX) was very steep (Γ_SX = 2.79 ± 0.04), and the photon index of the power-law fit to the hard X-ray component (HX) was quite flat (${\mathrm{\Gamma }}_{\mathrm{HX}}=1.{28}_{-0.10}^{+0.08}$). This model had χ²/d.o.f. = 1630.75/1631, with the two model components crossing over at ∼2 keV. In XMM2, we find that the soft X-ray component is still present, but has now hardened, as we can fit this component with a power law of photon index ${\mathrm{\Gamma }}_{\mathrm{SX}}=2.{19}_{-0.06}^{+0.08}$; the 0.5–2.0 keV flux of this component has dropped by a factor of 2.5 compared to during XMM1. The hard X-ray component had ${\mathrm{\Gamma }}_{\mathrm{HX}}=1.{12}_{-0.18}^{+0.15}$. This model had χ²/d.o.f. = 958.66/1014 with the two components crossing over at ∼4 keV. In both spectra, the soft component is very power law-like with no obvious strong curvature and it extends up to roughly 2 and 4 keV in XMM1 and XMM2, respectively. However, it is important to note that it has varied dramatically in both normalization and spectral shape from XMM1 to XMM2. We note that these values of photon index for the soft excess are similar to those for many other Seyfert 1s’ soft excesses when modeled with simple power laws (e.g., Petrucci et al. 2018).

We check for narrow emission from an Fe-Kα line at E = 6.4 keV by adding ZGAUSS to the best-fit model, with Gaussian width σ fixed to 1 eV. We do not detect the line with high significance in either observation; χ²/d.o.f. fell by less than 1/1 for both observations. The upper limit on the equivalent width for the Gaussian component using the model (TBABS × [ZPOWERLW + ZPOWERLW + ZGAUSS]) was found to be 45 eV and 155 eV in XMM1 and XMM2, respectively.

Best-fit model parameters for the double power-law model are given in Table 2. We note that the values of the hard X-ray power-law photon index in both observations are unusually hard. It is possible that this very hard continuum is due to a strong contribution from a Compton reflection continuum. However, given the limited bandpass in which to study the hard component, it is not possible to test such a scenario (although the lack of a very strong Fe Kα emission line would argue against it).

Physically motivated models. We tried both phenomenological and physically motivated models to explain the soft excess. Some of the tested models and the resulting fit values are tabulated in Table 2. We first fit the spectra with a redshifted hard X-ray (HX) power law and a single blackbody for the soft excess. For XMM1, this model yielded a poorer fit compared to the double power-law model, with χ²/d.o.f. = 1788.32/1631, and with worse fit residuals, as illustrated in Fig. 2g. For XMM2, fitting a blackbody also yielded a poorer fit compared to the double power-law model, although the discrepancy was not as bad as for XMM1. Here, χ²/d.o.f. was 967.82/1014, and residuals are plotted in Fig. 2h. We also tried fitting the soft excess with DISKBB and DISKPBB components; both are multiblackbody accretion disk emission with temperature as a function of radius T(r)∝r^p, where p is fixed to 3/4 for DISKBB and is a free parameter for DISKPBB. However, neither component could fit the soft band in a satisfactory way, yielding poor data/model residuals and values of χ²/d.o.f. > 1.2 for both.

Fig. 2. XMM-Newton EPIC spectra and best-fitting models. The left column shows the data and models for the XMM1 observation in August 2021, taken during the flaring of the event. The right column shows the data and models for the XMM2 observation, taken six months later as the flare subsided. Panels a and b show the counts spectra, with the best-fitting double power-law model plotted as a solid line. Panels c and d show the “unfolded” best-fitting double power-law models and data; the individual power-law components are plotted with dashed lines, and the summed model is shown as a solid line. Panels e and f show the residuals to double power-law model fits, which are our best-fitting phenomenological models. Panels g and h show the results from fitting a blackbody component to the soft band plus a hard X-ray power law. Panels i and j show the results from fitting a Comptonization model (compTT) to the soft band plus a hard X-ray power law.

Next, we tried the Comptonization model COMPTT, described in Titarchuk (1994), for the soft energy band. We fixed the seed photon temperature T₀ to a value of 20 eV, and the plasma temperature k_BT_e to a low energy value of 0.3 keV for both observations. This model produced a good fit to both observations, with values of χ²/d.o.f. ∼ 1.00 (XMM1) or 0.93 (XMM2) and data/model residuals virtually identical to those for the double power-law fit. The optical depths were τ = 12.3 ± 0.3 (XMM1) and $\tau =14.9_{-1,1}^{+0.9}$ (XMM2), and with hard X-ray photon indices Γ_HX = 1.76 ± 0.04 (XMM1) and ${\mathrm{\Gamma }}_{\mathrm{HX}}=1.{67}_{-0.05}^{+0.06}$ (XMM2).

We also tested if the soft excess could be attributed to relativistically blurred emission from the ionized, inner disk region. We used RELXILL, with emissivity index frozen to 3, iron abundance frozen to solar, outer radius frozen to 400 R_g, power-law cutoff frozen to 300 keV, and inclination angle frozen to 45°. We obtained reasonable fits for both observations: For XMM1 (XMM2), we obtained Γ_HX ∼ 1.45 (∼1.17) and ionization values of log(ξ[erg cm s⁻¹]) = 4.7 (3.7), and values of χ²/d.o.f. ∼ 0.99 (XMM1), 0.94 (XMM2). In both cases, black hole spin a^* pegged at 0.998, as a response to the extreme smoothness of the observed soft excess. However, in the relxill model, relatively higher disk ionization should correspond to a flatter soft excess (García et al. 2013), contrary to what is observed. Consequently, we hold that thermal Comptonization models thus provide the most plausible description of the soft excess.

4. Modeling the optical to X-ray SED

We construct SEDs for both flare peak and post-flare times, close to eRASS4/XMM1 and eRASS5/XMM2, respectively. However, we caution that it is not clear that eRASS4, XMM1, and the start of the ATLAS window in August 2021 sampled the true luminosity peak of the flare, (what we call ‘peak’ is the highest luminosity that we measure using XMM-Newton; the intrinsic luminosity may have been higher and we missed it). Furthermore, it is not clear that the flare had completely concluded by XMM2 and eRASS5. For simplicity and brevity, we henceforth refer to the two SEDs as “peak” and “post-flare”, bearing these caveats in mind.

We construct both SEDs from XMM-Newton EPIC, XMM-Newton OM, and ATLAS data. For ATLAS, we use data taken as close as possible to the two XMM-Newton observations, which occurred on MJD 59447 (peak) and 59656 (post-flare). For the flare peak, we average data taken within ±8 days of XMM1: c-band: 810 ± 18 μJy and o-band: 1150 ± 24 μJy. For the post-flare SED, ATLAS data only go up to MJD 59603 (o-band) and 59586 (c-band). We thus took data from the second half of the observing season (after MJD 59500), performed a linear regression, and extrapolated it to MJD 59656. Uncertainties on flux density are estimated from the most recent (closest in time to MJD 59656) cluster of data points: c-band: 581 ± 17 μJy, o-band: 901 ± 20 μJy. All SED fits were conducted in XSpec.

In both SED fits, we account for the host galaxy starlight contribution by using the spectral template of an Sb galaxy from the SWIRE template library (Polletta et al. 2007) to model it. The Galactic dust reddening is modeled using REDDEN with the value of E(B − V) fixed to 6.86 × 10⁻³ (Schlafly & Finkbeiner 2011)¹. The luminosity distance to the galaxy is 243 Mpc, determined using Ned Wright’s cosmology calculator (Wright 2006)²; throughout this paper, we assume a standard flat cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.286, and Ω_vac = 0.714.

For the multiwavelength SED fitting, we applied the AGNSED model (Kubota & Done 2018), which includes three emission regions: a standard outer disc region, a warm Comptonizing region, and an inner hot Comptonizing region. In the outer disc region, the emission thermalizes locally, forming a blackbody as in a standard disc, having the modified color temperature-corrected blackbody set to unity, but only in the outer regions, over radii extending from R_warm to R_out. Inward, from R_hot to R_warm, the disk is still optically thick, geometrically thin, but produces a warm Comptonized spectrum to explain the steep UV downturn and the soft X-ray upturn (soft excess) observed in AGN spectra (e.g., Davis et al. 2007). The authors consider a truncated disc geometry; interior to R_hot and extending down to the innermost stable circular orbit (R_ISCO) lies the inner hot Comptonization region. The flow does not have an underlying disc component, is geometrically thick, has low optical depth, and produces the primary hard X-ray emission. To model the reprocessing due to the hot corona as an extended source illuminating the warm Comptonization and the outer disc, the model assumes the lamppost geometry to calculate the reprocessed hard X-ray emission by having a spherical hot inner flow from a height H above the black hole, on the spin axis.

We perform a simultaneous fit to both SEDs using the model ZTBABS × REDDEN × (AGNSED + SWS0TEMPLATE). The main parameters of this model are the mass of the SMBH, fixed to M_BH = 2.1 × 10⁷ M_⊙, estimated via the FWHM value of Hβ emission line (discussed in Sect. 5), the distance to the source (fixed to D = 243 Mpc), the Eddington ratio (log(ṁ) ≡ log L_Bol/L_Edd), the dimensionless spin parameter of the black hole (a_*) fixed to zero, the electron temperature of the hot corona k_BT_e, hot fixed at the default value of 100 keV, the electron temperature of the warm corona k_BT_e, warm, and the radial size of the hot and warm coronae, R_hot and R_warm respectively. The log of the outer radius of the disc in units of R_g is set to a negative value (logrout = −1), to use the self-gravity radius as calculated from Laor & Netzer (1989). We fix the hot corona’s scale height to H = 10 R_g (Gardner & Done 2017). The inclination angle of the warm Comptonizing region and the outer disc is fixed at 30°. The only parameters allowed to vary between the two SEDs are log(ṁ), k_BT_e, warm, Γ_hot, R_hot, and Γ_warm; all other parameters are tied between the two SEDs. The model produces best-fitting values of χ²/d.o.f. = 495.10/485 and 339.34/380 for the flare peak and post-flare SEDs respectively. All best-fit parameter values are given in Table 3. From the best-fit parameters, we see that the mass accretion rate ṁ has decreased from ∼0.10 (flare peak) to ∼0.04 (post-flare), that is, ṁ has decreased by a factor of roughly 2.5. Meanwhile, the radius of the hot corona region (R_hot) has increased from the flaring state to the post-flaring state by a factor of roughly 4. The best-fitting spectral index slopes and temperature values for the hot and warm Comptonizing region during the flare observation are: Γ_hot = 1.75 ± 0.06, Γ_warm = 2.71 ± 0.03, and k_BT_e, warm = 0.33 ± 0.04 keV. For the post-flare SED, we obtained ${\mathrm{\Gamma }}_{\mathrm{hot}}=1.{60}_{-0.08}^{+0.06}$, ${\mathrm{\Gamma }}_{\mathrm{warm}}=2.{35}_{-0.05}^{+0.04}$, and k_BT_e, warm = 0.30 ± 0.03 keV. These values of the soft X-ray photon index are consistent with those obtained from the phenomenological double power-law model fits to XMM-Newton data in Sect. 3. These values of the hard X-ray photon indices are consistent with those obtained from the fits to the XMM-Newton data using a model consisting of COMPTT for the soft band a power law for the hard band. The best-fitting AGNSED models are plotted in Fig. 3.

Fig. 3. Flare peak and post-flare SEDs and the best-fitting AGNSED model. In the right panel, we present the (observed and absorbed/reddened) data and the best-fitting models. In the left column, we have corrected both data and models for Galactic absorption and reddening. In both panels, black and purple data points denote, respectively, flaring and post-flaring state data. XMM-Newton pn data have been rebinned by a factor of 4 for plotting purposes only. In both panels, red and blue lines denote, respectively, the best-fitting models to the flaring and post-flaring states. The dashed lines denote the best-fitting AGNSED components, the grey dashed line represents the host galaxy template, and the solid lines represent the total models.

As discussed below in Sect. 5, we assign an uncertainty of 0.45 dex to the estimate of M_BH. We thus refit AGNSED with values of M_BH factors 2.8 higher and lower than above. Values of log(ṁ) are shifted by roughly 0.5–0.6 compared to the 2.1 × 10⁷ M_⊙ case. Meanwhile, photon indices usually shift by roughly 0.1 to 0.3, and k_BT_e, warm usually changes by 0.1 or less. R_hot changes by factors ranging from 1.3 to 2.0 depending on mass. However, for all masses tested, the relative changes in best-fitting model parameters between the high- and low-flux states remain robust: log(ṁ) always decreases by 0.35–0.45; Γ_warm always decreases by roughly 0.4, and R_hot always increases by a factor of roughly 3–4 from the high-flux state to the low flux state.

5. Optical spectroscopic analysis

We display our five optical spectra in Fig. 4; the most visually prominent line emission features in the 4000–7000 Å (rest-frame) range include those associated with Hα, Hβ, Hγ, the [O III] λλ5007,4959 doublet, and He IIλ4686. When one is fitting an AGN spectrum, one typically needs to account for the following continuum components: a blue continuum approximated by a simple power law, Fe II emission, and the host galaxy. We use the empirical template presented in Bruhweiler & Verner (2008) for the Fe II emission, and an Sb spectrum generated by Bruzual & Charlot (2003) for the host galaxy template. However, in all our fits, we found no significant contribution from Fe II emission and thus excluded this component from our fitting process.

Fig. 4. Optical spectra of J0408−38 taken with VLT, SAAO 1.8m, and SALT. The flux has been re-scaled for all the spectra using the [O III] λ5007 emission line. Prominent emission lines are marked. In the inset figure, we show the window range of 4550–4800 Å; spectra are scaled to have matching continuum levels. The He IIλ4686 emission line fades dramatically over the course of the observations.

To account for variations in wavelength resolution and wavelength offsets between different instruments, as well as differing seeing conditions between different observations, we scale (grey-shift) each spectrum by fitting the integrated [O III] λ5007 line flux and wavelength centroid to a reference spectrum following van Groningen & Wanders (1992). Specifically, we use the FORS2 spectrum as the reference spectrum to flux scale all the other spectra since it has the highest S/N. We use the Goodman Weare Markov Chain Monte Carlo algorithm (Goodman & Weare 2010) for fitting. We do not correct for different levels of host galaxy contamination that can result from varying slit dimensions. We de-redden the spectra to account for Galactic extinction following Schlafly & Finkbeiner (2011). We performed the spectral analysis using the lmfit python program, which implements the least-squares fitting for the continuum and phenomenological line fitting. The fitting was performed for three spectral windows: 3800–4400 Å, 4400–5050 Å, and 5200–6800 Å in the rest frame. To fit the emission lines in the spectrum, we fit a local continuum and used Gaussian profiles to model the various emission line features.

We use residuals from the continuum-subtracted spectra to model the emission line features around Hα and Hβ in the wavelength ranges of 6450–6700 Å and 4800–4920 Å, respectively. We fit the narrow [O III] λλ5007,4959 lines. For all other narrow-line components, we fixed their velocity widths to have the same velocity width as [O III] λ5007.

In the detailed spectral fitting for the window range 4400–5050 Å, we fit He IIλ4686, Hβλ4861, and [O III] λλ5007,4959; the spectra and fits are shown in Fig. 5. Table 4 lists our best-fitting model parameters. Of particular note is the flux of the broad He IIλ4686 emission line: it fades dramatically over the course of the observations, concurrent to the fading in the X-ray/optical/UV continua (as shown in Fig. 4). We note that the S/N of spectrum #4, taken at MJD 59679 using the SAAO 1.9 m telescope, was too poor to obtain reliable constraints on Hβ profile parameters, and the He IIλ4686 line was also not detected within the noise limits. More specifically, the He II line flux decreases from (100.1 ± 22.6)×10⁻¹⁶ erg cm⁻² s⁻¹ to < 4.1 × 10⁻¹⁶ erg cm⁻² s⁻¹ from spectrum #1 to spectrum #3 and then to being undetected in spectrum #5, with an upper limit of 1.9 × 10⁻¹⁶ erg cm⁻² s⁻¹.

Fig. 5. The fit residuals after the subtraction of the AGN continuum and host galaxy contribution are represented by the continuous black line. Here, we model the He IIλ4686, Hβλ4861, and [O III] λλ5007,4959 emission profiles. The magenta line marks the fit of the total spectrum including all emission components.

In addition, our fits to spectrum #1 required an extra redshifted broad Gaussian component, as plotted in Fig. 5. The best-fitting energy centroid of this extra broad redshifted component is 4873 ± 3 Å, corresponding to a best-fit velocity redshift of 710 ± 200 km s⁻¹. The width of this component is σ = 28.4 ± 1.6 Å, which corresponds to a FWHM velocity of 4118 ± 232 km s⁻¹.

We consider the possibility that this feature redward of Hβ could be associated with a broad He Iλ4922 “red shelf”, as identified in optical spectra of Seyferts by Véron et al. (2002). As the ionization potential of He I is 24.6 eV, it is possible that such a feature would, like He II, respond to a bluer continuum than that powering Balmer lines and thus exhibit stronger variability than that displayed by Hβ. However, such an identification is tentative at best, due to blending with other lines. Moreover, the best-fit value of the red wing’s centroid, 4873 ± 3 Å, would seem to argue against an identification in He Iλ4922.

Overall, the intensity of the broad Hβ component (excluding the redshifted component detected in spectrum #1) decreases from spectra #1 and #2 to spectrum #5 (we ignore the low-signal/noise spectrum #4 in this statement), concurrent with the X-ray, UV, and optical continuum flux decreases. More specifically, the Hβ intensity decreases by a factor of ∼3 (roughly the same as the UV M2 and X-ray fluxes). However, He II displays a much sharper drop: a factor of ∼30. In the Discussion (Sect. 6), we discuss BLR diagnostics based on these measured line and continuum variability behaviors.

We can obtain virial estimates for both the radial location of the Hβ-emitting region of the BLR, R_BLR, and the black hole mass, M_BH, using the empirical relation between optical luminosity and BLR radius R_BLR in nearby Seyferts (Bentz et al. 2013). From the optical spectral fits, we find an average flux density of λF_5100 Å to be 4.18 × 10⁻¹⁶ erg cm⁻² s⁻¹ Å⁻¹, which for a luminosity distance of 243 Mpc corresponds to a monochromatic luminosity of $\lambda L_{5100\mathrm{\AA }}=2.0\times {10}^{43}$ erg s⁻¹. We use the average value of Hβ FWHM velocity across all spectral fits, 2514 km s⁻¹. We apply the $\lambda L_{5100\mathrm{\AA }}-R_{\mathrm{BLR}}$ relation of Bentz et al. (2013), log(R_BLR, lt-dy) = K + αlog($\lambda L_{5100\mathrm{\AA }}$/(10⁴⁴ erg s⁻¹)), with their best-fitting values of K = 1.56 and α = 0.55. We obtained R_BLR = 3.7 × 10¹⁶ cm = 14.1 lt-days. Assuming a virial factor f of 1.0, we obtained M_BH = fR_BLRv_FWHM²G⁻¹ = 2.1 × 10⁷ M_⊙. We assign an uncertainty of ∼0.13 dex to R_BLR, M_BH, and L_Edd, based on the scatter in the $\lambda L_{5100\mathrm{\AA }}-R_{\mathrm{BLR}}$ relation of Bentz et al. (2013). As a caveat, this use of single-epoch spectroscopy to estimate black hole mass is dependent on the assumption that the Hβ-emitting region is virialized, even during this flaring event. To be conservative, we assigned an additional systemic uncertainty based on the maximum and minimum optical/UV fluxes being a factor of two more or less than the average flux, which yielded an uncertainty of ∼0.17 dex from Bentz et al. (2013) relation. In addition, the scatter in measured values of Hβ FWHM, on the order of ±20% from the average value, translates into an additional uncertainty of ∼0.15 dex. The total uncertainty on M_BH is thus ±∼0.45 dex.

6. Discussions

6.1. Summary of main results

We summarize our key findings from the X-ray data analysis, the modeling of the optical-to-X-ray SED, and the optical spectroscopic analysis of the flare observed in J0408−38 detected by eROSITA below:

• The source was observed to be flaring wherein the soft X-ray flux increased within roughly six months by a factor of roughly 5 (compared to eighteen months earlier, the flux increase was a factor of roughly 18). Concurrently, the total optical continuum emission (AGN + host galaxy) increased by∼25–30%, as observed via ATLAS photometry.

• As the source started to fade again, the X-ray decreased by a factor ∼3 over a period of roughly six months. The UV M2 continuum decreased by a factor of 3.8 during this time. As inferred from SED fits, the optical continuum emission concurrently decreased by roughly 2.5.

• The spectral shape of the soft X-ray excess was found to be strongly varying between the bright and the faint state of the source. XMM-Newton EPIC spectra reveal that the power-law photon index of the soft X-ray band was very steep during the peak of the flare (Γ_SX ∼ 2.8) and flattened considerably by the end of the flare six months later (Γ_SX ∼ 2.2). Thermal Comptonization models are preferred over ionized disk reflection and blackbody models to describe the soft excess.

• We model the broadband optical/UV-X-ray SEDs both near the flare peak, and during post-flare times using the AGNSED thermal Comptonization model to describe the broadband emission, including a warm Comptonized component to describe the soft X-ray excess.

• The SED modeling implies that from near the flare peak to post-flare times, the bolometric luminosity L_Bol/L_Edd dropped by a factor of roughly 2.5, from roughly 10% to 4%.

• From the optical spectral observations, we see that the broad He IIλ4686 emission line was present during the flaring state of the source but faded dramatically (a factor of at least 30) over the next few months.

• We detect a redshifted (by 11 Å) broad component to Hβ line emission at flare peak. It disappeared over the successive two months as the continuum faded.

Based on the observations mentioned above, we explored the possible physical scenarios in the sections that follow.

6.2. Origin of the soft X-ray excess

The soft X-ray excess has been identified in some Seyferts to vary independently and sometimes more slowly than the hard X-ray coronal power law (e.g., Markowitz et al. 2009, NGC 3227; Turner et al. 2001, Ark 564; Edelson et al. 2002, Ton S180). The soft excess can even vary independently of the hard X-ray variability on timescales of years (Mehdipour et al. 2011; Petrucci et al. 2013). Moreover, strong variations in the strength of the soft X-ray excess have been seen, and it has even been observed to disappear between observations (Markowitz & Reeves 2009; Rivers et al. 2012; Noda & Done 2018). Below we discuss the flare in the context of several models, with a focus on the strong spectral variation observed in the soft X-ray excess.

In many Seyferts, the soft excess is modeled well by the “warm” Comptonization model wherein the optical/UV thermal photons are up-scattered in a Comptonizing medium that is optically thicker (τ ∼ 10 − 40) and lower in temperature (k_BT_e ≲ 1 keV) compared to the hot corona, which is responsible for the primary X-ray emission (τ ≲ 1 − 2; k_BT_e ∼ several tens to ∼100 keV). As one typical example, Petrucci et al. (2013) used ten simultaneous XMM-Newton and INTEGRAL observations of Mrk 509 to find that a hot (k_BT_e ∼ 100 keV), optically thin (τ ∼ 0.5) corona is producing the primary hard X-ray continuum and that the soft excess can be modeled well by a warm (k_BT_e ∼ 1 keV), optically thick (τ ∼ 10 − 20) plasma. This model has been successfully applied to many other sources such as NGC 5548 (Magdziarz et al. 1998), RE J1034+396 (Middleton et al. 2009), RX J0136.9−3510 (Jin et al. 2009), Ark 120 (Matt et al. 2014), 1H 0419−577 (Di Gesu et al. 2014), Mkn 530 (Ehler et al. 2018), and Zw229.015 (Tripathi et al. 2019), with warm corona temperatures and optical depths similar to those obtained for J0408−38. The Comptonization model is supported by correlations between optical-UV and soft X-ray variability, as observed for instance by Edelson et al. (1996), Marshall et al. (1997), Kaastra et al. (2011). In the case of J0408−38, the Compton up-scattering of a burst of optical/UV thermal emission from the disk could have produced a corresponding burst of X-ray emission. This could have effectively cooled the warm corona, leading to a steepening in the soft X-ray spectral slope. As the disk emission fades, the corona may not cool as effectively (for example, the optical depth could increase), leading to the soft X-ray slope getting flatter again.

From the best-fitting parameters to our source J0408−38 using the AGNSED model discussed in Sect. 4, we detected a varying soft-excess from the outburst phase to the declining phase with the photon index varying from Γ_warm = 2.71 ± 0.03 to $2.{35}_{-0.05}^{+0.04}$. The temperature of the warm corona was roughly 0.3 keV in both states. It is difficult to state why there would not be major changes in k_BT_warm given the change in seed photon luminosity on strictly physical grounds given that in general, we do not know the physical origin for the warm corona in Seyferts. However, in the context of the geometry of the AGNSED model, then it is possible that the power dissipated in the warm corona remained roughly constant (Kubota & Done 2018).

We estimated the optical depths of the warm and hot Comptonizing regions using their respective best-fit electron temperature and the X-ray power-law photon-index values using the following equation (Zdziarski et al. 1996):

$$\begin{array}{c}\hfill \tau =\sqrt{\frac{9}{4}+\frac{3}{{\theta }_{\mathrm{e}}[{(\mathrm{\Gamma }+\frac{3}{2})}^2-\frac{9}{4}]}}-\frac{3}{2},\end{array}$$(1)

where θ_e = k_BT_e/m_ec² and Γ is the photon index. The optical depth values are both estimated for the warm corona using Eq. (1): (a) peak: τ_hot ∼ 0.5, τ_warm ∼ 16; (b) post-flare: τ_hot ∼ 0.6, τ_warm ∼ 23. Due to poor spectral constraints, the values from our best fits should not be taken at face value, just qualitatively, but the SEDs seem to be consistent with a moderate increase in τ_warm. All of the inferred parameter changes are associated with, and likely caused by, variations in the Eddington ratio, L_Bol/L_Edd, which decreased from 10% during the highest flux state we sampled to 4% over a period of roughly 6 months.

6.3. Characteristic variability timescales

In this section we review the variability timescales associated with a standard Shakura-Sunyaev disk (a reasonable geometry to adopt given the inferred accretion rate); we compare them to the observed timescales to constrain the possible origin of the flaring observed in this source as done in Noda & Done (2018). We calculate the light-crossing time (t_lc), the fastest possible variability timescale for an isotropically emitting region, the orbital timescale (t_orb) associated with Kepler’s third law of motion, the thermal timescale (t_th), which governs heating and cooling, and the viscous timescale (t_vis), which describes the inward radial drift of matter. We used Eq. (1) of Edelson et al. (2019) to estimate the radii of peak flux emission for the optical (600 nm) and UV (220 nm) bands as roughly 55 and 160 R_g, respectively. Here, R_g is the gravitational radius of the black hole given as GM_BH/c², G is the gravitational constant, and c is the speed of light in vacuum. For calculating the timescales, we thus adopted a “typical” radial distance of r ∼ 100 R_g. We assume the viscosity parameter to be α = 0.03 and set H/R (the ratio of the height of the disk to radial distance) to 0.001. The mass of the black hole is M_BH = 2 × 10⁷ M_⊙.

$$\begin{array}{cc}\hfill t_{\mathrm{lt}}& \sim r/c\approx 2.6\mathrm{hours}.\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill t_{\mathrm{orb}}& \sim 2\pi {\left(\frac{r^3}{G{M}_{\mathrm{BH}}}\right)}^{1/2}\approx 7.4\mathrm{days}.\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill t_{\mathrm{th}}& \sim {(2\pi \alpha )}^{-1}\times t_{\mathrm{orb}}\approx 39\mathrm{days}.\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill t_{\mathrm{vis}}& \sim {\left(\frac{H}{R}\right)}^{-2}\times t_{\mathrm{th}}\approx {10}^5\mathrm{years}.\hfill \end{array}$$(5)

On the simplifying assumption that the highest flux point observed represents the maximum flux state of the outburst, we consider the rise and decay times of the outburst to be ∼170 and 190 days, respectively, estimated from Fig. 1. The light-crossing timescale is too short and the viscous timescale for a geometrically thin disk is too long to be comparable to the variability timescale observed in this source. The thermal timescale (alternately, the viscous timescale for a thick disk, wherein H/R approaches 1) best describes the flare. We note that if the value of α was lower than assumed above by a factor of ∼4–5, or if a radius of ∼250 R_g were adopted, then the thermal timescale would much more closely match the observed rise and fall times. We discuss below some of the scenarios or disk instability models below that could possibly cause variability at such timescales.

6.4. Disk instability as the origin of the flare in J0408–38

As noted previously, some accreting black holes display much stronger variability than that associated with the stochastic variability associated with persistent accretion. Instabilities in the accretion disk can potentially produce strong variations in accretion rate, yielding extreme variability in luminosity. We consider several models:

(i) A radiation pressure instability is a candidate for generating flares observed in several Galactic compact sources (e.g., Belloni et al. 1997, GRS 1915+105; Cannizzo 1996, GRO J1744−28; Pahari et al. 2013, IGR 17091−3624), and it is suspected of triggering flares in some Seyfert AGN as well (Grupe et al. 2015, IC 3599; Saxton et al. 2015, NGC 3599; Parker et al. 2019, NGC 1566). The mechanism is outlined in Lightman & Eardley (1974). More recently, simulations were conducted by Śniegowska et al. (2020, 2023), who model an unstable intermediate zone due to radiation a pressure instability operating between the outer cold stable disk and an inner, hot advection-dominated accretion flow.

As a brief overview, the mechanism involves four stages: a quiescent phase, a rising phase, the outburst phase, and the decay phase. Initially, the unstable zone is empty or filled with tenuous gas up to the truncation radius (R_trunc), beyond which the disk is stable. Gas from the outer disk flows into the inner disk, increasing its surface density, but during this stage, the accretion rate and consequently the luminosity do not increase appreciably. The relevant timescale is the viscous timescale for a geometrically thin disk. When the radiation pressure exceeds the gas pressure, a heat wave propagates from the innermost disk out to R_trunc, enhancing the local viscosity, the disk scale height, and local accretion rate. The luminosity increases quickly and so this phase is referred to as the rising phase; the thermal timescale (the time needed to heat the disk) dominates this activity. An outburst phase follows; emission is dominated by thermal emission from the hot and bright disk. The accretion mass starts to drain, and the viscous timescale for a geometrically thick disk dominates the activity. Then, the slow supply of gas from the outer disk fails to compensate for the pace at which the gas in this region is drained. Hence, the luminous phase fades as the disk is emptied and the local accretion rate rapidly decreases. However, the next flare can originate after the inner disk is refilled. Therefore, in this mechanism, recurring flares are possible. The time between successive flares can be approximated using the refilling time of the inner accretion.

In the case of J0408−38, the roughly months-long rise and decay times could be attributed to the rising and decay phases of this instability, in accordance with the thermal timescale (as discussed in Sect. 6.3). However, our observations may have missed the phase of peak emission, and thus we cannot comment on this phase. Additional long-term monitoring will be required to see if additional bursts occur, in which case constraints on the disk refilling time would be attained.

(ii) Another potentially relevant scenario that can cause luminosity flares is the hydrogen-ionization disk instability in the outer accretion disk (Noda & Done 2018). The mass accumulated at the circularization radius is initially quite cold, and the temperature is below that which ionizes hydrogen. However, the surface density and temperature build up over time. Hydrogen gets ionized when a temperature of 3000 K is attained; at this stage, photons get trapped by the free electrons, leading to additional ionization of hydrogen. Once all hydrogen gets ionized, the temperature increases rapidly to 10⁴ K. Such temperatures enhance the sound speed c_s and local accretion rate, causing a heating front to propagate inward.

When the outer disk temperature falls below 10⁴ K, a cooling front propagates that turns ionized hydrogen to neutral. However, the reduction of mass accretion is governed by the viscous time of the outer disc and occurs quite slowly. In the case of AGNs, it is plausible for the H-ionization instability to drive the changing-look AGN phenomenon on timescales of years. However, a caveat is that the H-ionization instability tends to operate on timescales of weeks–months in BHXRBs, and the translation of timescales to AGNs in general and to J0408−38, in particular, is not clear.

These disk instability models are still under development, and it is not immediately clear which models are best applicable to flaring/changing-state AGN in general, or to J0408−38 specifically. Additional observations, including denser monitoring of optical/UV thermal emission to identify the order in which different radial regions of the disk are impacted, may help distinguish among them.

6.5. On the possibility of a tidal disruption event as the trigger to the flare in J0408–38

Stars passing close to supermassive black holes can get tidally disrupted. For a < 10⁸ M_⊙ black hole, the tidal radius lies outside the black hole’s event horizon, and the debris produced from stellar disruption can produce emission, primarily in the soft X-ray and ultraviolet bands (Rees 1988; Drake et al. 2011). Such emissions could decay over timescales of weeks to months. The event rate is estimated to be about ≈10⁻⁵ galaxy⁻¹ yr⁻¹ (Wang & Merritt 2004). TDEs in quiescent galaxies tend to have asymmetric light curves with a fast rise and a slow decay, with optical/UV thermal emission from stellar streams usually peaking months before X-ray emission from circularized material appears. This scenario does not seem to apply to J0408−38, given the simultaneity of the flaring in all bands. However, there are other cases of TDEs occurring in already-active AGNs. Stellar debris can impact the accretion disk, heating it and thus causing excess thermal emission while also supplying extra material for accretion, temporarily increasing L_bol/L_Edd. Shocks in the disk emanating from debris/disk impacts could also temporarily increase the accretion rate (e.g., Chan et al. 2019) and the optical/UV continuum luminosity. For example, the transient PS16dtm is associated with the nucleus of a Seyfert 1 galaxy (Blanchard et al. 2017). This transient is suggested to originate from a tidal disruption event wherein the accretion debris powers the optical/UV emission and obscures the X-ray-emitting region. Homan et al. (2023a) reported on an SMBH transient event detected with eROSITA and Gaia: its optical/X-ray continuum variability properties (fast rise, slow decay) and ultra-soft X-ray spectrum (Γ ∼ 5) categorize it as a TDE. However, the optical spectral properties (broad Balmer lines, strong [O III], [N II], and [S II]) categorize the object as a low-luminosity AGN. A high-redshift example is that of the z = 1.1 quasar SDSS J014124+010306, wherein variability associated with a TDE is claimed to be separated from the AGN continuum variability (Zhang 2022). Finally, in the local Seyfert 1ES 1927+654, the X-ray emission disappeared following major changes in the optical spectral type (Ricci et al. 2020). After the optical and UV outburst, the power-law component produced in the X-ray spectrum vanished, and the blackbody component appeared. Ricci et al. (2020) hypothesized that the event of tidal disruption of a star by the accreting black hole caused depletion of the inner regions of the accretion flow, affecting the magnetic field that powers the X-ray corona. Additional examples of AGN in which TDEs are candidates for fueling temporary flares are Saxton et al. (2015), Grupe et al. (2015), and Parker et al. (2019).

Consequently, we should examine whether the multiwavelength behavior observed in J0408−38 is consistent with TDE-like accretion in an AGN. Speculatively, the transient broad He II emission line and/or the redshifted broad Hβ feature captured during the FORS2 spectrum could have emanated from the tail end of the TDE accretion-like debris falling into the disk. In addition, most optical/UV-selected TDEs have optical/UV luminosities peaking in the approximate range ∼10^{43.5 − 44.5} erg s⁻¹ (van Velzen et al. 2020). In J0408−38, the 3–30 eV luminosity peaks around 2 × 10⁴⁴ erg s⁻¹, which does fall into this range.

However, there are also difficulties in applying this explanation to J0408−38: The slow rise in optical/UV continuum emission in the months before the flare could argue against a TDE origin, although we may not have observed the true luminosity peak due to the data gaps, and that slow rise could also be normal stochastic disk variability.

In addition, our measured He II width during the FORS2 spectrum is only FWHM ∼ 3500 km s⁻¹; TDEs’ He II widths tend to be ∼10⁴ km s⁻¹ and larger (van Velzen et al. 2020).

Finally, we did not observe strong Bowen fluorescence lines, for example, O IIIλ3133 or N IIIλ4640. We conclude that TDE-like accretion in J0408−38 is unlikely.

6.6. Alternative models

There are some alternative models to explain large-amplitude variability at optical to X-ray wavelengths observed in some Seyferts and quasars. Such extreme variability is usually too rapid to be consistent with changes in the global accretion rate associated with the inflow in a standard geometrically thin accretion disc. However, extreme and rapid variations in accretion rate can be associated with a geometrically thick disc supported by magnetic pressure (e.g., Dexter & Begelman 2019). Here, the disk’s thickness is decoupled from the thermal properties, and extreme variations in M_⊙ can occur on ∼year timescales.

6.7. Diagnostics of the broad-line region

6.7.1. SED evolution as inferred from He II and Balmer emission properties

Here, we briefly discuss the potential diagnostics of J0408−38’s BLR from the Hβλ4861 and He IIλ4686 emission lines. We focus on potential reasons why the decrease in observed He II intensity from spectrum #1 to spectrum #5 (factor of at least ∼30 as it fades to becoming undetected by spectrum #5) is so much greater than that for Hβ (factor of ∼3) during the same time.

First, we consider the possibility that the divergence in behavior between the two lines is driven by evolution in the extreme UV (EUV) continuum, since hydrogen and He II emission are driven by continuum emission above 1 Ryd (13.6 eV) and 4 Ryd (54.4 eV), respectively. We hypothesize that as the flare fades, the > 54.4 eV continuum emission fades more sharply than the continuum emission above 13.6 eV, leading to an optical spectrum that becomes increasingly softer and redder. For example, one possibility is that the > 13.6 and > 54.4 eV continua are both dominated by thermal disk emission and that the AGN continuum emission follows a “bluer-when-brighter” and “redder-when-fainter” behavior, as has been suggested to occur in Seyfert AGNs by, for example, Paltani & Courvoisier (1994) and Vanden Berk et al. (2004).

We take the fluxes corresponding to the best-fitting AGNSED models from Sect. 4. The 13.6 eV–10 keV continuum flux drops from 3.72 × 10⁻¹¹ erg cm⁻² s⁻¹ in the high-flux state to 4.54 × 10⁻¹² erg cm⁻² s⁻¹ in the low-flux state, a drop of 8.2. Meanwhile, the 54.4 eV–10 keV continuum drops from 1.76 × 10⁻¹¹ to 3.97 × 10⁻¹² erg cm⁻² s⁻¹, a drop of 4.4. That is, the ratio of > 54.4 eV to > 13.6 eV model flux actually increases as the flare subsides. This increase occurs because as one increases in energy from the optical/near-UV range to the EUV range, in the context of our best-fitting AGNSED fits, the continuum emission becomes dominated by the Comptonized component above ∼10 eV as it connects the EUV regime to the soft X-ray excess. It is thus not clear that the evolution in the EUV continuum slope can be responsible for the diverging behavior in the He II and Hβ lines. However, we must caution that our observations of J0408−38 completely lack rest-frame energy coverage between 2185 Å (5.7 eV) and 66 Å (211 eV), and this exercise is consequently a strongly model-dependent extrapolation.

Moreover, such an exercise is valid only in the context of our AGNSED model fits. An alternate explanation is that, given the direct link between He II flux and the number of ionizing photons in the far-UV part of the SED (e.g., Ferland et al. 2020), it is likely that the ≳54.4 eV continuum did in fact vary by a factor of roughly 30 across our campaign, while the ≳13.6 eV and X-ray continua each only varied by factors of a few. In this scenario, the far-UV-emitting region of the disk experienced a spike in luminosity of at least a factor of 30. Such a statement may be in conflict with the more modest change in L_bol/L_Edd and the far-UV flux indicated by the AGNSED fits. However, our modeling using AGNSED was constrained only by optical, UV (5.4 eV), and soft and hard X-ray constraints; in addition, we may have missed part of the flare due to insufficient data sampling in the time domain. However, finding a suitable broadband emission model that describes SED variations wherein the far-UV varies by a factor of 10 more than the X-ray, mid-UV, and optical bands remains a challenge. One speculative possibility is that a disk instability occurred preferentially concentrated toward the annular region of the disk where the bulk of the far-UV emission is emitted. Following Eq. (1) of Edelson et al. (2019) and Sect. 6.3, this radius is expected to be on the order of a few R_g, that is, very close to the innermost stable circular orbit.

We also consider that the two emission lines are likely probing different radial regions of the BLR, by a radial factor of at least a few. Many works support that both He II emission and Hβ emission arise in regions of mostly virialized, and relatively dust-free gas (e.g., Peterson 2006). However, reverberation-mapping campaigns on several nearby Seyferts have revealed that reverberation lags for He II are generally shorter than those for Hβ, thus supporting radial ionization stratification in the BLR. For example, Kollatschny et al. (2001) and Grier et al. (2012) found He II lags to be shorter than those for Hβ by factors of roughly 6–8 in Mkn 110 and Mkn 335, respectively. In addition to empirical support, photoionization calculations support that the emission-weighted radii of these two emission lines can differ by roughly an order of magnitude (Korista & Goad 2004). In the case of J0408−38, the FWHM velocities for He II and Hβ, averaged across all observations, imply that if both lines originate in gas in Keplerian motion, then He II originates from a radial location roughly one-half that for Hβ. Moreover, the observation of He II emission varying much more strongly than Hβ emission could, speculatively, indicate diverging physical conditions (geometry, density, or illumination) between the two radially distinct emission regions. As a purely speculative example, there could conceivably have existed some temporary, vertically extended structure in the inner, highly ionized region of the BLR only. Such a structure could have created a larger covering factor for the He II-emitting radii than for the larger radii where Hβ is created, as it is most prominent during optical spectrum #1.

Finally, however, a simpler explanation that does not require different covering factors at different radii is based on He II and Hβ emission having divergent “responsivities”. Following Korista & Goad (2004), responsivity is defined qualitatively as how much a given line’s intensity varies for a given variation in ionizing continuum flux. It depends on the spectral shape of the ionization continuum, and on the line-emitting gas’ ionization parameter, density, and radial distance from the central engine (additional details in Ferland et al. 2020). It is defined quantitatively as η ≡ Δlog(F_line)/Δlog(Φ), where F_line denotes line intensity and Φ denotes the number of ionizing photons (Eqs. 1 and 2 in Korista & Goad 2004).

Responsivities generally increase toward larger radii for a given emission line. However, different atomic species can yield vastly different responsivities, depending on local conditions such as density, and on how variable the input continuum is. Consequently, there can exist physical conditions where He II’s responsivity is substantially stronger than that for Hβ. Multiwavelength observations comparing the variability of He II and optical fluxes are somewhat rare. In one recent example, however, multidecade monitoring of Mkn 110 by Homan et al. (2023b) revealed the He II flux to vary by ∼40 while the optical continuum varied by a factor of only a few, similar to variability factors observed during J0408−38’s flare. In addition, Homan et al. (2023b) found line responsivities to evolve between observing epochs, illustrating the extremely dynamic nature of the BLR and the difficulties in finding suitable models.

In our campaign on J0408−38, we found that from spectrum #1 to spectrum #5, the hydrogen-ionizing continuum drops by a factor of roughly 8 (0.9 in the log), considering the 13.6 eV–10 keV continuum emission fluxes from the best-fitting AGNSED model fits to the broadband SED. Meanwhile, Hβ intensity drops by roughly 3 (0.5 in the log), yielding η ∼ 0.5. Comparison to Fig. 4 of Korista & Goad (2004) indicates that such responsivity is typical for BLR number densities (log(n, cm⁻³) ≳ 10). He II intensity drops by at least 30 (1.5 in the log), yielding η ≥ ∼1.7. Such responsivity in He II is reasonably achieved toward relatively higher values of n, for example, log(n, cm⁻³) ∼ 11–12.

6.7.2. On the transient red wing to broad Hβ emission

Finally, we briefly discuss the broad, redshifted emission component detected in the Hβ profile of spectrum #1 only. This feature can be modeled with a Gaussian with an energy centroid suggesting gas moving away from us with a bulk velocity component along the line of sight on the order of 710 km s⁻¹. Such a red wing qualitatively resembles those observed in the Hβ profiles of some nearby reverberation-mapped Seyferts, including Arp 151 (Pancoast et al. 2014) and 3C 120 (Grier et al. 2013). It is not clear if there is any direct connection between J0408−38’s continuum flare and the (apparently transient) presence of this kinematic component, or if the two are merely coincidental. It is also not clear if this red wing arises from a kinematic component that is physically separate and distinct from the gas that yields the (nonshifted) broad emission lines. However, the transient nature of this red wing combined with the relative stability of the other broad Hβ component (“broad component 1”) seems to favor a component associated with bulk radial motion with respect to the rest frame of the BLR. Finally, it is not clear if this feature indicates that the BLR itself contains a kinematic component associated with bulk infall, as measured in several nearby reverberation-mapped Seyferts (Grier et al. 2013; Pancoast et al. 2014). For face-on disks, a redshift can indicate infall. However, for highly inclined disks and for certain viewing angles, a localized outflow driven radially outward could be redshifted.

The broad redshifted kinematic component could, speculatively, be associated with a “Failed Radiatively Accelerated Dusty Outflow” (FRADO; Czerny & Hryniewicz 2011; Naddaf et al. 2021). In this model, the BLR is a dusty wind flowing upward from the disk and pushed radially outward by radiation pressure. However, after the gas is heated by the radiation to the point where dust sublimates, the driving force subsides, and the cloud falls ballistically back down to the accretion disk, forming a “failed” wind. The dynamical simulations of Naddaf et al. (2021) predict that wind heights and velocities depend on source accretion rate, but infall velocities of hundreds to thousands of km s⁻¹ are plausible. In J0408−38, the energy centroid of the broad redshifted component suggests an velocity component along out line of sight of ∼710 km s⁻¹. It is (speculatively) possible that the FORS2 spectrum (#1) caught emission from such a ballistic component, which then faded or impacted into the accretion disk before spectrum #2 was taken, 19 days later.

7. Conclusions

The window into time domain X-ray astronomy of eROSITA enables identification of new transient extragalactic X-ray events associated with major changes in accretion rates in AGNs. When multiwavelength follow-up programs are employed, one can track the evolution of individual X-ray spectral emission components, the optical/UV-emitting accretion disk, and the BLR in response to significant changes in accretion rate and luminosity.

In this work, we have presented a study on a Seyfert 1 AGN, J0408−38 that underwent a flaring event detected using eROSITA. We triggered multiwavelength follow-up observations with XMM-Newton (EPIC and OM) and NICER near the flare peak and during the flare’s fading. We also obtained optical spectroscopic monitoring using VLT-FORS2, SALT-RSS, and SAAO 1.8m-SpUpNIC. The source also flared in the optical and infrared bands, as determined by public photometry from ATLAS and WISE’s NEOWISE-R scans and from photometry taken at CTIO. The X-ray flux increased by a factor of approximately five over six months (and a factor of 18 over 18 months) and had decreased by a factor of about three over the following six months. Meanwhile, the UV flux faded by a factor of almost four from flare peak to post-flare times.

We modeled and compared the X-ray spectra taken during the outburst and during the declining flux state of the source, and we found significant spectral variability in the soft X-ray excess. We argue against simple blackbody and ionized disk reflection models in describing the soft excess. A simple power-law model was found to be the best phenomenological fit for the soft excess. Its photon index flattened from Γ_SX = 2.79 ± 0.04 at the flare peak to ${\mathrm{\Gamma }}_{\mathrm{SX}}=2.{19}_{-0.06}^{+0.08}$ toward the end of the flare. However, modeling the soft excess via a more physical, warm Comptonizing corona also yielded a good fit, with optical depth increasing slightly (for a fixed electron temperature) as the flare faded.

We demonstrated that the broadband SED of J0408−38, from the optical band to the X-rays (and including the soft X-ray excess), could be modeled with the thermal Comptonization model AGNSED such that the soft excess arises in warm Comptonized emission. From the best-fitting parameters, the radial distance of the hot corona region increased from R_hot ∼ 23 to ∼88 R_g from the outburst to the declining flux state. Meanwhile, the optical depth increased slightly from ∼16 to ∼23. In the context of this model, the multiwavelength flaring was the result of a sudden change in the accretion rate ṁ, that is, the accretion rate decreased by a factor of ∼2.5 between the two XMM-Newton observations, as the UV and optical fluxes tracked the disk output. The increase in disk luminosity drove an increase in the luminosities of both the hot and warm X-ray coronae. As the flare faded, the optical depth increased, leading to strong spectral variability in the soft X-ray power-law slope.

Our optical spectroscopic monitoring revealed a strongly dimming broad He IIλ4686 line. As the optical/UV/X-ray continuum fluxes and the flux of the broad Hβ emission line each decreased by factors of a few, the He II flux decreased by a factor of at least 30. Generally, He II is expected to originate closer to the central engine than Hβ, and the diverging flux behavior observed in J0408−38 could, speculatively, indicate diverging physical conditions between the two line-emitting regions – for instance, geometry (covering factor) or illumination. An alternate explanation, however, could be that the He II and Hβ emission have divergent responsivities. A high value of the BLR cloud density could yield a comparatively stronger response in He II than for Hβ for a given variation in the ionization continuum flux.

We also observed a broad redshifted emission component to Hβ in the optical spectrum taken during the flaring state (spectrum #1; VLT-FORS2). Its energy centroid suggests a line of sight bulk velocity on the order of 710 km s⁻¹ moving away from us and may represent a transient kinematic component that is physically distinct from the BLR clouds that produce the broad, nonshifted Balmer line emission.

It is not clear if the bulk flow indicates an inflow or outflow relative to the BLR given the inherent uncertainties in the geometry, viewing angle, etc. However, such a transient component could speculatively be associated with a “failed radiatively accelerated dusty outflow”; (FRADO; Naddaf et al. 2021). In this scenario, the VLT-FORS2 spectrum caught emission from such a ballistic component that then faded or impacted into the accretion disk before spectrum #2 was taken, 19 days later.

As to the underlying mechanism of the flare, the timescales for the increase and decrease of the flare are broadly consistent with those that could be produced by an accretion disk instability, such as the radiation-pressure instability. The timescales for the rising and falling times of the flare are broadly consistent with the thermal timescale near the optical continuum-emitting region of the disk or the viscous timescale for a geometrically thick disk (which approaches the thermal timescale). However, the extreme variation in He II flux could indicate that the far-UV continuum did in fact experience a luminosity spike a factor of at least approximately ten stronger than the optical/UV and X-ray regions of the SED. One speculative possibility is that a disk instability preferentially affected the far-UV emitting region of the accretion disk (near the innermost circular stable orbit).

Future campaigns – either on J0408−38 in the case that another flare occurs, or on other flaring Seyferts – require optical/UV/X-ray photometric monitoring with a sufficiently dense cadence to potentially identify any inter-band lags or leads. Such a result could place constraints on the specific flaring mechanism, disk properties, and the radial direction of propagating fronts.

Acknowledgments

The authors would like to thank Matteo Guainazzi for the thorough reading and for providing valuable comments that helped enhance the manuscript. The authors also thank Chris Done for useful scientific discussions. SK, AM, and TS acknowledge partial support from Narodowe Centrum Nauki (NCN) grants 2016/23/B/ST9/03123 and 2018/31/G/ST9/03224. AM also acknowledges partial support from NCN grant 2019/35/B/ST9/03944. DH acknowledges support from DLR grant FKZ 50 OR 2003. MK is supported by DLR grant FKZ 50 OR 2307. SH is supported by the German Science Foundation (DFG grant number WI 1860/14-1). This work is based on data from eROSITA, the soft X-ray instrument aboard SRG, a joint Russian-German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors and is operated by NPOL with support from the Max Planck Institute for Extraterrestrial Physics (MPE). The development and construction of the eROSITA X-ray instrument were led by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the University of Tübingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of the University of Bonn and the Ludwig Maximilians Universität Munich also participated in the science preparation for eROSITA. The eROSITA data shown here were processed using the eSASS/NRTA software system developed by the German eROSITA consortium. Some of the observations reported in this work were obtained with the Southern African Large Telescope (SALT) under programs 2020-2-MLT-008 and 2021-2-MLT-003. Polish participation in SALT is funded by grant No. MEiN nr 2021/WK/01. This research has made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC. The Skynet Robotic Telescope Network is supported by the National Science Foundation, the Department of Defense, the North Carolina Space Grant Consortium, and the Mount Cuba Astronomical Foundation. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This publication also makes use of data products from NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology, funded by the Planetary Science Division of the National Aeronautics and Space Administration.

References

Appendix A: eROSITA observations

We bin each eRASS spectrum to a minimum of 15 cts bin⁻¹ and modeled them using χ² statistics. All spectra are fit moderately well by simple steep power laws due to the limited photon counts (due to the low exposure times), over the energy range of 0.2−3.0 keV. Results are presented in Table A.1. These spectra are consistent with being completely dominated by the soft excess observed in the XMM-Newton EPIC spectra. In Fig. A.1, we present the data, best-fitting model, and fit residuals. Values of photon indices are consistent with those obtained using XMM-Newton and NICER.

Fig. A.1. Five eRASS spectra, fit with a simple power-law model (top panel), along with the residuals (bottom panel).

Appendix B: NICER observations

We fit the August 2021 (NICER1) and November 2021 (NICER2) spectra separately; the data are shown in Fig. B.1. We first applied a simple power law, but in each case, we obtained values of χ²/dof above 1.2, with broadband concave data/model residuals similar to those seen in the EPIC data when a single power law was applied. For each spectrum, we instead obtained an excellent fit with a double power-law model, qualitatively consistent with the EPIC fits. Results are summarized in Table B.1. We do not include any systematic uncertainties associated with incomplete background modeling; typically, for background systematics of ±2%, Γ_SX, redshifted power-law flux f_{0.5 − 2}, and f_2 − 5 are shifted by up to 0.01, 0.03 × 10⁻¹² erg cm⁻² s⁻¹, and < 0.01×10⁻¹² erg cm⁻² s⁻¹, respectively. Γ_HX is shifted by up to 0.03 (August 2021) or 0.10 (November 2021).

We note that the 0.2–5.0 keV flux measured from NICER1 using the double power-law model fit (χ²/dof = 275.54/281) extrapolated from 0.4 down to 0.2 keV is $f_{0.2-5.0}=5.{83}_{-0.06}^{+0.08}\times {10}^{-12}$ erg cm⁻² s⁻¹. The soft and hard power-law photon indices are ${\mathrm{\Gamma }}_{\mathrm{SX}}=2.{78}_{-0.27}^{+0.48}$ and ${\mathrm{\Gamma }}_{\mathrm{HX}}=1.{49}_{-0.49}^{+0.31}$ respectively. These values are roughly consistent with the fit values from XMM1. In addition, the flux from NICER2 using the double power-law model fit (χ²/dof = 485.88/430) is $f_{0.2-5.0}=3.{99}_{-0.25}^{+0.02}\times {10}^{-12}$ erg cm⁻² s⁻¹ and is roughly consistent with the interpolation of fluxes between XMM1 and XMM2. The soft and hard power-law photon indices during NICER2 are ${\mathrm{\Gamma }}_{\mathrm{SX}}=2.{27}_{-0.12}^{+0.42}$ and ${\mathrm{\Gamma }}_{\mathrm{HX}}=1.{08}_{-0.92}^{+0.63}$, respectively.

Fig. B.1. NICER X-ray spectra of J0408−38 and model fits. The left column shows the data and models for the NICER1 observation in August 2021, taken during the flaring of the event. The right column shows the data and models for the NICER2 observation in November 2021, during the fading of the flare. Panels (a) and (b) show the counts spectra, with the best-fitting double power-law model plotted as a solid line. Panels (c) and (d) show the "unfolded" best-fitting double power-law models and data; the individual power-law components are plotted with dashed lines, and the summed model is shown as a solid line. Panels (e) and (f) show residuals to a simple power-law fit; the fits are poor and suggest the presence of additional components in the spectra. Panels (g) and (h) show residuals from the double power-law model; the fits are good, and in qualitative agreement with fits to the XMM-Newton EPIC spectra.

All Tables

All Figures

Fig. 1. Multiwavelength light curves of J0408−38. In the top panel, we plot the 0.2–5 keV X-ray flux from the XMM-Newton observations, the five eRASS scans, and the two NICER observations. In the middle panel, we plot the optical light curves from the ATLAS photometric monitoring, showing data taken in the cyan (c, covering 420–650 nm) and orange (o, 560–820 nm) filters, as well as data from our B-band photometric monitoring at CTIO, shown in purple. In the bottom panel, we plot the MIR light curve from WISE’s W1 and W2 bands. The source was flaring around the same time in multiple wavelengths. The epochs of optical spectroscopy are denoted by red vertical dotted lines.

In the text

Fig. 2. XMM-Newton EPIC spectra and best-fitting models. The left column shows the data and models for the XMM1 observation in August 2021, taken during the flaring of the event. The right column shows the data and models for the XMM2 observation, taken six months later as the flare subsided. Panels a and b show the counts spectra, with the best-fitting double power-law model plotted as a solid line. Panels c and d show the “unfolded” best-fitting double power-law models and data; the individual power-law components are plotted with dashed lines, and the summed model is shown as a solid line. Panels e and f show the residuals to double power-law model fits, which are our best-fitting phenomenological models. Panels g and h show the results from fitting a blackbody component to the soft band plus a hard X-ray power law. Panels i and j show the results from fitting a Comptonization model (compTT) to the soft band plus a hard X-ray power law.

In the text

Fig. 3. Flare peak and post-flare SEDs and the best-fitting AGNSED model. In the right panel, we present the (observed and absorbed/reddened) data and the best-fitting models. In the left column, we have corrected both data and models for Galactic absorption and reddening. In both panels, black and purple data points denote, respectively, flaring and post-flaring state data. XMM-Newton pn data have been rebinned by a factor of 4 for plotting purposes only. In both panels, red and blue lines denote, respectively, the best-fitting models to the flaring and post-flaring states. The dashed lines denote the best-fitting AGNSED components, the grey dashed line represents the host galaxy template, and the solid lines represent the total models.

In the text

	Fig. 4. Optical spectra of J0408−38 taken with VLT, SAAO 1.8m, and SALT. The flux has been re-scaled for all the spectra using the [O III] λ5007 emission line. Prominent emission lines are marked. In the inset figure, we show the window range of 4550–4800 Å; spectra are scaled to have matching continuum levels. The He IIλ4686 emission line fades dramatically over the course of the observations.
In the text

	Fig. 5. The fit residuals after the subtraction of the AGN continuum and host galaxy contribution are represented by the continuous black line. Here, we model the He IIλ4686, Hβλ4861, and [O III] λλ5007,4959 emission profiles. The magenta line marks the fit of the total spectrum including all emission components.
In the text

	Fig. A.1. Five eRASS spectra, fit with a simple power-law model (top panel), along with the residuals (bottom panel).
In the text

Fig. B.1. NICER X-ray spectra of J0408−38 and model fits. The left column shows the data and models for the NICER1 observation in August 2021, taken during the flaring of the event. The right column shows the data and models for the NICER2 observation in November 2021, during the fading of the flare. Panels (a) and (b) show the counts spectra, with the best-fitting double power-law model plotted as a solid line. Panels (c) and (d) show the "unfolded" best-fitting double power-law models and data; the individual power-law components are plotted with dashed lines, and the summed model is shown as a solid line. Panels (e) and (f) show residuals to a simple power-law fit; the fits are poor and suggest the presence of additional components in the spectra. Panels (g) and (h) show residuals from the double power-law model; the fits are good, and in qualitative agreement with fits to the XMM-Newton EPIC spectra.

In the text