© The Authors 2025

1. Introduction

A fundamental discovery in galaxy evolution is the tight correlation observed between the supermassive black holes (SMBHs) that reside at the center of most galaxies and the properties of the spheroidal component of the host galaxy (Kormendy & Ho 2013). This correlation implies a possible co-evolution between the SMBH and its host. The study of active galactic nuclei (AGNs) can provide pivotal information about this co-evolution, as it reveals the SMBH actively accreting material.

Stochastic variability is a defining feature of AGN (e.g., Kozłowski 2016). Variability has been observed on timescales ranging from minutes to days to decades (e.g., Ulrich et al. 1997; Padovani et al. 2017) in all wavebands. The variability timescales provide a powerful insight into the physical mechanisms beneath them. The minimum timescale observed at a specific band provides an estimate of the linear size of the coherent AGN component emitting at that waveband. Fast X-ray variability is thought to be produced in the innermost regions of the accretion flow. Longer timescales (days to months) are associated with UV and optical variability, which are highly correlated, indicating that they arise from similar physical regions (from annuli of the accretion disk; e.g., Krolik et al. 1991; Arévalo et al. 2008).

An overall active accretion phase can last up to 10^7 − 9 yr (Marconi et al. 2004). However, simulations and observations (Novak et al. 2011; Schawinski et al. 2015) indicate that this phase can be broken down into shorter (∼10⁵ yr) cycles, causing the AGN to “flicker” and implying that SMBHs grow episodically. These timescales are well beyond what can be directly observed in a human lifetime. To fill in this timescale gap, we must look for probes of past activity, potentially imprinted on large physical scales.

Extended emission-line regions (EELRs) can be observed on some AGNs, stretching tens of kiloparsecs from the nucleus. In some, the EELR exhibits evidence of photoionization by a powerful AGN, while the small-scale (10 pc) signatures of activity indicate a nuclear source that is several dex too faint to ionize the large-scale clouds. Considering the light-crossing time from the nucleus, this suggests an AGN fading on timescales of ∼10^4 − 5 yr. Altogether, ∼30 fading AGN candidates have been observed in recent years (e.g., Lintott et al. 2009; Keel et al. 2012; Schawinski et al. 2015; Gagne et al. 2014; Sartori et al. 2016). AGN variability across a wide range of timescales has been uncovered through models (Sartori et al. 2018, 2019) and simulations (Novak et al. 2011). Sartori et al. (2018) presents a framework to test whether variability can be explained based on the distribution of Eddington ratios among galaxy populations. Their analysis proposes that the observed AGN diversity is the result of generating light curves from the same set of statistical properties. Using a compilation of variability measurements from the literature, covering an extensive range of time lags and variability amplitudes, they create a structure function (SF) that provides a valuable overview of the AGN variability on a wide range of timescales.

In this paper, we analyze five nearby active galaxies that present extended emission-line regions. We created ionization models that we compared to our data to estimate the bolometric luminosity needed to explain the photoionization state of the gas at different radii. Considering the light travel time from the central source to the ionized clouds, we can connect this to the AGN luminosity at different moments in the past. With this, we were able to infer the luminosity histories of the nuclear source.

This paper is structured as follows. In Sect. 2, we provide details of the sample analyzed, the Multi Unit Spectroscopic Explorer (MUSE) observations, and the data-reduction process. In Sect. 3, we present our data analysis and results. In Sect. 4, we discuss the reach and implications of our results, and conclusions are drawn in Sect. 5. Standard cosmological parameters (H = 70 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.714, and Ω_m = 0.28) were adopted.

2. Observations and data reduction

Our sample consists of five nearby active galaxies showing EELRs taken from the voorwerpjes sample (VPs; described in Keel et al. 2012). The total VP sample was visually selected from an AGN catalog formed of a combination of optically selected AGNs from SDSS and the Veron & Veron-Cetty (1995) AGN catalog. This sample was complemented by serendipitously found objects from the Galaxy Zoo project (Lintott et al. 2008). The most robust candidates from this sample were followed up with Hα and [OIII]λ5007 imaging, as well as optical long-slit spectra to confirm that the EELRs observed were mainly AGN photoionized.

The five targets analyzed in the present paper were selected for further integral field spectroscopy observations based on: (a) their redshift (z < 0.1); (b) the spatial extension of each galaxy and the associated EELRs, which should fit within the field of view (FOV) of the Multi Unit Spectroscopic Explorer (MUSE) instrument at the Very Large Telescope (VLT); and (c) the surface brightness of the clouds is high enough to achieve a sufficient signal-to-noise ratio in a reasonable exposure time. The galaxies of our sample were observed with the integral field spectrograph MUSE (Bacon et al. 2010) installed on the VLT using the wide field mode (WFM), which covers a 1′×1′ FOV, has a wavelength range between ∼4600 and 9300 Å, resolving power of R = 1770–3590, and spatial sampling of 0$\stackrel{″}{.}$2 per spaxel. The observations were taken as part of program ID 0102.B-0107 (PI: L. Sartori), carried out between March 2019 and April 2020, for a cumulative on-target time of 11 hours and 27 minutes. NGC 5972 and the Teacup galaxies were analyzed in depth by Finlez et al. (2022) and Venturi et al. (2023), respectively. Therefore, we focus solely on processing the remaining sources here. The observations for SDSS J151004.01+074037.1 (hereafter SDSS1510+07), UGC 7342, and SDSS J152412.58+083241.2 (hereafter SDSS 1524+08) were carried out in three observation blocks (OBs), with each OB consisting of three on-target observations and a 950-second exposure. The observations were dithered by 1″ (4 native pixels) in each direction. Data were reduced using the ESO VLT-MUSE pipeline (v2.8 Weilbacher et al. 2020), under the ESO Reflex environment (Freudling et al. 2013). Briefly, the reduction process consists of the following steps. The first is the creation of a master bias and master flat and the processing of the arc-lamp exposures. The calibrations are then applied to a standard star to generate a sensitivity function. Afterward, we applied the bias correction, flux, and wavelength calibrations to the raw science exposures. Given that the sources do not cover the entire FOV, a sky model was created from selected spaxels free of source emission. Then, this model was subtracted from the on-target exposures. To combine the dithered observations of each system, we identified point sources in the FOV to align the individual exposures, and then we generated a stacked final data cube.

3. Results

Color composite images of SDSS1510+07, UGC 7342, and SDSS1524+08 are shown in Fig. 1. These images were created by collapsing three pseudo-bands of the data cubes, encompassing [OIII] (green pseudo-color, bandwidth of 40 Å), Hα (in red, bandwidth of 20 Å), and the stellar continuum (in white) obtained by integrating the emission in the 7000–9000 Å range. In this figure, we can see the morphology of the EELRs. In the case of SDSS1510+07, the ionized gas extends N of the nucleus toward the W, and S of the nucleus towards the E. The clouds appear to be somewhat symmetric, presenting a complex filamentary structure. UGC 7342 shows EELRs extending almost the entire FOV from SE to NW, also showing symmetry. As for SDSS524+08, the ionized cloud is observed mainly in the SE with no NW counterpart and, hence, is highly asymmetric. At least partially, the EELR seems connected to a tidal arm illuminated by the nuclear ionizing source. NGC 5972 shows two large structures extending toward the N and S with a helix-type shape, as discussed extensively by Finlez et al. (2022). The Teacup shows a loop of ionized gas extending 15 kpc to the NE, nicknamed “the handle”. This source has been studied in detail by several authors, including Gagne et al. (2014), Keel et al. (2017), Lansbury et al. (2018), and Venturi et al. (2023).

Fig. 1. Composite color images for targets analyzed in this work with red, white, and green pseudo-colors. Images were created using three bands extracted from the data cube. Green was created with a spectral window that encloses the [OIII]5007 emission lines to represent the highly ionized gas; red was created with a spectral window that encloses Hα 6563 Å, while white represents the stellar continuum, with a spectral window of 7000–9000 Å. In all images, north is up and east is left.

3.1. Stellar populations

In this section, we describe the process we followed to fit and subtract the optical continuum, aimed at isolating the stellar component morphology and kinematics from the ionized gas component. Additionally, we outline the observed properties of the flux distribution, velocity, and velocity-dispersion maps. This process was applied to all the galaxies in the sample.

To ensure a minimum signal-to-noise ratio (S/N) in the entire FOV for our analysis, we bin our data cube using a Python implementation of the Voronoi binning method (vorbin; Cappellari & Copin 2003), requiring a minimum S/N of 20. To study the stellar kinematics and extract the stellar continuum from the data cube for the posterior emission-line analysis, we used the Python version of the penalized spaxel fitting algorithm (pPXF; Cappellari & Emsellem 2004).

A linear combination of template spectra and additive or multiplicative polynomials was used to fit the resulting spectrum for each bin. For the template spectra, we used the E-MILES (Vazdekis et al. 2016) library, covering the entire MUSE wavelength range. The line-of-sight velocity distribution (LOSVD) is parameterized by a Gauss-Hermite function. The outcome of this routine is a model cube for the stellar continuum, which is then subtracted from our original data cube to obtain a continuum-free spectrum for every spaxel.

SDSS 1510+07: In Fig. 2, we show the optical continuum flux by collapsing the data cube around 8000 Å. We carried out an isophotal profile ellipse fitting, from which we infer an ellipticity (ϵ) and position angle (PA) of 0.1 and 50°, respectively, for the central stellar component. A structure reminiscent of a tidal tail is observed to the NW and SE. The velocity and velocity-dispersion maps are consistent with a relatively unperturbed rotating component.

Fig. 2. Flux, velocity, and velocity dispersion maps for the stellar component of SDSS1510+07. The left panel shows the integrated flux, obtained from collapsing the data cube around 8000 Å. Overplotted in white are the isophote profile fits, which show the extension and inclination of a stellar component; this is consistent with a disk. The velocity map (middle panel) shows that the kinematics of the stellar component are rotation-dominated. The velocity-dispersion map (right panel) shows higher velocities concentrated in the center, as expected from a rotation-dominated disk.

UGC 7342: Fig. A.1 shows the flux distribution of the stellar component, which is dominated by a large tidal structure toward the NW, while a fainter structure is also seen to the SE. A component extending 20″ is fit with isophotal ellipses, which indicates an ϵ that changes between 0.2 and 0.1 in the inner 10″ and a PA that increases steadily from 20° to 75°. The velocity and velocity-dispersion maps reveal a small (10″ along PA ∼ 120°) stellar component, reminiscent of a rotating disk, outside which the kinematics are largely non-rotational.

SDSS 1524+08: The flux distribution (Fig. A.2) shows a clearly perturbed stellar component, with tidal arms towards the NE and SE, and a loop-like feature in the SW curving toward the east. The isophotal ellipse fitting in the stellar component shows a PA of 24° and an ϵ of 0.2, changing into 0.1 when reaching 10″ in radius. The kinematics show a highly complex distribution, as evidenced by the velocity and velocity-dispersion maps. There is little evidence of rotational support, as discussed in more detail in Sect. 3.3.

3.2. Ionized gas components

Using our now continuum-subtracted data cube, we fit the emission lines using Gaussian components for every line. We used a custom Python routine based on Bayesian statistics. This routine creates a model spectrum where every emission line is modeled as either one Gaussian or the sum of two components. The model was then fit to the observed spectra (per bin) using an implementation of the Markov chain Monte Carlo (MCMC) adaptive Gibbs sampler. We used three chains to avoid local minima, and sample 10 000 times per spaxel. The following emission lines are fit simultaneously: He II 4685 Å, Hβ 4861 Å, [OIII] 4959,5007 Å, [OI] 6302,6365 Å, [NII] 6548,6583 Å, Hα 6563 Å, and [SII] 6717,6731 Å. [OIII] 5007 Å is hereafter referred to as [OIII]. To reduce the number of degrees of freedom for the fitting process, we tied the following parameters: the fluxes from [NII]6548,6583 and [OIII]4959,5007 with ratios of 2.98:1 and 2.97, respectively, and the width of the Gaussian profiles of Hα and Hβ.

We ran the fit with one Gaussian component for every emission line and then repeated the fit, adding a second component. To identify the regions on every source where a second component is needed while not overfitting the spectrum, we created W80 maps for the [OIII] emission line. This line is used as a probe because it is an isolated bright emission line. This non-parametric descriptor is defined as W80 = V₉₀ − V₁₀, where V₉₀ and V₁₀ are the velocities at which 90% and 10% of the line flux accumulate, respectively. Hence, W80 describes the velocity range that encloses the central 80% of the flux within the spectral window. We added a second Gaussian component in the areas where W80 >  v_lim km/s, where v_lim is 320, 400, and 350 km/s for SDSS1510+07, SDSS1724+08, and UGC 7342, respectively, based on a fit to a 2D histogram of W80 for each galaxy. The contours of this fit above v_lim are shown in the flux map for every galaxy. The second Gaussian component is used to better recover the total flux from every emission line, and not for the kinematical analysis. From these fits, we created flux, velocity, and velocity-dispersion 2D maps.

SDSS1510+07: In Fig. 3, we show the flux, velocity, and velocity-dispersion maps for the [OIII] emission line. The flux distribution is largely dominated by the ionized gas, which seems to follow the stellar tidal tails to the NW and SE. These features are clear in [OIII], and a similar distribution is observed in the other emission lines. The velocity field shows clear rotation within the inner 4″, reaching +162/–243 km/s. This nuclear component is co-rotating with the stellar component, along a different PA. Outside this radius, the velocity field’s major axis seems to shift from the NE to the NW. The tidal tails show redshifted (blueshifted) emission that coincides with the redshifted (blueshifted) side of the disk, reaching projected velocities of –200/+100 km/s. The velocity field for [OIII] shows twisted zero-velocity iso-contours, which can indicate the presence of non-rotating components or an inner kinematically decoupled disk.

Fig. 3. Normalized flux (left), velocity (middle), and velocity-dispersion (right) maps for SDSS1510+07 based on the measurements of the [OIII] emission line. The black contours represent the values (in km/s) of the parameter W80, above the velocity limit used to fit a secondary Gaussian component to the emission-line fit.

UGC 7342: Fig. A.3 shows the flux, velocity, and velocity-dispersion maps for the [OIII] emission line as representative of the other emission lines, which show very similar behavior. The flux distribution shows very extended [OIII] emission from SE to NW, reaching 32″ (33 kpc at the redshift of the galaxy) from the center. The distribution is complex and highly filamentary. The ionized gas does not seem to have the same morphology as the tidal tail observed in the NW of the stellar distribution. Furthermore, there is gas in the opposite direction (SE) of this tidal tail, where we do not observe a strong stellar counterpart. The velocity field shows a NW-to-SE velocity gradient along PA 310°, which could be an indication of rotation. However, it is highly perturbed. We further explore these perturbations in Sect. 3.3.

SDSS1524+08: Fig. A.4 shows the flux, velocity, and velocity-dispersion maps for the [OIII] emission line as being representative of the behavior of the other emission lines. We observe ionized gas mainly in the region of the centrally focused stellar feature and 21″ (16 kpc) to the SE of it, perpendicular to the stellar component extension. The ionized gas coincides with the tidal tail observed to the E, indicating that it may have been swept along the tail and then photoionized. Contrary to the other galaxies in this sample, there is no symmetric counterpart to the EELR in this galaxy. The velocity field shows blueshifted velocities for the EELR. A negative-to-positive velocity gradient is observed along PA 320 degrees. This is perpendicular to the major axis of the stellar component distribution. The velocity dispersion is also larger along this PA. Keel et al. (2024) finds that SDSS1524+08 has a companion galaxy to the SE, along the PA of the ionization cone. This companion shows extended ionized tidal tails, which were photoionized by the nuclear source in SDSS1524+08. This phenomenon is referred to as cross-ionization. The distance from the ionized gas in the companion galaxy to the AGN in SDSS1524+08 suggests longer timescales than the one traced in this work (see Sect. 3.6), reaching ∼10⁵ yr.

3.2.1. Dust extinction

We estimate the total extinction in the V band (A_V) from the Balmer decrement, using the Hα/Hβ ratio. We use the PyNeb (Luridiana et al. 2012) software, assuming an Hα/Hβ = 2.86 intrinsic ratio, with an electron temperature of T_e = 10⁴ K (Osterbrock & Ferland 2006). With this, we derive a reddening map using the extinction curve of Cardelli et al. (1989), for a galactic diffuse ISM (R_V = 3.12). We then corrected the flux maps of our target galaxies for each bin using the derived extinction values.

3.2.2. BPT classification maps

Using the reddening-corrected flux maps for the emission lines, we calculated the spatially resolved [OIII]/Hβ versus [NII]/Hα Baldwin-Phillips-Terevich (BPT; Baldwin et al. 1981) diagrams in order to identify the dominant contribution to the gas ionization on every bin, for every galaxy. This diagram is used to classify the main ionization mechanism into star forming (Kauffmann et al. 2003), Seyfert, LINER (Kewley et al. 2001), or Composite (Schawinski et al. 2007), depending on the position of the spaxel’s emission-line ratios in the diagram. In Fig. 4, we show the BPT classification maps based on the [OIII]/Hβ versus [NII]/Hα BPT diagram for all the galaxies in the sample. All three EELRs indicate that they have been mainly photoionized by the central AGN.

Fig. 4. Spatially resolved BPT classification maps based on [OIII]/Hβ versus [NII]/Hα BPT diagrams for SDSS1510+07, UGC 7342, SDSS1524+08 (top from left to right), NGC 5972, and the Teacup (bottom from left to right). The color bars indicate the different BPT regions, where green indicates the spaxels fall in the Seyfert section of the diagram, red for LINER, blue for composite, and yellow for star forming. Every color is weighted with the [OIII] emission-line flux to more clearly show the bright photoionized regions.

SDSS 1510+07 shows emission classified as composite along PA ∼ 90°, covering most of the centrally concentrated stellar component. In the inner 2″ there is Seyfert-like photoionization along PA ∼ –20°, and outside 5″ we can observe Seyfert photoionization along the tidal arms toward the NW and SE. UGC 7342 only shows Seyfert-like emission along the entire EELR, with some Voronoi bins perpendicular to the ionization cone indicating LINER and composite-like photoionization. In SDSS 1524+08, we observe that in the inner 5″ region – extending in the direction of the stellar component – there is LINER-like photoionization. Some Seyfert-like photoionization is observed along PA ∼ 140° from the center. The entire EELR toward the SE indicates that it is mainly AGN-photoionized. The counterpart toward the NW shows a similar behavior, with a couple of Voronoi bins indicating composite emission.

All galaxies in the sample show that most of the gas is consistent with Seyfert-like ionization. In the cases of SDSS1510+07 and SDSS1524+08, the AGN-like ionization presents evidence for bi-conical structures.

3.3. Kinematics

To analyze the kinematics of these galaxies, we generated position-velocity diagrams (PVDs) using the 3D-Barolo package (Di Teodoro & Fraternali 2015). This code fits a rotation model using rings propagating from the center of the map. Based on the 3D-Barolo 2D fit routine, we fit tilted rings with different PAs, inclinations, and rotation velocities to the stellar component velocity map of each galaxy. With this information, we constructed models of the line-of-sight velocity (V_LOS), which at first order, not considering noncircular motions, can be expressed as

$$\begin{array}{c}\hfill V_{\mathrm{LOS}}(R)=V_{\mathrm{sys}}+V_{\mathrm{rot}}(R)cos\theta sini,\end{array}$$

where V_sys corresponds to the systemic velocity, θ to the azimuthal angle measured in the plane of the galaxy, and i to the galaxy inclination. The disk rotation model is then subtracted from the velocity map in order to reveal noncircular motions. To analyze if the stellar component is supported by rotation, we used the λ_R parameter, which was defined by Emsellem et al. (2007) as a proxy of the observed projected angular momentum of the stellar component and can indicate whether the galaxy kinematics exhibit clear large-scale rotation or not. This parameter is defined as

$$\begin{array}{c}\hfill {\lambda }_R=\frac{{\mathrm{\Sigma }}_{i=1}^N{F}_i{R}_i|V_i|}{{\mathrm{\Sigma }}_{i=1}^N{F}_i{R}_i\sqrt{V_i^2+{\sigma }_i^2}},\end{array}$$

to be measured in 2D spectroscopy. Where F_i, R_i, V_i, and σ_i are the flux, distance to the center, velocity, and velocity dispersion of the stellar component in each ith spaxel (of a total of N in each galaxy), respectively.

SDSS 1510+07: In Fig. 5, we show the rotating disk model fit to the stellar velocity map with the 2D fit routine. The best-fit model shows mean values of 73° and 55° for the PA and inclination, respectively. These values remain almost constant along the different radii. The kinematics are consistent with a rotating stellar disk. For the emission lines, however, we used the 3D fit task, which fits the 3D data cube centered and trimmed to contain only the [OIII] line. For the emission lines, we opted for the 3D fit, considering that information from complex or multiple-component kinematics can be lost in a velocity map. Therefore, a fit directly to the emission line and not just its median velocity can provide more accurate results. We performed this fit for [OIII] and Hα emission lines, using the values obtained from the posited stellar disk as first guesses. In Fig. 6, we show the resulting model subtracted from the corresponding velocity maps. We find that the PA for both lines changes noticeably on every ring after a 0$\stackrel{″}{.}$7 radius. The PA inside 4.5″ radius is fit to be 65°, resulting in an 8° offset from the stellar disk. The continuously changing PA outside this radius indicates that the disk is likely warped. As can be seen from the BPT classification maps (Fig. 4) and velocity maps (Fig. 3), at a 5.5″ radius, we find the EELR toward the NW and SE. The kinematics in these tidal arms appear to be non-rotating. This implies that if they remained in the plane of the galaxy, there would be outflowing gas reaching 220 km/s. The second scenario corresponds to the tidal tail extending outside the plane of the disk; therefore, depending on the inclination angle, this could reach up to ∼490 km/s in deprojected velocity. In the inner 1.2″ we find a component with a blueshift-to-redshift velocity gradient along PA 135°, reaching velocities of 100 km/s. This feature can correspond to a kinematically decoupled component (KDC) rotating on a major axis perpendicular to the rotation of the large-scale disk. The PA offset between the ionized gas and the stellar disk, the tidal arms, and the KDC can all be indicators of past interactions, implying SDSS 1510+07 is likely a post-merger.

Fig. 5. Velocity map for stellar component of SDSS1510+07 (left). Best-fit kinematic Barolo model (middle). Residual map; i.e., velocity map subtracted from the model map (right).

Fig. 6. Velocity maps for emission lines [OIII] (top row) and Hα (bottom row). Velocity map for the ionized gas component of SDSS1510+07 (left). Best-fit kinematic 3D-Barolo model (middle). Residual map; i.e., velocity map subtracted from the model map (right).

UGC 7342: The stellar component shows an inner component that extends over 10″, with two tidal tails extending from it, forming loops toward the NW and the NE (Fig. A.1). The NE tail is fainter and closer to the centrally concentrated component. The main stellar component shows largely non-rotating kinematics. A gradient of ±50 km/s is observed in the inner 4″ along PA ∼ 115° (counterclockwise from N; see Fig. A.1). Outside this nuclear region, there are signs of rotation (blueshift on the SE and redshift on the NW) along PA ∼ 290°; i.e., misaligned 180° from the nuclear component. This could indicate a nuclear disk that is counter-rotating and kinematically decoupled from the remainder of the stellar component. This type of feature can often be observed on disturbed galaxies (e.g., Emsellem et al. 2007); this also seems to be suggested by the presence of two tidal tails. We fit the possible rotation on a disk with up to a 4″ radius, following the same method described above. We obtain mean values of 115° and 35° for the PA and the inclination, respectively (Fig. A.5). The velocity maps for the ionized (Fig. A.3) gas show a highly disturbed gradient from blue to redshifted velocities along a similar PA (310°) to the large-scale stellar distribution. We fit a rotating model to the inner 6″ to cover the most undisturbed section of the gaseous disk (Fig. A.6). We find that part of the observed velocities can be fit by a rotating disk. The PA varies from 300–320° indicating a warping of the disk. The residuals denote an important contribution of non-rotating components.

The highest-velocity residuals are redshifted toward the E and blueshifted toward the W. The spectra in these regions are complex. In the former, we find a broad wing, while in the latter we see a separated blueshifted component, as well as a broad wing. The deviation from the systemic velocity reaches 100 km/s for the W feature and up to ∼400 km/s for the E one. These non-rotating components suggest the presence of an outflow if the SW side of the galaxy corresponds to the near side. Beyond the central kiloparsecs, the ionized gas retains the same sense of rotation along the tidal tails. In the blueshifted region (SE), there are areas with higher velocity dispersion (Fig. A.3) and double-peaked profiles. This implies that kinematics are more complex than pure rotation. We probed this region, extracting flux along different slits. The associated PVDs show the presence of an additional (fainter) component in the areas where the velocity dispersion reaches ∼200 km/s. This component could be interpreted as outflowing gas. There is no observable counterpart to this possible outflow in the redshifted (NW) side of the galaxy; this can be due to more entrained gas or an asymmetric morphology. Another possibility is that it corresponds to the presence of multiple components along the line of sight. Considering that the gas in the tidal tails has likely left the plane of the disk, the velocities observed are harder to interpret. However, the velocity dispersion is mostly around ∼100–150 km/s. Kinematic studies on interacting galaxies with nuclear activity indicate that tidally induced forces can produce dynamically hot systems characterized by velocity dispersions of ∼100 km/s, while systems with σ > 300 km/s show signs of gas flows (e.g., Rich et al. 2015; Bellocchi et al. 2013). Therefore, it is likely that in this source, the tidal features display kinematic turbulence that can be attributed to past interactions.

SDSS 1524+08: We calculated the λ_R parameter along the major axis and found a resulting profile consistent with a slow rotator classification. This denotes that the kinematics of the large-scale galactic distribution are not dominated by rotation. However, in the inner 10″ of the stellar velocity field, there are possible signs of rotation along PA ∼ 20°. The velocity field of the large-scale ionized gas component (Fig. A.4) shows a strong velocity gradient at PA ∼ 321°. This gradient reaches velocities of –220 km/s toward the SE and 360 km/s to the NW. A second gradient can be observed along PA ∼ 20–30°, similar to the stellar component. This gradient reaches velocities of –60 km/s from systemic velocity toward the NE and 78 km/s to the SW.

Considering that the ionization cone aligns with the stronger velocity gradient – as shown by the spatially resolved BPT classification maps – and that the perpendicular velocity gradient aligns with the stellar rotation major axis, it is clear that, on large scales, the stellar and ionized gas components are offset. We consider a scenario in which both the stellar component and the ionized gas remain partially rotating along PA ∼ 22°, in the inner 15″ and 25″ for the stellar and gas components, respectively. To test this scenario, we used the 2D routine of Barolo to fit the velocity maps of the stellar component and the [OIII] emission line. For the latter, we used an [OIII] velocity map derived from a new Gaussian fit, which was carried out in a data cube binned considering the S/N in a spectral window that contains only the [OIII] line. This was done to increase the spatial resolution of this map. We fit the inner (< 10″) of the stellar component leaving PA, V_sys, and V_rot as free parameters (Fig. A.7). The residual map of this shows some agreement with a gradient caused by rotation in a plane along PA ∼ 22°. Considering the possibility that this rotation corresponds to a component of the ionized gas, we used the same PA to fit the velocity map of the [OIII] map, leaving only V_sys and V_rot as free parameters (Fig. A.8). Considering that the gradient along PA ∼ 321° heavily dominates the velocity map, the fits at all radii will be strongly biased by this component. Therefore, the residual map shows there was not a clean subtraction of the gradient along PA 22°. To better visualize this, we extracted the velocities of both components along this PA and fit a rotation curve (Fig. A.9) that shows the rotation gradient is similar in the inner ∼8″, but the ionized gas reaches higher velocities. We interpret the strong gradient along PA ∼ 321° as gas outflowing along the ionization bicone if we consider the NW side as the far side of the galaxy. The velocity dispersion along the bicone axis (Fig. A.4) reaches values of ∼400 km/s. Additionally, the emission-line profiles along this axis reveal the presence of a wing. These features – combined with the high velocities observed – suggest the presence of an outflow, likely along the ionization bicone.

Finally, a feature of bright ionized gas towards the SE also shows high blueshifted velocities. Considering that this feature is mostly disconnected from the main galactic distribution and physically coincides with the tidal arms, we infer that this gas was located along the tidal arm and was ionized by the AGN and likely pushed by the outflow along the bicone. However, given the low velocity dispersion of this feature and the unknown 3D position of the tidal arms, it is not possible to confirm that the gas is indeed outflowing as well. The kinematics of the other two galaxies from this sample, NGC 5972 and the Teacup, were already described in detail by Finlez et al. (2022) and Venturi et al. (2023), respectively.

3.4. Extended emission-line regions

The spatially resolved BPT classification maps established that the EELRs observed in this sample have been mainly ionized by the nuclear AGN emission (see Sect. 3.2.2). It follows that the central source has illuminated and ionized the entire extension of the observed EELRs. Considering the light travel time from the nuclear source to the extension of the EELR, we can trace the luminosity evolution over this time range, thus generating the AGN light curve over timescales of 10^4 − 5 yr. To probe the AGN luminosity at different radii, we compare our observations with photoionization models, following the procedure described in detail by Finlez et al. (2022) and briefly summarized below.

3.5. Description of the photoionization models

These models were constructed using the photoionization code CLOUDY (version 17.02, Ferland et al. 2017). This software solves the thermal and statistical equilibrium equations on a plane-parallel slab of gas that is being ionized by a central source. To model the ionization source, we considered observations of local AGNs (as described in Elvis et al. 1994; Kraemer & Crenshaw 2000) and defined a consistent spectral shape of the central source by a broken power law (L_ν ∝ ν^α), with α = −0.5 for E < 13.6 eV, α = −1.5 for 13.6 eV < E < 0.5 keV, and α = −0.8 for E > 0.5 keV, where E = hν corresponds to the photon energy. The He II/Hβ ratio is sensitive to the shape of the ionizing continuum (e.g., Penston & Fosbury 1978; Kraemer et al. 1994; Gagne et al. 2014). More specifically, it is related to the ionizing UV power-law index α in the range 228–912 Å. In Fig. 7, we show the spatially resolved ratio for our objects. The ratio is mainly constant along the EELRs in each galaxy, albeit different from source to source, and it is lower (suggesting a flatter ionizing spectrum) in the nuclear region of all galaxies. This region is not considered in the following analysis. Considering that the ratio does not change with radius along the EELR, we keep the spectral shape of the continuum fixed for all models, as further discussed in Sect. 4.5.

Fig. 7. Spatially resolved He II/Hβ ratio maps. The ratio remains largely constant along the EELRs for all sources, with lower values in the central region. This area is excluded from our luminosity history analysis. Since this ratio reflects the continuum shape and is mostly constant along the ionized clouds, we assume a fixed spectral shape of the continuum with radius.

The ionization state of the gas depends on the intensity of the ionizing photon flux relative to the gas density. This can be encompassed in the dimensionless ionization parameter, U, which at the illuminated face of the cloud is defined as

$$\begin{array}{c}\hfill U=\frac{Q(H)}{4\pi r_0^{2}n(H)c},\end{array}$$

where r₀ is the distance between the nuclear source and the cloud, n_H is the hydrogen number density in cm⁻³, and c is the speed of light.

The number of ionizing photons per unit of time (s⁻¹) is then defined as

$$\begin{array}{c}\hfill Q(H)={\displaystyle {\int }_{{\nu }_0}^{\infty }}\frac{L_{\nu }}{h\nu }d\nu ,\end{array}$$

with L_ν representing the AGN luminosity as a function of frequency, h the Planck constant, and ν₀ = 13.6 eV/h the frequency corresponding to the hydrogen ionization potential (Osterbrock & Ferland 2006).

We then created a grid of models varying the value of the ionization parameter, while assuming a solar metallicity. Following Keel et al. (2012), we considered a fully ionized gas scenario. Hence, we assumed that the hydrogen density (n_H) is the same as the electron density (n_H = n_e). We then estimated the electron density using the [SII]λ6716/6730 Å emission-line ratio (hereafter referred to as the [SII] ratio) with the help of the PyNeB Python package. This software can compute the emission-line emissivities (Luridiana et al. 2012). We further assumed a temperature of 10⁴ K (Osterbrock & Ferland 2006) and finally obtained an electron density from the [SII] ratio for every aperture used. The models were run with five different values of column density, N(H), ranging from 10¹⁹ cm⁻² to 10²³ cm⁻². We describe the effects of the N(H) changes in Appendix B in more detail. In Sect. 3.6 we show the case with a lower mean chi-square value for all objects, where N(H) = 10²⁰ cm⁻².

In every source, we defined a sequence of apertures, starting from the nucleus and reaching the largest extension of the EELR (shown in Fig. 8). We filled the vicinity of each aperture, centered up to 1.5 times the radius of the aperture, with 30 new apertures of the same size. From these new apertures, we chose the one with the largest S/N to carry out the analysis. We extracted the integrated spectrum of each aperture and fit two Gaussian components to every emission line. The fit emission lines are the following: Hβλ4681, [OIII]λ5007, He IIλ4685, [OI]λ6365, [NII]λ6548,6583, Hαλ6563, and [SII]λ6716,6730. We computed the integrated emission-line fluxes and their ratios relative to Hβλ4861 in order to compare them to the line ratios produced by our photo-ionization models, as CLOUDY delivers the emission line fluxes as ratios relatively to Hβλ4861.

Fig. 8. Apertures used for luminosity history analysis marked as black circles over [SII] emission line normalized flux maps. Left to right: SDSS1510+07, UGC 7342, and SDSS1524+08.

We consider a model to be acceptable when the difference between the observations and model ratios is less than a factor of two for every line ratio. However, considering that [OIII] is the brightest emission line, we require its ratio to match within 50%. The models that do not follow these requirements are rejected before the χ² calculation. In this manner, we ensure that all the lines fall within an acceptable range. We compute the reduced χ² for every model and choose as our best fit the model that matches our criteria and has the smallest reduced χ².

We then adopted the obtained best-fit ionization parameter, the electron density from the [SII] ratio, and the distance to the cloud as the projected distance from the nucleus to each aperture. With this, we were able to obtain a value of Q(H) for each aperture. We then integrated our SED, together with the value for Q(H), in order to estimate the bolometric luminosity required to ionize the gas to its observed state at different radii from the galaxy nucleus.

We estimated uncertainty levels for these luminosity estimations using Monte Carlo (MC) simulations, where random noise was added to the observed emission-line ratios. The simulated ratios were then fit, as described above, and the standard deviation of the best-fit values obtained are considered as the associated uncertainties.

We consider the projected distance between the nucleus and the cloud as an approximation of the true distance between the AGN and the VP. However, the true spatial distribution of the clouds is unknown. Considering the angle between the cloud and the line of sight (θ), the true distance is given by r_proj/sin θ, where r_proj is the projected in-the-sky distance. Applying this correction, the estimated light-travel time difference is given by

$$\begin{array}{c}\hfill \mathrm{\Delta }t=\frac{r_{\mathrm{proj}}}{csin\theta }(1-cos\theta ),\end{array}$$

which should be considered in our estimate of the light travel time for every cloud.

3.6. Luminosity history

In Fig. 9, we present the results of the process described above. For all the sources analyzed here, the bolometric luminosity increases with increasing radius from the nucleus. Considering the projected light travel time associated with each radius, we converted the distance from the nucleus to a derived timescale.

Fig. 9. Derived bolometric luminosity as function of projected distance from the galaxy center. Color circles correspond to the derived luminosity for each aperture along a source, using different colors for each system, as presented in the figure. The different estimates of current LBol are marked in the corresponding color for every galaxy; a small offset from the nucleus is added to improve the clarity of the different symbols. Shaded areas of corresponding colors represent the error for each luminosity value derived from MC simulations, as described in the text. There is a clear trend of increasing luminosity with increasing projected distance. This implies a significant change between ∼1–3 dex over ∼104 yr timescales.

SDSS 1510+07: We created 16 apertures along the extension of the arm-like structures formed by the ionized gas in the NS direction. The derived AGN luminosity from 2 to 14 kpc ranges from 2 × 10⁴⁴ to 10⁴⁶ erg/s. Considering that the projected distance corresponds to the distance traveled by the light from the nucleus; this implies a change of AGN luminosity of ∼1.7 dex on a timescale of 5 × 10⁴ yr.

UGC 7342: We defined 10 apertures along the ionized gas, SE to NW. The regions further from the nucleus (∼25–27 kpc) are consistent with a bolometric luminosity between 2 − 6  ×  10⁴⁷ erg/s. At ∼5 kpc, the luminosity is 5 × 10⁴⁵ erg/s. The physical distance from the nucleus reached by the ionized gas (25 kpc) in this galaxy is the largest in our sample and corresponds to ∼8 × 10⁴ yr, in which the nucleus has faded 2 dex in luminosity.

SDSS 1524+08: We generated 12 apertures, most associated with the SE ionized cloud, covering distances from 2 to 13 kpc. The bolometric luminosity at these limits corresponds to 1 × 10⁴⁴ and 3 × 10⁴⁶ erg/s, respectively. This implies that the AGN luminosity has faded 2.5 dex over 4 × 10⁴ yr.

4. Discussion

The spatial resolution provided by our VLT-MUSE observations allowed us to spatially resolve the EELRs of these sources and study their luminosity histories together with their morphologies and kinematics in detail.

4.1. Origin of the EELRs

The presence of tidal tails in the stellar component of all three objects, as well as the kinematic perturbations observed in the ionized gas, all point toward a scenario where recent mergers or galaxy interactions have been part of the sources’ histories. These interactions have likely hauled gas from the galactic disk along the tidal arms, where it has been illuminated and photoionized by the AGN. The presence of this gas can cover a larger volume, as seen from the AGN, in contrast with a galaxy where all, or most, of the gas remains within the confines of the disk. This volume filling percentage increases the likelihood of the AGN opening angle intersecting with gas and, consequently, ionizing it. This could explain the high incidence of merger or interaction features within the known VPs.

We explored the environmental properties of the sources in our sample using the group and cluster catalog of Tempel et al. (2017). This catalog was obtained using the friends-of-friends algorithm. We found that SDSS1510+07 is part of a group containing 16 sources. UGC 7342 lies within a lower density environment in a small group of three galaxies, and SDSS1524+08 is part of a group containing 107 galaxies. Both SDSS1510+07 and SDSS1524+08 can be found at lower radius and velocity in the line of sight (LOS) on a phase-space diagram, which combines information on clustocentric velocity and clustocentric radius (see Rhee et al. 2017, for a detailed description). This could indicate that these sources are in the “infalling” region of the phase-space diagram and, therefore, recently arrived in the denser regions of their respective groups.

4.2. Luminosity history

To derive estimates of the current AGN bolometric luminosity, we employ the values obtained from the luminosities measured from the [OIII] emission line, and the IR and X-ray emission. We then compare these current luminosity estimates with the bolometric luminosity derived from CLOUDY in the aperture closest to the center.

In Table 2, we show the current AGN luminosity indicators in comparison to the values previously obtained from the CLOUDY computations. For the [OIII]-derived luminosity, the flux is taken from a fit to an aperture in the central region. We assume a correction factor of L_bol/L_[OIII] ≃ 3500, with a variance of 0.38 dex, to estimate the bolometric luminosity, following Heckman et al. (2004). For the bolometric luminosities estimated from L (2–10 keV), we used an X-ray bolometric correction of κ_bol = L_bol/L_2 − 10 assuming a range of 10 < κ_bol < 100, based on typical AGN values (e.g., Vasudevan & Fabian 2009; Lusso et al. 2012). The far-infrared (FIR) luminosity from Keel et al. (2012) was used as a conservatively high estimate bolometric luminosity indicator, assuming potentially obscured AGNs; and thus, for these cases, they can be considered as upper limits.

For the sources UGC 7342, SDSS1510+07, and SDSS1524+08, we used available archival X-ray data to estimate the intrinsic AGN flux. UGC 7342 was observed with NuSTAR and XMM-Newton, but in both cases the galaxy was only marginally detected. The spectrum can be well fit by an absorbed power law, with a very high value of photon index ($4.4_{-1.3}^{+2.1}$), and a low column density ($2.0_{-1.5}^{+2.5}\times {10}^{21}$ cm⁻²). Uncertainties on the parameters are at the 90% confidence level. The emission may be consistent with star formation, primarily reflecting host-galaxy X-ray emission. SDSS1510+07 has been observed with NuSTAR and XMM-Newton, and its spectrum can be well fit by an absorbed power-law plus host-galaxy emission contributing to the soft X-rays. The derived column density is $8.9_{-2.8}^{+4.1}\times {10}^{22}$ cm⁻², while the photon index is $1.5_{-0.2}^{+0.3}$. SDSS1524+08 was observed by Chandra, and its spectrum can be well fit with a simple power law (no absorption besides Galactic) model, with $\mathrm{\Gamma }=1.8_{-0.7}^{+0.8}$. The bolometric luminosity estimated from the intrinsic 2-10 keV fluxes are listed in Table 2.

We compared the estimated bolometric luminosity with the highest values ($L_{\mathrm{Bol}}^{\max }$), reached at the largest distances from the nucleus, with the different current AGN luminosity estimators. For SDSS1510+07, $L_{\mathrm{Bol}}^{\max }$ reaches 10⁴⁷ erg s⁻¹. Comparing this value with the current estimators derived from L_FIR (from Keel et al. 2012) and L_X, and L_[OIII] suggests fading factors of 2.6, 3.9, and 4 dex, respectively. This occurs over a timescale of 4.6 × 10⁴ yr. For UGC 7342, $L_{\mathrm{Bol}}^{\max }$ corresponds to 4 × 10⁴⁶ erg s⁻¹, indicating a fading of 2.5, and 3.5 dex from the FIR, and [OIII] L_Bol estimators, respectively, over a timescale of 8 × 10⁴ yr. For this source, the X-ray luminosity estimator is likely reflecting the host galaxy X-ray emission rather than the AGN, which may be heavily obscured. For SDSS1524+08, $L_{\mathrm{Bol}}^{\max }=4\times {10}^{46}$ erg s⁻¹, comparing this luminosity with the FIR, X-ray, and [OIII] indicators, we find fading factors of 2.4, 4.7, and 2.4 dex, respectively, on timescales of 4.5 × 10⁴ yr. For NGC 5972, $L_{\mathrm{Bol}}^{\max }=2\times {10}^{46}$ erg s⁻¹, indicating a fading of 2.6, 2.5, and 2.4 dex fading over timescales of 5.5 × 10⁴ yr, as derived from the FIR, X-ray, and [OIII] L_Bol indicators. For the Teacup, $L_{\mathrm{Bol}}^{\max }$ corresponds to 9 × 10⁴⁵ erg s⁻¹ at a distance corresponding to 5.4 × 10⁴ yr. The current AGN indicators derived from L_FIR (from Keel et al. 2012), L_X and L_[OIII] indicate a difference with $L_{\mathrm{Bol}}^{\max }$ of 2.6, 1.1 and 0.7 dex. This source has the largest dispersion on the derived fading factor, ranging from 1.5 to 1000 times the current L_Bol. A significant contribution to this scatter comes from the uncertainty on the L_X, as can be seen in Fig. 9.

The derived variability on timescales of ∼10⁴ yr for our sample is consistent with previous studies of other similar sources (e.g., Lintott et al. 2009; Gagne et al. 2014; Sartori et al. 2016; Keel et al. 2017). These results could be interpreted in the context of a scenario in which the AGN duty cycle is broken down into shorter phases of 10^4 − 5 yr, resulting in a total of 10^7 − 9 yr “on” lifetime (Schawinski et al. 2015).

4.3. Structure function

We can interpret our results in the context of the framework for AGN variability presented by Sartori et al. (2018). In this work, the variability of AGNs is linked over a wide range of timescales, from days to 10⁹ yr, as derived from simulations of Novak et al. (2011). The total variability of every source is estimated as the luminosity difference between the closest and the farthest aperture from the nuclear source. This is then translated into a magnitude difference at a given time lag (τ; as defined in Eq. (1) of Sartori et al. 2018). We repeated this computation for the entire light-time travel along the EELR for every source. We obtain a value of δm (mean luminosity difference) for every bin in a time lag, which is computed from the intermediate points in the recreated luminosity curves (Fig. 9). These values are shown in Fig. 10, where they can be contrasted with observations and the structure function (SF) from hydrodynamics simulations. Furthermore, we included models BPL1, BPL2, and BPL3 from Fig. 2 in Sartori et al. (2018). These models represent the SF from light-curve simulations, derived from different power spectral densities (PSDs) of the Eddington ratio distribution. The total luminosity variability for the sources in the sample is consistent with the values obtained for the Voorwerpjes in the SF (as derived by Keel et al. 2017; Sartori et al. 2018). Our estimated δm values cover the time-lag gap between changing-look QSOs and Voorwerpjes. The observed distribution is more consistent with the SF model, BPL3, which considers a broken power law for the PSD, with α_low = −1, α_high = −2, and ν_break = 2 × 10⁻¹⁰ Hz ∼160 yr. Adding to the analysis carried out in Sartori et al. (2018), it can help to further constrain the underlying PSD and the ability to link AGN variability on different timescales.

Fig. 10. SF for magnitude difference at every measured radius (colored circles as described in the legend) and the total magnitude difference (colored triangles). We add in gray the variability data summarized in Fig. 2 of Sartori et al. (2018) for comparison. Gray circles represent the quasars’ structure function. Squares represent the δm for changing-look quasars. Downward triangles show the data for Voorwerpjes from Keel et al. (2017, filled triangles), and from X-ray data (empty triangles). The solid lines represent the BPL1, BPL2, and BPL3 models in yellow, orange, and red, respectively. The green and blue solid lines indicate SF from hydrodynamic simulations from Novak et al. (2011) and Bournaud (adapted from Sartori et al. 2018), respectively (see Sartori et al. 2018 for more details).

4.4. Physical timescales

Our analysis indicates that all the AGN in this sample have experienced a luminosity decline of 1–4 dex over the past 7 − 8 ∼ 10⁴ yr. Similar results were observed by Keel et al. (2017) for SDSS1510+07, UGC 7342, NGC 5972, and the Teacup galaxy. Their technique made use of the HST Hα surface brightness of the brightest pixel at each radius. Both procedures track each other within a constant, reconstructing similar luminosity histories. Naturally, these results raise the question of the possible physical mechanisms connected to this decline, which we explore below.

4.4.1. The role of mergers

The presence of merger or galaxy interaction signatures in all the galaxies of this sample suggest the possibility of a link between the strong AGN activity and the merger. In this scenario, the galaxy goes through a merger or interaction, which disturbs the galactic potential, hauling gas outside the disk and possibly funneling gas toward the central region; this triggers a powerful AGN phase, which then ionizes the clouds in its path, reaching tens of kiloparsecs. Following this phase, the AGN would fade over time until reaching its current luminosity. However, it is estimated that a major merger could accumulate a large fraction of galactic gas in the inner region on roughly a dynamical timescale (∼10⁸ yr). Considering the discrepancy between the merger timescale and the one derived for the AGN fading (i.e., 10^4 − 5 yr), the connection between the merger and the AGN triggering is not evident, and they could be unrelated phenomena. Therefore, the merger increases the volume-filling fraction for the AGN to illuminate distant clouds, which may translate into a selection bias.

4.4.2. Accretion-disk instabilities

Accretion-disk instabilities play a crucial role in regulating the flow of matter onto the SMBH, influencing the luminosity and activity of the AGN. The time it takes for the accretion disk to reach hydrostatic equilibrium, referred to as the dynamical timescale, can be expressed as follows:

$$\begin{array}{c}\hfill t_{\mathrm{dyn}}\sim \frac{1}{\mathrm{\Omega }},\end{array}$$

where Ω = (GM_BH/r³)^1/2 corresponds to the angular frequency for a Keplerian orbit.

The timescale for the cooling of the disk – known as the thermal timescale – is expressed as

$$\begin{array}{c}\hfill t_{\mathrm{th}}\sim {\mu }^{-1}{t}_{\mathrm{dyn}}.\end{array}$$

Viscous disturbances propagate over a distance l through the disk on timescales of l²/v_k, where v_k is the kinematic viscosity and μ is the viscosity parameter. LaMassa et al. (2015) parameterized this timescale as

$$\begin{array}{c}\hfill t_{\mathrm{vis}}\sim {\left(\frac{\mu }{0.1}\right)}^{-1}{\left(\frac{{\lambda }_{\mathrm{Edd}}}{0.03}\right)}^{-2}{\left(\frac{\eta }{0.1}\right)}^2{\left(\frac{r}{10r_g}\right)}^{7/2}\left(\frac{M}{{10}^8{M}_{\odot }}\right),\end{array}$$

where η corresponds to the radiation efficiency, r to the radius, and M to the black-hole mass.

To obtain representative timescales for these processes, we assumed μ ∼ 0.03 (Davis et al. 2010) and η ∼ 0.1 (Soltan 1982) and considered representative values for our sample of M ∼ 10⁸ M_⊙ and λ_Edd ∼ 0.17. For the gravitational radius, we considered r ∼ 100r_g for a UV-emitting region (Kato et al. 2008; LaMassa et al. 2015) and obtain the following:

$$\begin{array}{cc}\hfill t_{\mathrm{dyn}}(r/r_g=100)& \sim 16\mathrm{days},\hfill \\ \hfill t_{\mathrm{th}}(r/r_g=100)& \sim 1.5\mathrm{yr},\hfill \\ \hfill t_{\mathrm{vis}}(r/r_g=100)& \sim 2\times {10}^4\mathrm{yr},\hfill \end{array}$$

while for optical emission we considered r ∼ 200r_g, which translates into the following timescales:

$$\begin{array}{cc}\hfill t_{\mathrm{dyn}}(r/r_g=200)& \sim 45\mathrm{days},\hfill \\ \hfill t_{\mathrm{th}}(r/r_g=200)& \sim 4\mathrm{yr},\hfill \\ \hfill t_{\mathrm{vis}}(r/r_g=200)& \sim 2\times {10}^5\mathrm{yr}.\hfill \end{array}$$

This estimate indicates that both the dynamical and thermal timescales are too short to account for the suggested change in luminosity timescale from our results. Instead, the viscous timescale is more consistent with our inferred timescales. This suggests an accretion rate change in the accretion disk, following the inflow timescale of gas accretion.

4.5. Potential biases and selection effects

In Finlez et al. (2022), we discussed some caveats to take into consideration when the bolometric luminosity is derived from the best-fit CLOUDY models. These include, among others, the geometry, as we assume the real distance traveled by the light from the nuclear source corresponds to the projected distance (r_proj) from the center to the ionized cloud. A correction for this will include the inclination angle (i) of the cloud with respect to the galaxy plane as r_proj/sin(i). Additionally, it is important to consider that the density, as derived from the [SII]λλ6716, 6731 Å doublet ratio, can easily fluctuate in the 100–1000 cm⁻³ range. This is mostly explained by the steepness of the slope in the conversion curve between [SII] line ratio and derived electron density (Osterbrock & Ferland 2006). Furthermore, this tracer can saturate at high densities or in gas with a large matter-bounded fraction.

To test the impact of changing the shape of the continuum based on the He II/Hβ ratio, we tested a series of CLOUDY models focused on the median of this ratio along the EELR of every galaxy. We converted the He II/Hβ to an α in the 13.6–500 eV range following Penston & Fosbury (1978). Higher values of this ratio translate into higher α values. We ran CLOUDY models following the method described in Sect. 3.6. At higher values of α, the grid of simulations moves toward lower values of [NII]/Hα and [OIII]/Hβ in the BPT diagram. In this direction, the value of U remains mainly constant. Along the same direction, the density changes. Considering that the ionization parameter stays mainly invariable with changes in the shape of the continuum in the 13.6–500 eV range and that we used the [SII] emission-line ratio to derive the density of the gas, we chose to keep the same shape of the continuum with radius and for every object to avoid adding more degeneracies to the model fit. The suite of models that consider the dispersion of the He II/Hβ ratio along the EELRs indicates that the variations in the shape of the continuum do not translate into a large difference in the luminosity curve. Therefore, based on the He II/Hβ ratio, the corresponding changes in the shape of the continuum do not largely affect our results. More details can be found in Appendix B.

The sources analyzed were selected from the sample described in Keel et al. (2012). This sample was obtained from the Galaxy Zoo project, combining targeted and serendipitous approaches based on the visual identification of extended ionized clouds that were near or connected to a galaxy. These were later filtered in order to be consistent with AGN ionization based on BPT classification. The main bias of this selection is driven by the extension and brightness of the EELRs to be flagged as VP candidates. This would cause our method to be preferentially sensitive to the most extreme cases, as a powerful AGN in the past would illuminate up to the greatest distances from the nuclear source and strongly ionize large regions of the available gas. Furthermore, the sample selection is not based on confirmed AGNs to discover large AGN-ionized clouds with no AGN presence, which can be due to strong obscuration or acute variability during the light travel time from the nucleus to the cloud extension.

All the objects in our sample indicate fading AGNs. Considering the sample selection characteristics described above, the sources in this work may represent very extreme cases of VPs. However, in principle, both fading and rising AGNs could be identified from their EELR ionization state. Rising AGNs would show large-scale signatures of a weak nuclear source, while the small-scale (within ∼10 pc) signatures would denote more powerful activity. To further understand if the detection of fading AGNs reflects AGN evolution or is the result of a bias in our sample selection, the study of a larger, less biased sample is required. We will explore this in future work (Finlez et al., in prep.).

Binette & Robinson (1987) examined how the properties of a photoionized cloud change over time after a sharp reduction in the intensity of the ionizing photon field. They sought to determine how long the nebula could continue to glow and maintain its excitation signatures after the central AGN becomes inactive. They found that [OIII] is rapidly quenched by an efficient charge transfer combination, which proceeds as soon as the gas begins to recombine and the fraction of neutral hydrogen increases. However, hydrogen recombines relatively slowly, and its timescale can be thousands of years in EELRs. If the AGN is turned off suddenly, the [OIII] luminosity of the EELR would fade quickly, leading to a low [OIII]/Hβ ratio. However, the fading progresses through the ionized cloud on a timescale that depends on the light crossing time and, thus, on the size of the EELR. Fig. 11 shows the [OIII]/Hβ for all the sources, presenting an overall trend of a higher [OIII]/Hβ ratio along the EELR, which is in contrast with lower ratios near the nucleus. This is consistent with the interpretation of the EELRs being “ionization echoes”. Closer to the nucleus, the [OIII] emission has already declined, but this fading has not yet reached the outermost regions. Therefore, the EELR farthest from the center still retains information about the ionizing power before the shutdown.

Fig. 11. Spatially resolved [OIII]/Hβ ratio maps. A general trend of higher values of this ratio along the EELR is observed, while the ratio is lower closer to the nucleus. This suggests that near the nucleus, the [OIII] emission has already declined, but this fading has not yet reached the outer regions because of the large extension of the clouds. This is consistent with the hypothesis that EELRs are “ionization echoes”.

5. Conclusions

We present integral-field spectroscopic VLT-MUSE observations for a small sample of nearby galaxies presenting prominent EELRs at large physical extensions. We used the method outlined in Finlez et al. (2022) to analyze these sources in detail. Our main results can be summarized as follows.

1. We detect EELRs extending to 14, 27, and 14 kpc for SDSS1510+07, UGC 7342, and SDSS1524+08, respectively. The presence of tidal tails in the stellar component and the morphologies and kinematics of the ionized clouds indicate the presence of recent galaxy interactions in the sources’ histories.

2. The main ionizing mechanism, as shown through a BPT diagnostic, is consistent with AGN photoionization. In SDSS1510+07 and SDSS 1524+08, the Seyfert-like photoionization seems to follow a clear biconical distribution, with a combination of LINER and composite emission in the central region and perpendicular to the possible bicone. For UGC 7342, effectively all the ionized gas is consistent with Seyfert ionization.

3. The kinematics of the ionized gas and stellar components show complex and different characteristics for each source. SDSS1510+07 seems to show the least perturbed kinematics of the sample, with a clearly rotating stellar component. The ionized gas shows signs of rotation on a different major axis than the stellar component, while the gaseous disk appears to be warped. Ionized gas has been hauled along the stellar tidal tails, where it was likely illuminated by the nuclear source. UGC 7342 shows weaker signs of rotation and a high degree of perturbation. Signatures of a biconical outflow are observed. The gas outside the galactic disk seems to be following the tidal tails and, therefore, is located outside the galactic plane. The low velocity-dispersion values along these areas of the EELR, further from the nucleus, indicate that the gas may not be outflowing and is kinematically turbulent due to the past galaxy interaction. For SDSS1524+08, there are weak signs of rotation only from a nuclear disk, both in the stellar and ionized gas components. However, on larger scales, the galaxy lacks rotational support. The velocities in the LOS and the velocity dispersion indicate a possible outflow along the bicone formed by the NLR. The main feature of the EELR, a large ionized cloud at the SE, shows low velocity dispersion and is located at the same spatial location as the tidal arms. Given that this cloud is in the direction of the bicone, this gas was possibly hauled by the tidal arms or even pushed by the biconical outflow and later illuminated by the AGN.

4. We used the CLOUDY photoionization code to generate a grid of models covering a range of ionization parameters to derive the luminosity required to ionize the gas at every aperture to the observed state. For all the sources in our sample, we see an increase in luminosity with increasing radius. This suggests that the AGN has faded from ∼0.2 − 3 dex on timescales of 4 − 8 × 10⁴ yr. These variability amplitudes and related timescales are in good agreement with previous luminosity history studies and AGN variability models (e.g., Lintott et al. 2009; Keel et al. 2012; Gagne et al. 2014; Sartori et al. 2018).

5. The dramatic changes in luminosities described above occur on timescales that are in agreement with a scenario in which the AGN duty cycle is broken down into shorter phases of 10^4 − 5 yr, resulting in a total 10^7 − 9 yr “on” lifetime (Schawinski et al. 2015). This “flickering” of the nuclear source could be responsible for the observed timescale variability. Considering that there remains significant ionized gas (and likely neutral and molecular as well) in the vicinity, it is reasonable to expect that a portion will eventually reach the nucleus within time frames of < 10^7 − 9 yr.

Future observations of a larger and potentially less biased sample are necessary to further consolidate the results reported here regarding large-timescale AGN variability. Considering the biases and selection effects discussed in 4.5, different regions of the AGN parameter space and their host galaxies must be explored to understand the reach of our conclusions. We aim to explore this in a subsequent study, using a sample selected from the BAT AGN Spectroscopic Survey (BASS; Koss et al. 2017), which is the first large all-sky survey of hard X-ray-selected AGNs with over 1,000 sources, providing the most detailed and complete census of SMBHs in the local Universe. One of our sources, the Teacup, is part of this sample. This large sample of bright AGNs, selected uniformly and robustly, provides an unbiased sample of currently luminous AGNs in the local Universe with well-measured SMBH masses and Eddington ratios.

Data availability

The derived line maps are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/703/A63

Acknowledgments

We would like to thank the anonymous referee for their comments that improved the quality of this article. We acknowledge support from FONDECYT through Postdoctoral grant 3220751 (CF), Regular 1190818 (ET, FEB), 1200495 (ET, FEB), 1250821 (ET, FEB), and 3230310 (YD); ANID grants CATA-Basal AFB-170002 (ET, FEB, CF), FB210003 (ET, FEB, CF) and ACE210002 (CF, ET, and FEB); Millennium Nucleus NCN19_058 (TITANs; ET, CF, NN); Millennium Science Initiative Program – ICN12_009 (FEB); and European Union’s HE ERC Starting Grant No. 101040227 - WINGS (GV).

References

Appendix A: Stellar and gas component maps for all targets

In this section we present the figures for the objects UGC 7342 and SDSS1524+08, see Sect. 3.1 and 3.2 for corresponding text.

Fig. A.1. Flux (left), velocity (middle), and velocity dispersion (right) maps for the stellar component of UGC 7342. The integrated flux maps show the extension of the stellar component, which is highlighted with the ellipses fitted to the isophotal profiles. Fitted ellipses are shown as white contours in all the maps. Beyond the centrally concentrated high surface brightness stellar structure, tidal arms dominate the distribution of the stellar component. The velocity map in the middle panel shows a gradient that spatially coincides with the centrally concentrated component, indicating the possible presence of rotation. Beyond the centrally concentrated distribution, the kinematics of the stellar component are clearly perturbed due to the streaming motions. The velocity dispersion map shows a higher velocity concentration in the nuclear region, indicating some possible rotational support.

Fig. A.2. Flux (left panel), velocity (middle panel), and velocity dispersion (right panel) maps for the stellar component of SDSS1524+08. Overplotted in gray are the ellipses fitted to the isophotal profiles. This fit indicates a close-to-face-on disk. Outside the main centrally-concentrated stellar distribution, there are clear tidal tails that dominate the stellar distribution. The velocity map indicates little to no sign of rotation, with an apparent small gradient in velocity in the central region. The velocity dispersion map shows a high velocity dispersion in the redshifted region NE of the nucleus and seemingly along the tidal tail that goes E to N. This indicates very low rotational support in the stellar distribution.

Fig. A.3. Normalized flux (left), velocity (middle), and velocity dispersion (right) maps for UGC 7342 based on the [OIII] emission line. The black contours represent the values (in km/s) of the parameter W80, above the velocity limit used to fit a secondary Gaussian component to the emission line fit.

Fig. A.4. Normalized flux (left), velocity (middle), and velocity dispersion (right) maps for SDSS1524+08 obtained from measurements of the [OIII] emission line. The black contours represent the values (in km/s) of the parameter W80, above the velocity limit used to fit a secondary Gaussian component to the emission line fit.

Fig. A.5. Velocity map for the stellar component of UGC 7342 (left). Best fit kinematic Barolo model (middle). Residual map, that is, velocity map subtracted from the model map (right)

Fig. A.6. Velocity map for the ionized gas component of UGC 7342 (left). Best fit kinematic 3D-Barolo model (middle). Residual map, that is, velocity map subtracted from the model map (right)

Fig. A.7. Velocity map for the stellar component of SDSS1524+08 (left). Best fit kinematic 3D-Barolo model (middle). Residual map (right). The residual map indicates possible rotation along PA 22° in the inner 10″.

Fig. A.8. Velocity map for the ionized gas component of SDSS1524+08 (left). Best fit kinematic 3D-Barolo model (middle). Residual map (right). The residual map shows that rotation along PA 22° can explain some of the features along this PA, indicating there could be a rotating component in the same PA of the stellar component. However, the strong components along PA 320° clearly dominate the fit.

Fig. A.9. Rotation curve for SDSS1524+08, extracted along the major axis for the stellar component (red) and the [OIII] emission velocity map (black). Following the same color scheme, we show the fitted rotation curves for the stars and the gas.

Appendix B: CLOUDY models parameters

To assess the impact of the column density (N_H) on our CLOUDY models, we generated a suite of models with N_H values spanning 10¹⁹ − 10²³ cm⁻². This range was chosen to cover the expected values for both radiation-bounded and matter-bounded clouds, allowing us to evaluate the influence of this parameter on our results.

To illustrate the effect of N_H on the derived bolometric luminosity (L_bol), Figure B.1 presents a modified version of Figure 9, displaying five lines per object corresponding to log(N_H) values of 19, 20, 21, 22, and 23.

We find that while variations in N_H influence the ionization parameter (U) and, consequently, the estimated L_bol, the overall trend in L_bol remains consistent across all N_H values. Specifically, the estimated L_bol increases with radius, and the difference in L_bol between the nearest and furthest radii remains largely unchanged. However, differences in N_H do affect the fading factor when comparing the earliest L_bol estimate with the current AGN L_bol, derived from observations at different wavelengths. As a result, we consider the difference between the L_bol estimated for N_H = 10¹⁹ cm⁻² and our earliest L_bol estimate to represent a lower limit on the discrepancy between the current and earlier L_bol values.

To estimate the impact of the shape of the continuum chosen, we have run a suite of models with different values of α in the range 13.6–500 eV, based on the values of the He II/Hβ ratio observed along the VPs. In Fig. B.2 we show the values of the He II/Hβ along the VP of every source, at different projected distances from the galaxy nucleus. In Table B.1 we tabulate the median values and the 1σ dispersion of these values. We calculate the α values corresponding to the median and 1σ values of the He II/Hβ ratio, following Penston & Fosbury (1978). In Fig. B.3 we plot the luminosity curves corresponding to the median and ±1σ values of the He II/Hβ ratio, for every galaxy. We find that the variations observed in the He II/Hβ ratio translate into mostly minor variations to the luminosity curve. The exception is SDSS 1510+07, for which the variations are larger. However, the trend of increasing luminosity with radius remains.

Fig. B.1. Derived bolometric luminosity as a function of the projected distance from the galaxy nucleus, as described in Fig. 9. The model outlined in Sect. 3.5 was run with five different values of N(H), ranging from 1019cm−2 (darker color shades) to 1023cm−2 (lighter color shades). The trend of increasing luminosity with projected distance persists across all values of N(H).

Fig. B.2. He II/Hβ flux ratio over projected distances from the center. Colors are the same as Fig. 9.

Fig. B.3. Derived bolometric luminosity as a function of the projected distance from the galaxy nucleus, as described in Fig. 9. The model outlined in Sect. 3.5 was run with different values of alpha in the range 13.6–500 eV. This change in the shape of the continuum correspond to the dispersion of values in the He II/Hβ ratio, showing as a solid line the median value, and in dotted lines (shaded in between) the ±1σ difference. For most of the galaxies in the sample, the change in the shape of the continuum translates into a small difference in the luminosity curve. SDSS 1510, which despite showing a larger spread in Lbol, maintains the trend of increasing luminosity with projected distance.

All Tables

All Figures

Fig. 1. Composite color images for targets analyzed in this work with red, white, and green pseudo-colors. Images were created using three bands extracted from the data cube. Green was created with a spectral window that encloses the [OIII]5007 emission lines to represent the highly ionized gas; red was created with a spectral window that encloses Hα 6563 Å, while white represents the stellar continuum, with a spectral window of 7000–9000 Å. In all images, north is up and east is left.

In the text

Fig. 2. Flux, velocity, and velocity dispersion maps for the stellar component of SDSS1510+07. The left panel shows the integrated flux, obtained from collapsing the data cube around 8000 Å. Overplotted in white are the isophote profile fits, which show the extension and inclination of a stellar component; this is consistent with a disk. The velocity map (middle panel) shows that the kinematics of the stellar component are rotation-dominated. The velocity-dispersion map (right panel) shows higher velocities concentrated in the center, as expected from a rotation-dominated disk.

In the text

	Fig. 3. Normalized flux (left), velocity (middle), and velocity-dispersion (right) maps for SDSS1510+07 based on the measurements of the [OIII] emission line. The black contours represent the values (in km/s) of the parameter W80, above the velocity limit used to fit a secondary Gaussian component to the emission-line fit.
In the text

Fig. 4. Spatially resolved BPT classification maps based on [OIII]/Hβ versus [NII]/Hα BPT diagrams for SDSS1510+07, UGC 7342, SDSS1524+08 (top from left to right), NGC 5972, and the Teacup (bottom from left to right). The color bars indicate the different BPT regions, where green indicates the spaxels fall in the Seyfert section of the diagram, red for LINER, blue for composite, and yellow for star forming. Every color is weighted with the [OIII] emission-line flux to more clearly show the bright photoionized regions.

In the text

	Fig. 5. Velocity map for stellar component of SDSS1510+07 (left). Best-fit kinematic Barolo model (middle). Residual map; i.e., velocity map subtracted from the model map (right).
In the text

	Fig. 6. Velocity maps for emission lines [OIII] (top row) and Hα (bottom row). Velocity map for the ionized gas component of SDSS1510+07 (left). Best-fit kinematic 3D-Barolo model (middle). Residual map; i.e., velocity map subtracted from the model map (right).
In the text

	Fig. 7. Spatially resolved He II/Hβ ratio maps. The ratio remains largely constant along the EELRs for all sources, with lower values in the central region. This area is excluded from our luminosity history analysis. Since this ratio reflects the continuum shape and is mostly constant along the ionized clouds, we assume a fixed spectral shape of the continuum with radius.
In the text

	Fig. 8. Apertures used for luminosity history analysis marked as black circles over [SII] emission line normalized flux maps. Left to right: SDSS1510+07, UGC 7342, and SDSS1524+08.
In the text

Fig. 9. Derived bolometric luminosity as function of projected distance from the galaxy center. Color circles correspond to the derived luminosity for each aperture along a source, using different colors for each system, as presented in the figure. The different estimates of current LBol are marked in the corresponding color for every galaxy; a small offset from the nucleus is added to improve the clarity of the different symbols. Shaded areas of corresponding colors represent the error for each luminosity value derived from MC simulations, as described in the text. There is a clear trend of increasing luminosity with increasing projected distance. This implies a significant change between ∼1–3 dex over ∼104 yr timescales.

In the text

Fig. 10. SF for magnitude difference at every measured radius (colored circles as described in the legend) and the total magnitude difference (colored triangles). We add in gray the variability data summarized in Fig. 2 of Sartori et al. (2018) for comparison. Gray circles represent the quasars’ structure function. Squares represent the δm for changing-look quasars. Downward triangles show the data for Voorwerpjes from Keel et al. (2017, filled triangles), and from X-ray data (empty triangles). The solid lines represent the BPL1, BPL2, and BPL3 models in yellow, orange, and red, respectively. The green and blue solid lines indicate SF from hydrodynamic simulations from Novak et al. (2011) and Bournaud (adapted from Sartori et al. 2018), respectively (see Sartori et al. 2018 for more details).

In the text

Fig. 11. Spatially resolved [OIII]/Hβ ratio maps. A general trend of higher values of this ratio along the EELR is observed, while the ratio is lower closer to the nucleus. This suggests that near the nucleus, the [OIII] emission has already declined, but this fading has not yet reached the outer regions because of the large extension of the clouds. This is consistent with the hypothesis that EELRs are “ionization echoes”.

In the text

Fig. A.1. Flux (left), velocity (middle), and velocity dispersion (right) maps for the stellar component of UGC 7342. The integrated flux maps show the extension of the stellar component, which is highlighted with the ellipses fitted to the isophotal profiles. Fitted ellipses are shown as white contours in all the maps. Beyond the centrally concentrated high surface brightness stellar structure, tidal arms dominate the distribution of the stellar component. The velocity map in the middle panel shows a gradient that spatially coincides with the centrally concentrated component, indicating the possible presence of rotation. Beyond the centrally concentrated distribution, the kinematics of the stellar component are clearly perturbed due to the streaming motions. The velocity dispersion map shows a higher velocity concentration in the nuclear region, indicating some possible rotational support.

In the text

Fig. A.2. Flux (left panel), velocity (middle panel), and velocity dispersion (right panel) maps for the stellar component of SDSS1524+08. Overplotted in gray are the ellipses fitted to the isophotal profiles. This fit indicates a close-to-face-on disk. Outside the main centrally-concentrated stellar distribution, there are clear tidal tails that dominate the stellar distribution. The velocity map indicates little to no sign of rotation, with an apparent small gradient in velocity in the central region. The velocity dispersion map shows a high velocity dispersion in the redshifted region NE of the nucleus and seemingly along the tidal tail that goes E to N. This indicates very low rotational support in the stellar distribution.

In the text

	Fig. A.3. Normalized flux (left), velocity (middle), and velocity dispersion (right) maps for UGC 7342 based on the [OIII] emission line. The black contours represent the values (in km/s) of the parameter W80, above the velocity limit used to fit a secondary Gaussian component to the emission line fit.
In the text

	Fig. A.4. Normalized flux (left), velocity (middle), and velocity dispersion (right) maps for SDSS1524+08 obtained from measurements of the [OIII] emission line. The black contours represent the values (in km/s) of the parameter W80, above the velocity limit used to fit a secondary Gaussian component to the emission line fit.
In the text

	Fig. A.5. Velocity map for the stellar component of UGC 7342 (left). Best fit kinematic Barolo model (middle). Residual map, that is, velocity map subtracted from the model map (right)
In the text

	Fig. A.6. Velocity map for the ionized gas component of UGC 7342 (left). Best fit kinematic 3D-Barolo model (middle). Residual map, that is, velocity map subtracted from the model map (right)
In the text

	Fig. A.7. Velocity map for the stellar component of SDSS1524+08 (left). Best fit kinematic 3D-Barolo model (middle). Residual map (right). The residual map indicates possible rotation along PA 22° in the inner 10″.
In the text

	Fig. A.8. Velocity map for the ionized gas component of SDSS1524+08 (left). Best fit kinematic 3D-Barolo model (middle). Residual map (right). The residual map shows that rotation along PA 22° can explain some of the features along this PA, indicating there could be a rotating component in the same PA of the stellar component. However, the strong components along PA 320° clearly dominate the fit.
In the text

	Fig. A.9. Rotation curve for SDSS1524+08, extracted along the major axis for the stellar component (red) and the [OIII] emission velocity map (black). Following the same color scheme, we show the fitted rotation curves for the stars and the gas.
In the text

	Fig. B.1. Derived bolometric luminosity as a function of the projected distance from the galaxy nucleus, as described in Fig. 9. The model outlined in Sect. 3.5 was run with five different values of N(H), ranging from 1019cm−2 (darker color shades) to 1023cm−2 (lighter color shades). The trend of increasing luminosity with projected distance persists across all values of N(H).
In the text

	Fig. B.2. He II/Hβ flux ratio over projected distances from the center. Colors are the same as Fig. 9.
In the text

Fig. B.3. Derived bolometric luminosity as a function of the projected distance from the galaxy nucleus, as described in Fig. 9. The model outlined in Sect. 3.5 was run with different values of alpha in the range 13.6–500 eV. This change in the shape of the continuum correspond to the dispersion of values in the He II/Hβ ratio, showing as a solid line the median value, and in dotted lines (shaded in between) the ±1σ difference. For most of the galaxies in the sample, the change in the shape of the continuum translates into a small difference in the luminosity curve. SDSS 1510, which despite showing a larger spread in Lbol, maintains the trend of increasing luminosity with projected distance.

In the text