The Close AGN Reference Survey (CARS): IFU survey data and the BH mass dependence of long-term AGN variability

[Abridged] AGN are thought to be intimately connected with their host galaxies through feeding and feedback processes. A spatially resolved multiwavelength survey is required to map the interaction of AGN with their host galaxies on different spatial scales and different phases of the ISM. The goal of CARS is to obtain the necessary spatially resolved multiwavelength observations for an unbiased sample of local unobscured luminous AGN. We present the overall CARS survey design and the associated wide-field optical IFU spectroscopy for all 41 CARS targets at z<0.06 randomly selected from the Hamburg/ESO survey of luminous unobscured AGN. This data set provides the backbone of CARS and allows us to characterize host galaxy morphologies, AGN parameters, precise systemic redshifts, and ionized gas distributions including excitation conditions, kinematics, and metallicities in unprecedented detail. We focus our study on the size of the ENLR which has been traditionally connected to AGN luminosity. Given the large scatter in the ENLR size-luminosity relation, we performed a large parameter search to identify potentially more fundamental relations. Remarkably, we identified the strongest correlation between the maximum projected ENLR size and the black hole mass, consistent with an $R_\mathrm{ENLR,max}\sim M_\mathrm{BH}^{0.5}$ relationship. We interpret the maximum ENLR size as a timescale indicator of a single BH radiative-efficient accretion episode for which we inferred log(t_AGN) = (0.45+- 0.08)log(M_BH)+1.78 using forward modeling. The extrapolation of our inferred relation toward higher BH masses is consistent with an independent lifetime estimate from the HeII proximity zones around luminous AGN at z~3. While our proposed link between the BH mass and AGN lifetime might be a secondary correlation itself or impacted by unknown biases, it has a few relevant implications if confirmed.


Introduction
Active galactic nuclei (AGN) have drawn a lot of attention over the last decades because they have been beacons for the existence and demographics of super-massive black holes (BHs) throughout the history of the Universe (e.g., Soltan 1982;Kollmeier et al. 2006;Greene & Ho 2007;Schulze & Wisotzki 2010;Kelly & Shen 2013;Kormendy & Ho 2013;Schulze et al. 2015;Bañados et al. 2018). The release of gravitational binding energy via accretion of matter onto these BHs is expected to have a profound impact on the evolution of their host galaxies (e.g., Silk & Rees 1998;Granato et al. 2004;Di Matteo et al. 2005; Hopkins et al. 2008;Somerville et al. 2008;Fabian 2012;Harrison 2017;Harrison et al. 2018;Gaspari et al. 2019;Nelson et al. 2019). Large spectroscopic surveys such as the Sloan Digital Sky Survey (SDSS, York et al. 2000;Abazajian et al. 2009), the 2df QSO redshift survey (2QZ, Croom et al. 2004), the VIMOS VLT Deep Survey (VVDS, Le Fèvre et al. 2013) or the VIMOS Public Extragalactic Redshift Survey (VIPERS, Scodeggio et al. 2018) in combination with several deep X-ray surveys taken with ROSAT (Voges et al. 1999), Chandra (Elvis et al. 2009;Xue et al. 2011), XMM-Newton (Pierre et al. 2016), and eROSITA (Predehl et al. 2021) as well as large radio surveys (e.g., Becker et al. 1995;Condon et al. 1998;Smolčić et al. 2017;Shimwell et al. 2019;Lacy et al. 2020;Gordon et al. 2021) have provided an enormous data set to characterize the AGN population and its evolution with redshift in great detail. While the standard model for the AGN central engine has been successful in unifying the various classes of AGN appearance (Antonucci 1993;Urry & Padovani 1995;Padovani et al. 2017), some aspects such as changing-look AGN (CLAGN, e.g. MacLeod et al. 2016;Ruan et al. 2016;Graham et al. 2017;Noda & Done 2018) and tidal-disruption events (e.g., Gezari et al. 2009;Merloni et al. 2015;Auchettl et al. 2017) are just being explored more extensively in the time domain.
It remains challenging to understand the impact of AGN on their host galaxies as a function of AGN luminosity and dominant mode of accretion, that is radiatively efficient or inefficient, across cosmic time. Modern cosmological galaxy evolution simulations implement energy coupling between the AGN and its host galaxy to prevent excessive star formation in massive halos (e.g., Vogelsberger et al. 2014;Crain et al. 2015;Steinborn et al. 2015;Anglés-Alcázar et al. 2017;Weinberger et al. 2017). However, the different physical subgrid prescriptions of AGN feedback in those simulations yield similar results for the galaxy population despite different assumptions. This suggests that the necessity of the impact is well established, but not the accurate physical process at play. Observations of AGN and their host galaxies are therefore required to reveal the various physical processes operating on different spatial and temporal scales.
All these processes are known to exist, but how they play together, what relative impact they have on galaxy evolution, and how frequent various process occur, remains elusive.
A major complication in understanding the impact of AGN on galaxies lies in the complexity of the galaxy itself. Galaxies exhibit a multiphase ISM with various different temperatures that cannot be easily observed together (e.g., Cicone et al. 2018); their evolution is not only driven by galaxy mass but also linked with their specific environment (e.g., Peng et al. 2010b); and the dynamical time of galaxies (several hundred million years) is probably much longer than the short outbursts of bright AGN phases (e.g., Hickox et al. 2014). Hence, a comprehensive study of AGN host galaxies probing all the different aspects of their evolution across the underlying galaxy population is needed to quantify the impact of the various coupling mechanisms.
In this paper, we present the first set of optical IFU observations taken as part of CARS, which represents the backbone of our first public data release. We describe the reduction and processing of the data including the adopted AGN-host galaxy deblending method, the combined stellar continuum and emission line fitting procedure, and a morphological characterization of the host galaxies. This paper is one in a series of three core papers from the CARS survey and focuses on the ENLR properties and its scaling relations. As a key result we interpret the size of the ENLR as a long-term AGN variability indicator and report a BH mass dependence of the AGN variability on 10 4 -10 6 yr timescale.
The paper is organized as follows. In Sect. 2 we present the survey aims and basic sample characteristics for CARS. The details of the IFU observations and data reduction are described in Sect. 3 which is followed by the data analysis in Sect. 4 including AGN-host galaxy deblending, host galaxy characterization, emission-line mapping and AGN parameter determination. We present and discuss the results in Sect. 5 and close the paper with our conclusions in Sect. 6. Ancillary information on our AGN variability model and figures for the entire sample can be found in the Appendix. In this paper, we assume a concordance cosmology with H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7.

Survey strategy and main aims
Galaxies are composed of stars and a complex multiphase interstellar medium that is also linked to the larger gas reservoir of the intergalactic medium. It has been argued that investigating only a single gas-phase provides a very limited or biased picture about AGN outflows, feedback and feeding (e.g., Cicone et al. 2018). In addition, the gas-phases often need to be spatially resolved to properly differentiate between possible AGN feeding and feedback processes (e.g., Husemann & Harrison 2018). While large statistical studies of AGN with integrated properties can address some of the questions from a statistical point-of-view, it is impossible to discern all internal galaxy processes based solely on integrated measurements.
Various surveys have started to obtain detailed spatially resolved multiphase observations of AGN host galaxies to complement large statistical investigations. For example, the KMOS AGN Survey at High-z (KASHz, Harrison et al. 2016;Scholtz et al. 2020), Survey for Unveiling the Physics and the Effect of Radiative feedback (SUPER, Circosta et al. 2018;Kakkad et al. 2020;Circosta et al. 2021) and the WISE/SDSS-selected hyperluminous quasar survey (WISSH, Bischetti et al. 2017;Vietri et al. 2018;Bischetti et al. 2021) map the properties of outflows and host galaxies of luminous high-redshift AGN by combining NIR spectroscopy with molecular gas and dust mapping with ALMA. At low redshifts, the BAT AGN Spectroscopic Survey (BASS, Koss et al. 2017) provides a multiwavelength spectroscopic characterization for a large AGN sample, while the Nuclei of Galaxies survey (NUGA, García-Burillo et al. 2005;Haan et al. 2009), the Galactic Activity, Torus and Outflow Survey (GATOS, Alonso-Herrero et al. 2019), the Measuring Active Galactic Nuclei Under MUSE Microscope survey (MAG-NUM, Cresci et al. 2015;Mingozzi et al. 2019), the Local Luminous AGN with Matched Analogs survey (LLAMA, Rosario et al. 2018;Caglar et al. 2020), and the Quasar Feedback Survey (Jarvis et al. 2021) provide extensive spatially resolved data for smaller samples. Such low-redshift samples achieve much better spatial resolution and sensitivity than possible at high redshifts, but are often limited to low-luminosity AGN or selecting a specific AGN parameter space.
The aim of CARS is to provide a representative reference data set for a significant sample of the most luminous AGN in the nearby Universe that bridges the gap in AGN luminosity between the low-redshift and high-redshift surveys. This requires a wide area AGN survey as a parent sample which contains such rare AGN in the local Universe with sufficient numbers. The CARS AGN survey is limited to z < 0.06 to ensure a subkiloparsec spatial resolution across the host galaxies in most wavelength regimes. While the backbone of the survey is optical IFU spectroscopy of the ionized gas phase, we are combining the IFU data with ALMA, VLA and SOFIA observations to map the cold gas, Chandra for the hot gas, and deep panchromatic imaging from various resources. With such a comprehensive data set, CARS aims to tackle a few major open questions about the connection between AGN and their host galaxies in terms of feeding and feedback cycle such as the incidence and properties of multiphase gas outflows, relative role of radiation and radio jet-driven outflows, suppression or enhancement of star formation, timescale of AGN accretion and outflows, localized versus. global impact of AGN feedback, and signature for fueling mechanisms on host galaxy scales.
The value of this multiwavelength approach for addressing those questions has already been explored for the CARS galaxy HE 1353−1917(Husemann et al. 2019bSmirnova-Pinchukova et al. 2019). Furthermore, limits on detecting the hot gas in Xrays with current facilities has been highlighted in (Powell et al. 2018) for two objects (Mrk 1044and HE 1353-1917. In addition to our discovery of a prominent changing look (CLAGN) event in Mrk 1018 (McElroy et al. 2016;Husemann et al. 2016;Krumpe et al. 2017) within CARS, we performed a systematic study of star formation within the bars of the galaxies presented by Neumann et al. (2019), which indicates differences in star formation rates that are likely related to secular evolution rather than AGN feedback.

The unobscured AGN parent sample for CARS
The CARS sample is entirely focused on unobscured (type 1) AGN for which black hole masses and accretion rates can be easily estimated from the accretion disk brightness and broad line characteristics from a single spectrum (e.g., Peterson & Wandel 2000;Kaspi et al. 2000). In order to select the rare bright AGN in the nearby Universe we use the Hamburg/ESO survey of bright UV-excess sources (HES Wisotzki et al. 2000) which covers a large area of ∼9500 deg 2 in the Southern sky. The limiting magnitude of HES is B J < 17.3 mag with a dispersion of 0.5 mag across fields and allows a selection of bright type 1 AGN up to a redshift of z ≈ 3.2. Follow-up spectroscopy as part of HES is used to confirm AGN and determine accurate redshifts (Wisotzki et al. 2000;Schulze et al. 2009).
Applying a redshift cut of z ≤ 0.06 to the HES catalog leads to a sample of 99 AGN, which defines the CARS parent sample. The chosen redshift cut ensures a subkiloparsec resolution in seeing-limited optical and near-infrared observations. This is necessary for a detailed structural decomposition and in particular enables us to separate the AGN point-source emission from the host galaxy components (e.g., Busch et al. 2014). The exact value was chosen such that the CO(2-0) band-heads are still observable in NIR K-band in order to map the stellar kinematics and populations from high-angular NIR observations (e.g., Fischer et al. 2006;Busch et al. 2016).
Bolometric luminosities, derived from the X-ray luminosity, are lower than those of high-z QSOs by about one order of magnitude, but the Eddington ratios of CARS targets are similar to more high-z AGN (Laha et al. 2018). Molecular gas masses ranging from 0.4 × 10 9 M to 9.7 × 10 9 M (Bertram et al. 2007) and neutral atomic gas masses ranging from 1.1 × 10 9 M to 3.8 × 10 10 M (König et al. 2009) also fall between those of nearby low-luminosity AGN and high-z luminous QSOs (e.g., Kakkad et al. 2017;Circosta et al. 2018). Our sources can therefore serve as an important bridge between these AGN populations (Moser et al. 2012).

The CARS AGN sample and data release 1
From the parent sample of 99 objects described above, the CARS sample is a representative subset of 41 galaxies that were randomly chosen from the parent sample for follow-up single-dish CO(1-0) observations as presented by Bertram et al. (2007). A gallery of broad-band images for the entire sample is shown in Fig. 2. The availability of single-dish submillimeter observations for this sample was considered as an advantage for future followup spatially resolved observations with ALMA and the VLA. All CARS targets are listed in Table 2.3 together with the absolute B J -band magnitude from HES, the continuum radio flux from the FIRST or NVSS survey, the soft X-ray flux from ROSAT as Nuclear B J -band absolute magnitude distribution of the CARS sample in comparison to the local (z ∼ 0) type 1 AGN luminosity function. The data points represent measurements of the combined SDSS+HES luminosity function as inferred by Schulze et al. (2009) and the solid line is the corresponding best-fit double Schechter function as reported by Schulze et al. (2009). The vertical dashed line highlights the break luminosity M B, * of the Schechter function. The red histogram shows the nuclear B J -band absolute magnitude distribution of the CARS sample. well as the integrated H 2 (Bertram et al. 2007) and H i (König et al. 2009) line fluxes.
The overall AGN luminosity function at low redshifts was determined by Schulze et al. (2009). They combined the large SDSS and HES type 1 AGN samples and derived a double Schechter function to describe the observed luminosity function as shown in Fig. 1. The distribution in AGN luminosity for the CARS sample peaks around an absolute magnitude of M B = −20.5 mag which is one magnitude brighter than the break luminosity M B, * = −19.46 mag of the best-fit Schechter function (Schulze et al. 2009). Hence, the CARS sample is confirmed to probe the bright tail of the AGN luminosity function in the local Universe, but still lacks the bright QSOs due the limited volume imposed by the redshift cut. An exception in this regard is the faint AGN HE 1248−1356 with an absolute magnitude of M B J = −16.99 mag. It has only been detected by HES due to its low redshift of z = 0.0145.
As a first step in the CARS project we investigated optical IFU spectroscopy as well as panchromatic images collected from various resources. We release all those observations and highlevel data products to the community as part of the CARS data release 1 (DR1) that can be accessed at http://cars.aip.de. The details of the processing of the distributed data products are described in this paper for the IFU data and ionized gas measurements, in Singha et al. (2021) for the circumnuclear [O iii] kinematics and in Smirnova-Pinchukova et al. (2021) for the overall spectral energy distributions and star formation properties. One important result from the CARS IFU follow-up observations related to the sample is that two of the CARS targets, HE 0045−2145 and HE 0150−0344, turned out to be starburst galaxies which were initially identified as AGN by HES. The exceptionally blue continuum and broad line components from starburst-driven outflows could not be distinguished from type 1 AGN signatures in the low-resolution HES spectra. We can exclude changing-look AGN as a cause for this discrepancy after comparing with the original HES spectra (Schulze priv. communications). This misclassification in the parent sample is reducing the actual number of AGN targets from 41 to 39, but we keep the other two galaxies in the DR for completeness and determine all relevant parameters.

Multi-unit spectroscopic explorer
Follow-up optical IFU spectroscopy of the CARS sample was mainly taken with the Multi-Unit Spectroscopic Explorer (MUSE, Bacon et al. 2010Bacon et al. , 2014 at the Very Large Telescope (VLT) as part of two filler programs 094.B-0345(A) and 095.B-0015(A) (PI: B. Husemann). MUSE maps a large 1 × 1 fieldof-view (FoV) at a sampling of 0 . 2 with a spectral coverage of 4750-9300 Å and a spectral resolution of 1800 < R < 3500. 35 targets from the CARS sample have been obtained with our program. The total integration times range from 600 s up to 2400 s split into several exposures as listed in Table 2. Some observations were repeated due to nonoptimal observing conditions, but some of those exposures are added if data appeared usable after inspection. Because all those CARS targets are significantly smaller than the MUSE FoV, no dedicated sky exposures were taken and sky background is estimated from the blank areas within the science exposures.
In addition to our two dedicated programs, 2 CARS targets were observed by other programs and are publicly available in the ESO archives. HE 0212−0059 (Mrk 590) is a well-known changing-look AGN and was targeted as part of program 099.B-0249(A) (PI: Raimundo) for 8 × 1100 s with interleaved dedicated sky frame exposures due to the large size of the galaxy (Raimundo et al. 2019). HE 1237−0504 was observed as part of the MUSE Atlas of Disc (MAD) project (Erroz-Ferrer et al. 2019) under program 099.B-0242(B) (PI: Carollo) with 4 × 900 s exposures. Unfortunately, no dedicated sky frames were obtained for this galaxy although it is significantly more extended than the MUSE FoV.
We process all data with the official MUSE pipeline v2.8.1 (Weilbacher et al. 2012(Weilbacher et al. , 2014(Weilbacher et al. , 2020 following standard procedures for reducing the bias frames, the continuum lamp exposures to trace the slit, the arc lamp exposure to establish the wavelength solution, the standard star exposures for flux calibration, and the twilight flat frames. Those calibrations are applied to each individual exposure. A simple mean sky spectrum is created from the area within the FoV selected from the lower 10% of the flux distribution in the white light images. This mean sky spectrum is then subtracted from the full cube which leaves significant sky line residuals which we further suppress as described in more detail in the next section. Telluric absorption bands are approximately corrected by dividing with the normalized transmission from the standard star exposures closest in time. Given that standard stars are usually not taken close in time, some significant residuals can still be present. For some targets a bright star is in the field which we used to calibrate and subtract the residual telluric features whenever possible.

Suppression of sky line residuals
Due to significant variation of the line spread function (LSF) in MUSE it is very challenging to accurately subtract the sky line across the entire FoV. A principal component analysis (PCA) scheme has shown to be very efficient in those cases to characterize the spectral residuals (e.g., Kurtz & Mink 2000;Sharp & Birchall 2010;Bai et al. 2017). A minimal set of orthogonal basis spectra are created that can describe any sky line residual as a linear super-position. Fitting those linear contributions to individual spectra allows us to remove most of the systematic residual sky line emission patterns. While a general and flexible software tool has been designed for MUSE, called Zurich Atmosphere Purge (ZAP, Soto et al. 2016), the parameters still need to be optimized depending on specific content of the observations. For CARS we created a simplified sky residual suppression algorithm based on the PCA approach which requires less free parameters to be set and is therefore more robust with respect to the actual science content of the observations. The code termed CubePCA 1 is publicly available and we briefly describe the algorithm in the following. First, a 2D image with the same spatial dimension of the MUSE cube is created to flag empty sky regions (1 for sky and 0 otherwise). The corresponding Eigenspectra are computed from 1 https://git.io/cubepca this selected region of empty sky spectra which were limited to a maximum subset of 20,000 spectra. Any remaining continuum signal is subtracted beforehand by running a median filter of 25 pixel width. The derived orthogonal Eigenspectra are then used to clean the sky line residuals across the entire science cube. The science cube is temporarily cleaned from any continuum signal using a running median filter of the same width. A linear χ 2 fitting of the first 50 Eigenspectra is independently performed for each spaxel in the cube to find the best linear coefficients of the Eigenspectra to match the sky line residuals. Here, the χ 2 is only computed in predefined wavelength regions covering bright sky lines and the best-fit linear combination is only subtracted in those wavelength regions. Afterwards the previously subtracted continuum signal from the median filtering is added back to the A&A proofs: manuscript no. CARS_presentation_v3 data. We highlight the resulting significant improvements of the sky line residual suppression with CubePCA in Fig. 3.

VIsible Multi-Object Spectrograph
Archival  VIMOS. The data were reduced with the Py3D data reduction package, which was initially developed for the CALIFA survey (Sánchez et al. 2012;Husemann et al. 2013a) and successfully applied to this VIMOS data set as described in Husemann et al. (2014), where more details on the data reduction can be found. Briefly, the reduction process includes basic steps such as bias subtraction, fiber tracing, flexure correction, wavelength calibration, flat fielding, and flux calibration. Since the spectral resolution significantly varies in dispersion and cross-dispersion direction across the detector, we first characterized the variation based on the arc lamp exposure and then adaptively smooth the data to a common spectral resolution of 3 Å (FWHM). Since no dedicated sky frames were taken, a mean sky spectrum was created from the object free regions and subtracted from the individually calibrated exposures. Given the size of HE 1338−1423 this leads to a some over-subtraction of the stellar continuum but not of the emission lines. Each observation is mapped to a cube so that the position of the AGN is traced as a function of wavelength. All exposures are then drizzled (Fruchter & Hook 2002) to a common grid using the AGN position as a reference for each wavelength to align the different exposures.

Potsdam multi-aperture spectrophotometer
As the MUSE data was only a filler program and not fully completed, we observed the two missing sources, HE 0853−0126 and HE 0949−0122, from the CARS sample with the Potsdam Multi-Aperture Spectrophotometer (PMAS, Roth et al. 2005;Kelz et al. 2006) at the 3.5m telescope of the Calar Alto Observatory under program H18-3.5-010 in December 2018. We used the PMAS lens-array with a 16 × 16 FoV and 1 × 1 spaxel size sufficient to map nearly the entire galaxies. The employed V600 grating covers the rest-frame optical wavelength range from Hβ to [S ii] at a spectral resolution of R ∼ 1500.
Total integration times were 2400 s split-up into two exposures bracketed by three 300 s long dedicated sky exposures for both HE 0853−0126 and and HE 0949−0122. Continuum and arc lamps calibrations were taken together after or before the target observations sky flats and standard star were taken in twilight for the matching spectral setups. We also used Py3D to reduce this data set following the same basic reduction steps as for the VIMOS data. Here, we adaptively smooth the spectral resolution also to a common value of 3 Å across the entire field and along wavelengths.
For background subtraction, we measure the sky background as the mean of all spectra in the dedicated sky exposures and subtracted it from the science exposures closest in time. The skysubtracted cubes were then resampled to a cube in order to characterize the differential atmospheric refraction by tracing the position of the AGN position as a function of wavelength. Based on the positions measured for each observation at each wavelength, we resampled the data into one cube following the drizzle algorithm (Fruchter & Hook 2002) similar to the VIMOS data.

Galactic extinction and absolute photometry
Most of the observations are not conducted under photometric conditions and standard stars are not always obtained close in time, in particular for the ESO service mode observation. Hence, it is important to infer the absolute photometric correction factor for a given observations. We use calibrated Pan-STARRS DR2 (Chambers et al. 2016;Magnier et al. 2020) and SkyMapper DR2 (Onken et al. 2019) cut-out r and i band images as a reference for our absolute photometric reference system. We construct corresponding broad-band images from the IFU cubes and compare the AB magnitude of stars in the IFU FoV with the reference images. In cases where no bright stars are captured in the IFU FoV, we use bright off-nuclear galaxy components or the integrated galaxy to minimize the impact of AGN variability. In particular for the changing-look AGN HE 0203−0031 and HE 0212−0059 it is essential to exclude the nucleus in the comparison apertures. The derived photometric scale factor Table 2. All reported fluxes in this paper are therefore divided by this factor to ensure a common photometric scaling across the survey with an absolute photometric uncertainty of <10 per cent as inherited from the reference surveys.
Foreground Galactic extinction can significantly reduce the observed flux and alter the overall shape of the spectra recorded for our extra-galactic targets. We therefore correct all MUSE, VIMOS and PMAS data cubes from Galactic extinction by dividing with the Cardelli et al. (1989) Milky Way optical extinction curve as a final step. The extinction curve is scaled to the line-of-sight V-band extinction as reported by the NASA/IPAC Extragalactic Database (NED) which is based on the far-IR maps presented by Schlegel et al. (1998) or SDSS stars (Schlafly & Finkbeiner 2011).

AGN-host deblending
The bright point-like emission from the AGN is able to dominate the light out to a few arcsec in ground-based seeing-limited observation. It is therefore important to subtract its contribution from each spaxel before the host galaxy extended emission can be properly analyzed. Here, we employ the AGN-host deblending method described in Husemann et al. (2013b) and Husemann et al. (2014), but expand the algorithm to deal with the wavelength-dependent PSF.
The AGN-host deblending process consists of five steps, 1) estimation of the PSF from broad lines at different wavelengths using QDeblend 3D (Husemann et al. 2013b, available at https://git.io/qdeblend3d), 2) modeling of the PSFs with a 2D Moffat profile to suppress noise at large radii, 3) interpolating the PSF as a function of wavelength, 4) reconstructing the intrinsic host galaxy surface brightness profile from 2D image modeling, and 5) applying an iterative AGN-host deblending scheme combining the wavelength-dependent PSF with the host galaxy surface brightness profile. In the following subsections, we describe each of the steps in more details.

Creation of PSFs from broad emission lines
The unobscured AGN selection for CARS has the great advantage that the PSF for a given observations can be empirically determined from the intensity distribution of broad emission lines emitted by the AGN broad-line region (BLR). The size of the BLR is known to be a few light months at most based on reverberation mapping studies (e.g., Kaspi et al. 2000;Peterson et al. 2004;Bentz et al. 2009b). Considering the distance to the galaxies, the BLR is inherently spatially unresolved in seeing-limited optical data and represents a perfect point source for our purposes. As initially outlined by Jahnke et al. (2004b), IFU observations are well suited to obtain the PSFs for unobscured AGN by constructing narrow-band images centered on the broad line wings from which the continuum emission is subtracted based on narrow-band images from the adjacent continuum. Similar approaches have been widely used to create PSFs from IFU data to properly recover extended emission around luminous AGN (e.g., Christensen et al. 2006;Husemann et al. 2008;Herenz et al. 2015;Borisova et al. 2016;Cantalupo et al. 2019;Drake et al. 2019).
Based on the wavelength coverage of MUSE and the redshift range of our targets, we determined empirical PSFs from the broad Hβ, Hα and O i λ8446+Ca ii λ8498 (see Matsuoka et al. 2007) emission lines using the QDeblend 3D (Husemann The wavelength coverage of the PMAS observation only allows us to construct the empirical PSFs from the Hβ and Hα lines. The VIMOS observations are taken with two spectral setups, of which the HR blue grism covers the broad Hγ, Hβ and He i lines while the HR orange grism only covers the broad He i and Hα lines.

Modelling each PSF with a 2D Moffat
Although estimating empirical PSFs is easy in IFU observations of unobscured AGN, they are not free from measurement uncertainties. The confusion with the host galaxy spectrum dominates at radii a few times the FWHM of the seeing disk, depending on the brightness of the AGN with respect to the host galaxy. Applying the pure empirical PSF would therefore significantly degrade the S/N in the host galaxy due to the high uncertainties in the light profiles in the wings of the PSF. To suppress the degradation of S/N in the empirical PSF as a function of radius, we model the PSF with a 2D Moffat function which is a good representation of the PSF profile in seeing-limited observations (e.g., Moffat 1969;Racine 1996;Trujillo et al. 2001;Jahnke et al. 2004a;Gadotti 2008).
While we could in principle use the best-fit 2D Moffat as our PSF profile, the empirical PSF is superior close to the AGN position where the most accurate PSF determination is required. Hence, we create a hybrid PSF as a combination of the modeled and empirical PSF, where the modeled PSF is replaced by the empirical PSF within the a certain radius of the peak position. The replacing radius is individually set for each PSF at the point where noise starts to dominate.

Interpolating the PSFs with wavelengths
The PSF is wavelength dependent as can be seen from Fig. 4. Therefore, we need to interpolate the 2D PSFs across the entire wavelength range. Here, we are limited in precision by the sparse sampling of PSFs in wavelength space for all IFU data sets. However, we obtained PSFs close to the blue and red end of the covered wavelength range, so that no significant extrapolation of the measured PSFs is necessary. We use a very simple but effective approach to interpolate the PSFs along the wavelengths. After normalizing each hybrid PSF to a peak flux of 1, we describe the wavelength-dependence of the normalized PSF flux for each spaxel with a simple polynomial function. The order of the polynomial is set such that it passes through the provided data points, so that we use we use a 2nd order polynomial for three PSFs. This ensures that the PSF cube is equal to the actual measurements at the wavelength where the hybrid PSF are provided.
For all MUSE targets with broad lines we used the nominal three PSFs with the only exception of HE 0021−1810. For this target the broad O i+Ca ii blended line is too weak to allow a PSF determination, so that we only used the PSF from Hβ and Hα PSF to interpolated them with a 1st order polynomial. This limits the usable wavelength range for the AGN-host galaxy separation. For the VIMOS data we can use three PSFs for the HR blue grism, but only two PSFs for the HR orange grism. The PMAS data allows us to use three PSFs distributed across the wavelength range for the interpolation.

2D surface brightness modeling of the host
For the next step in the iterative IFU AGN-host deblending, we need to create a surface brightness model of the host galaxy to avoid over-subtraction of the host galaxy. We therefore created i-band images from the MUSE cube itself. A corresponding PSF image is created from the interpolated PSF cube. We show two examples of the re-constructed broad-band images for HE 0108−4743 and HE 0351+0240 in Fig. 5. For the PMAS and VIMOS targets we use the Pan-STARRS i band images as the sensitivity, sampling and FoV is limiting the quality of the re-constructed images for those instruments. Similarly, we prefer to use the Pan-STARRS image also for HE 1237−0505 as the galaxy extends significantly beyond the MUSE FoV.
We model the re-constructed broad-band image for all galaxies with a PSF-convolved 2D surface brightness model using galfit (v3.0.5) (Peng et al. 2010a). Generally, two Sersic components plus a point source for the AGN leads to a sufficient model of the host galaxy for our purpose. While a more sophisticated modeling has been performed for some of the CARS galaxies to investigate their bar properties in more detail (Neumann et al. 2019), we only need to infer the central surface brightness profile of the host galaxy within the central 1 for the IFU AGN-host galaxy deblending. Here, we manually masked out nearby stars, background galaxies and interacting companions to keep the input model as simple as possible.
The best-fit model and residuals are shown for HE 0108−4743 and HE 0351+0240 in Fig. 5 and we list the inferred model parameters for a single and double Sersić model for all galaxies in Table 3. Best-fit 2D surface brightness model determined with galfit, which consists of a PSF for the AGN and two Sersić components for the host galaxy. Bottom panels: Residuals of the best-fit model which highlights the substructure of the galaxies such a spiral arms, but also foreground stars, companions and background galaxies are visible. Those features are masked out during the fitting as shown by the shaded areas.

Iterative AGN-host deblending of the IFU data
The last step is the actual AGN-host deblending which creates an AGN cube and a host galaxy cube. In principle, a host galaxy cube can already be created by subtracting the previously determined wavelength-dependent PSF cube scaled to match the actual AGN spectrum at the brightest, highest S/N, spaxel. However, even the brightest spaxel at the position of the AGN will include host galaxy emission so that some over-subtraction will occurs. As outlined in Husemann et al. (2014), we therefore perform an iterative AGN subtraction where the AGN spectrum is corrected for host galaxy contribution. We briefly recap the iterative part of the algorithm in the following.
The first iteration starts with a pure PSF subtraction as outlined above. Here, the PSF cube is scaled to the average of the central 3 × 3 spaxels in the MUSE cubes. After subtracting the scaled PSF cube from the original data a host galaxy spectrum is created from a rectangular shell surrounding the previous AGN spectrum region with a 2 spaxels width. For the PMAS and VI-MOS data we use the brightest spaxel for the AGN region and the surrounding spaxels as the host galaxy rectangular shell region to be consistent given the coarser spatial sampling of those data sets. The AGN-subtracted host galaxy spectrum is then scaled up in surface brightness toward the central AGN spectrum extraction region based on the 2D surface brightness of the host as determined from the MUSE data (see Sect. 4.1.4). Thereby, we implicitly assume that the circumnuclear host galaxy spectrum is the same as the center of the galaxy and can be subtracted from the initial AGN spectrum used for scaling the PSF cube. The underlying assumption here is that the central part of the galaxy is more extended than a point source and its spectral signature will be visible after point-source subtraction. The iterative process is repeated until the AGN and host galaxy spectrum have converged to a stable solution. We have consistently chosen 6 iterations for all IFU data sets after we confirmed a convergence to better than 1 per cent toward the final solution. Of course, the quality of the iterative process depends on the physical sampling and resolution. Any stellar population feature that appears pointlike as the AGN will still be associated with the point-like AGN, but also does not lead to further over-subtraction across the resolved galaxy which is the main point of the iterative approach.
The result of the iterative algorithm is shown in Fig. 6 for two representative targets observed with MUSE. The central 3 aperture spectra are shown for the total light and separated into the deblended host galaxy (extended) and AGN (point-like) contribution. Only the AGN spectrum contains broad emission lines, while the absorption lines from the stellar continuum are associated with the host galaxy spectrum. Although narrow emission lines are certainly extended, they are visible in both spectra, because part of the NLR originates from pc scales and therefore appear point-like in our seeing-limited observations. The narrow emission line contribution is thus shared between both components depending on the flux ratio of the unresolved versus resolved emission. Furthermore, the 2D surface brightness distribution in the iterative algorithm has been determined from the stellar continuum which may not necessarily apply also for the 2D surface brightness distribution of emission line flux. The iterative AGN-host galaxy deblending provides only a first order subtraction of point-like emission of narrow emission lines and a more detailed analysis is required to map the line properties in the central <1 as performed by Singha et al. (2021). Given the angular size of tens of arcsec for the CARS AGN host galaxies, our AGN-host deblending algorithm is capable of subtracting any compact emission. It allows us a detailed 2D investigation of the stellar and ionized gas properties in the AGN host galaxy down to 100-600 pc from the nucleus across the CARS target redshift range.

Visual classification of galaxies
Due to the selection of our sample based solely on the presence of broad emission lines, the host galaxies of our targets exhibit a variety in morphologies from disks to ongoing major mergers. Visual classifications of the AGN host galaxy morphologies were determined independently by 9-10 CARS team members based on re-constructed MUSE and archival broad-band images. A visual morphology can be easily obtained even without AGN subtraction in ground-based seeing-limited observations due to the low-redshift of the targets. Here, we limit the visual classification to very basic parameters and do not try to apply the full Hubble sequence classification. We focused mainly on basic morphology, bar presence and environment of the galaxies. The classifiers were asked to decide whether a galaxy was 1) bulgedominated, disk-dominated or irregular, 2) unbarred, barred or uncertain bar presence, as well as whether they showed 3) the presence of tidal tails, nearby companions or if the galaxy appears isolated. The third category of classification are not necessarily mutually exclusive because a galaxy may appear isolated but exhibit tidal tails from a past interaction. We therefore combined the classifications into the two classes, interacting and isolated galaxy, where a decision for either tidal tails or companions led to an "interacting galaxy" classification and galaxies determined as isolated were classified as an "isolated galaxy". Since the detection of faint tidal features is strongly dependent on image depth, our classifications need to be considered relative across the sample. Because the classification images have sim-  Table 4 for each galaxy.
In all cases, except for HE 0429−0247, a classification was obtained by absolute majority of classifications. The resulting distribution of majority classifications is shown in Fig. 7 which reveals that the CARS sample is clearly dominated by diskdominated AGN host galaxies with a fraction of 74%. A minority of CARS targets is found in bulge-dominated or irregular host galaxies. Strongly disturbed systems are clearly not prevalent in the sample, which is unsurprising given that no excess of strongly disturbed AGN host galaxies have been reported in the literature (e.g., Cisternas et al. 2011;Schawinski et al. 2012;Mechtley et al. 2016;Villforth et al. 2017;Marian et al. 2019) at least in the AGN luminosity range of the CARS sample. About 40% of the targets show some signs of interactions including major mergers and galaxy pairs, but these are limited to close companions due to the limiting size of the images. A full environmental study with 1 Mpc around the CARS galaxies is in preparation, so we cannot identify potential wide pairs and loose group environments for which the host would still be classified as isolated by the current metric.
The overall bar fraction is ∼50% within the entire CARS sample and 64% among the disk-dominated galaxy subgroup. This fraction could even be slightly higher as a couple of diskdominated galaxies are too much inclined to identify bars and are counted as uncertain, in addition to most of the irregular galaxies. A similarly high bar fraction was reported for X-ray-selected low-redshift AGN by Cisternas et al. (2015), while a bar fraction of 28.5% was reported by Alonso et al. (2013)  fraction and so this fraction needs to be interpreted with caution. The bar fraction of the overall low-redshift galaxy population is a strong function of color, stellar mass and bulge-to-disk ratio (Masters et al. 2011) ranging from ∼10% to ∼60% with a mean fraction of 29.4%±0.5%. A much higher bar fraction reaching ∼ 70% is typically found in IR imaging observations (e.g., Buta et al. 2015;Erwin 2018), which implies that the recovered bar fraction may be wavelength dependent. It is not the purpose of CARS to investigate the role of bars or interactions in triggering or fostering BH accretion in a statistical sense, but the presence of bars or strong gravitational interactions with companions can have a strong impact on the galaxy dynamics and distribution of gas on host galaxy scales. This needs to be carefully taken into account in the interpretation of the stellar and gas kinematics and distribution of star formation across the galaxies in terms of AGN outflows and AGN feedback. Hence, the visual classifications of the CARS sample presented here will be essential for all investigations performed with this sample.

Host galaxy sizes and axis ratios
The size and inclination of a galaxy with respect to our line of sight are important basic parameters for the CARS AGN hosts. For example, the impact of the AGN radiation on the host galaxy can be significantly affected by the orientation of central AGN engine with respect to the stellar disk (e.g., Husemann et al. 2019b). We estimated the effective radius R e and axis ratios through the 2D surface brightness modeling as part of the QSO-host galaxy deblending process, but those measurements can be influenced by a strong bar structure and are more reliable for the central part of the galaxy. We therefore performed an isophotal ellipse fitting additionally on the archival Pan-STARRS i-band images. Exceptions are HE 1108−2813, for which we use the r-band due to insufficient data quality in the i-band, HE 0108−4743 and HE 2211−3903, for which we use our dedicated wide-field imager (WFI) B-band images as presented in Smirnova-Pinchukova et al. (2021) given the lack of Pan-STARRS coverage.
The fitting is performed using the python package Photutils, an Astropy package for detection and photometry of astronomical sources (Bradley et al. 2019). We use the tools in photutils.isophote to fit ellipses to the isophotes in the images, which are measured using an iterative method described in Jedrzejewski (1987). The procedure determines for each ellipse the coordinates of its center, the semi-major axis, the ellipticity and the position angle; based on visually informed initial parameter guesses. It iteratively increases the sizes of the ellipses up to a predetermined surface brightness limit. In case of significant contamination of nonelliptical features, we use sigma-clipping, which especially improves our ability to fit in low signal-to-noise regions.
From the ellipse fitting we determine the semi-major axis and axis ratio at a limiting i-band surface brightness limit Σ lim = 24.5 mag arcsec −2 + 2.5 log((1 + z) 4 ) to correct for the surface brightness dimming across the redshift range of our sample. We apply an average color correction for the limiting brightness of r−i ∼ 0.25 mag and B−r ∼ 0.65 mag for star-forming galaxies in the three cases without i band images. This ensures that the host galaxy sizes of these three galaxies are comparable to the other ones. The corresponding measurements are listed in Table 4 for the entire sample.

IFU continuum and emission-line fitting
From the IFU data we want to map the stellar and ionized gas properties across the host galaxies. The process to subtract the point-like QSO emission is described above so that we can work with QSO-subtracted cubes at this stage. Here, we model the stellar continuum and ISM emission lines subsequently using PyParadise (see Walcher et al. 2015;Weaver et al. 2018;Husemann et al. 2019b). The subsequent fitting of stellar continuum and line emission ensures that complex line shapes do not impact the continuum modeling. While pPXF (Cappellari & Emsellem 2004) has become the standard for the stellar continuum modeling in the IFU analysis of nearby galaxies (e.g Westfall et al. 2019;Bittner et al. 2019;Croom et al. 2021), the PSF interpolation along wavelength for the QSO subtraction leads to significant low-frequency continuum variation close to the nucleus that cannot be easily described with low-order polynomials used by pPXF.
PyParadise instead normalizes the stellar continuum and template spectral library before fitting and is therefore much less sensitive to global unphysical continuum variations. While some parameters, such as line-of-sight extinction for the stellar continuum, cannot be determined this way, it allows us to obtain robust stellar kinematics and absorption line equivalent width much closer to the nucleus. The emission lines are fitted to the continuum residuals and we choose single Gaussians coupled in radial velocity and velocity dispersion for all emission lines as a zeroorder approximation. A coupling of the lines enhances the robustness of the flux measurements for fainter lines and increases the precision on global gas kinematics. Follow-up analysis in specific regions with complex line kinematics can be performed later on with dedicated methods (Husemann et al. 2019b). PyParadise is publicly available at https://git.io/pyparadise and more details on the algorithm can be found in the user manual. The S/N of individual spaxels in the IFU data can be too low to properly model the stellar continuum and low surfacebrightness emission line regions. We therefore employ two different binning strategies as described in the following.
For bright emission line regions we want to retain the native MUSE spatial sampling of 0 . 2 for which the S/N in the stellar continuum may not be sufficient for the continuum modeling. We therefore perform an initial Voronoi binning to achieve a mini-

HeII
[SIII] λ9096 [NI] λλ5200,5202 mum continuum S/N of 20 per Voronoi cell per spectral pixel at around 6000Å. We then model the binned stellar continuum spectra with PyParadise using stellar spectra from the INDO-US spectra library (Valdes et al. 2004) and remap the inferred stellar kinematics back to the native MUSE sampling grid. Afterwards we model the stellar continuum in the full cube while keeping the stellar kinematics fixed, so that the emission lines can be fitted to the residual continuum spectra for each individual MUSE spectra. Additionally, we bin the QSO-subtracted datacube by 8×8 pixels for MUSE, 2×2 pixels for PMAS, and 4×4 pixels for VIMOS corresponding to a binned sampling of 1 . 6, 2 . 0, and 2 . 64, respectively. A full continuum and emission-line modeling with PyParadise is performed on the binned cubes without the intermediate Voronoi binning step for the stellar continuum. The analysis of the binned data allows us to trace the emission lines nearly an order of magnitude fainter in surface brightness than in the unbinned data.
An example of a single spaxel fitting with PyParadise is shown in Fig. 8. All masked regions for the emission-line modeling are indicated as gray shaded areas, whereas regions of strong atmospheric absorption and sky line residuals are highlighted as red shaded areas. We also label all emission lines that we incorporate in the automatic modeling. An important part of this process is to infer proper errors on all measured parameters to identify robust measurements and to set upper limits of undetected emission lines for diagnostics purposes. Errors on the stellar kinematics are inferred using an MCMC approach as part of PyParadise iterative continuum modeling. The errors on all emission-line parameters are inferred combining a bootstrapping and Monte-Carlo approach. The entire modeling of the stellar continuum (at fixed stellar kinematics) and emission-line modeling are repeated 50 times after randomly modifying each spectral pixel within its noise and randomly reducing the number of input template stellar library spectra to 60%. This way we incorporate some systematic uncertainty of the stellar continuum modeling also in the derived emission line parameters. This is most important for emission lines that are superimposed on stellar absorption lines, mainly Hβ and Hα.
An example of the inferred 2D kinematics and emission-line distribution is shown in Fig. 9. For the emission-line surface brightness maps we combined the results from the unbinned and binned QSO-subtracted data. Assuming a minimum line detection of S/N>5, we replace undetected regions in the unbinned data with the surface brightness of the binned data if the S/N limits were achieved for a given binned spaxel.

Ionized gas excitation conditions and ENLR properties
It is a common approach to map the excitation conditions of the ionized gas across the AGN host galaxies. The Baldwin-Phillips-Terlevich (BPT, Baldwin et al. 1981) diagram comparing the [O iii]/Hβ versus [N ii]/Hα line ratios has been extensively used at optical wavelengths for this task (e.g., Veilleux & Osterbrock 1987;Kewley et al. 2001;Kauffmann et al. 2003;Kewley et al. 2006;Cid Fernandes et al. 2011). We show the overall PyParadise results on stellar kinematics and emission line properties for the CARS data of HE 0853+0102 in Fig. 9 together with the derived BPT diagram and associated classification. different sources of ionized gas excitation, such as ionization by young stars in star-forming sides often referred to as H ii regions, the photoionization by the hard radiation field of an AGN and low-ionization nuclear emission-line regions (LINERs). Clearly distinguishing between such mechanisms is not trivial, but several demarcations lines have been proposed to discriminate between those sources. Given the redshift of our targets, the spatial resolution is limited to a few 100 pc, which does not allow to see individual line emitting clouds. Mixing of different ionization mechanisms along the line-of-sight is an additional complication. This is evident in Fig. 9 where the majority of spaxels are consistent with ionization by star forming regions, which has a smooth extension toward the AGN location in the BPT. The 2D spatial mapping of the excitation confirms that the higher ionization is found closer to the AGN. This is known as the SF-AGN mixing-sequence (e.g., Davies et al. 2014;Richardson et al. 2014;Davies et al. 2016) and is important to take into account when inferring the current star formation using the Hα line powered solely by H ii regions. Such a detailed analysis will be presented in the companion paper (Smirnova-Pinchukova et al. 2021) and has been applied for a subsample of galaxies (Neumann et al. 2019). Here, we focus our investigation on the characterization of clouds significantly influenced by the AGN photoionization, which can be unambiguously associated with the ENLR.
From the excitation maps as shown in Fig. 9, we determine the maximum distance of AGN-ionized regions (R ENLR,max ) with respect to the AGN location that we can cover with our depth and area. To suppress the impact of noise on misclassification close to the borders of the demarcation lines, we require that at least six surrounding unbinned spaxels share the same AGNionization classification or the classification is present in the spatially binned data with its intrinsically higher S/N. We need to exclude the LINER and the intermediate BPT regions, because those are not necessarily associated with AGN photoionization and can originate from post-AGN stars (e.g., Binette et al. 1994;Singh et al. 2013), shocked gas from stellar winds (e.g., Ho et al. 2016;López-Cobá et al. 2019), or the diffuse ionized gas inbetween star forming regions (e.g., Lacerda et al. 2018;Levy et al. 2019). We define the maximum ENLR size, R ENLR,max , as the maximum projected distance of a robust AGN region to the AGN location detectable within the IFU FoV at the observational depth. In general, the maximum ENLR extent can be biased due to surface brightness dimming and observational depth. This is not a concern for the large FoV of MUSE and the narrow redshift distribution of the CARS sample. Nevertheless, we also determine two additional ENLR radii, R ENLR,15 and R ENLR,16 , corresponding to the ENLR size out to a surface brightness limit of Σ 15 = 10 −15 erg s −1 cm −2 arcsec −2 (1 + z) −4 (e.g Liu et al. 2013a;Hainline et al. 2013;Liu et al. 2014;Hainline et al. 2014), or Σ 16 = 10 −16 erg s −1 cm −2 arcsec −2 (1 + z) −4 (Chen et al. 2019a), respectively. All these measurements are listed in Table 5. The ENLR size can become unresolved for the CARS targets, in particular, at high intrinsic ENLR surface brightnesses. For such nondetections we set the upper limits for the ENLR size to be the angular size of the binned spaxels. The error on the sizes are generally driven by the spatial resolution of a given observations, so that we assume the error on all sizes to be the FWHM of seeing for each observation as reported in Table 2.
As the blue part of the optical spectrum below <4700Å is not covered by the MUSE data, a full spectrum fitting of the optical AGN spectra usually performed for large AGN surveys (e.g., Shen et al. 2008;Park et al. 2012;Coffey et al. 2019) is not ideal in our case. In particular, the use of Fe ii templates (e.g., Boroson & Green 1992;Véron-Cetty et al. 2004;Kovačević et al. 2010) often leads to bad fits for strong Fe ii emitters in our high S/N MUSE spectra and quality was deemed insufficient for several individual cases. Furthermore, the Fe ii blends on the blue side of Hβ are not covered in most of our IFU observations which reduces the robustness of the template fitting. Therefore, we decided to model only the rest-frame wavelength from 4750Å to 5100Å and adopt a linear pseudo-continuum for this narrow wavelength range as described in Singha et al. (2021). For the broad and narrow emission lines we adopted super-positions of Gaussian line profiles. More specifically we adopted one or two Gaussians for the broad Hβ line and the two prominent isolated Fe ii λλ4923, 5018 lines originating both from the broad-line region of the AGN. The narrow lines of [O iii] λλ4960, 5007 and Hβ are commonly modeled as two Gaussian components to account for the typical asymmetry of those lines in AGN (e.g., Greene & Ho 2005a;Mullaney et al. 2013;Harrison et al. 2016;Bennert et al. 2018). To keep the modeling robust and reduce the number of free parameters, all broad and narrow line components are coupled in redshift and kinematics across line transitions. Additionally, the [O iii] doublet line ratio is constrained to be 3 (Storey & Zeippen 2000) and we fix the flux ratio Fe ii λ4923/Fe ii λ5018 = 0.81, which is the mean of the empirically measured ratios when leaving the fluxes unconstrained.
An example of the modeling is shown in Fig. 10 and entire sample fits are presented in Singha et al. (2021), where we focus on the line shape of [O iii]. For the purpose of this paper, we are only interested in measurements that are relevant to estimate bolometric luminosities, BH mass, etc., so we list only the total line fluxes of the broad Hβ and Fe ii lines together with their line widths as well as the total [O iii] flux in Table 6. Errors are determined through a Monte Carlo approach as described in Singha et al. (2021), but a 10% systematic uncertainty is added on the line fluxes to take into account the absolute photometric calibration uncertainty of the IFU observations. A&A proofs: manuscript no. CARS_presentation_v3

AGN parameter estimation
The main AGN parameters we are interested in are the bolometric luminosity, the BH mass and the Eddington ratio which are easy to obtain from the unobscured AGN spectra. Bolometric correction factors have been estimated for the AGN continuum light at 5100Å such as L bol ∼ (8-12) × λL 5100 (e.g., Kaspi et al. (1) from which we computed L bol , listed in Table 6, adopting a bolometric correction factor of 10 for L 5100 consistent with Richards et al. (2006).
Single-epoch BH masses are based on the empirical BLR size-luminosity relation and a virial factor f that captures assumptions on the detailed BLR geometry in order to convert the observed broad-line width into a gravitational potential. Different size-luminosity relations (e.g., Kaspi et al. 2000;Bentz et al. 2009aBentz et al. , 2013, virial factors (e.g., Onken et al. 2004;Woo et al. 2013;Pancoast et al. 2014), and BLR line width definitions, such as FWHM versus line dispersion (e.g., Collin et al. 2006), have been used in the literature which makes it hard to compare different samples. We adopt the relation of Greene & Ho (2005b) to compute BH mass which provides one specific single-epoch BH mass estimate. In Table 6, we provide all measurements necessary to re-calculate the BH mass using different calibrations when needed for comparison purposes. Beside the measurement errors, a systematic uncertainty of 0.3 dex is typically assumed for a single-epoch BH mass estimate of a single object. We also compute the corresponding Eddington luminosity as L Edd /[erg s −1 ] = 1.26 × 10 38 M BH /[M ] and define the Eddington ratio λ = L bol /L Edd .

AGN and Eigenvector 1 parameter space
One of the major reasons to focus on type 1 AGN for CARS is that the primary AGN parameters such as M BH , L bol and λ can be directly inferred from the AGN spectra, which are often difficult to constrain for type 2 AGN without additional assumptions. One intrinsic correlation of AGN parameters has been the famous Eigenvector 1 parameter space (e.g., Boroson & Green 1992;Sulentic et al. 2000), which shows that AGN spectra exhibit a striking correlation between [O iii] peak height and Hβ FWHM as well as an anticorrelation between the Fe ii band strength and Hβ FWHM. These trends have been confirmed to be tightly correlated with the Eddington ratio of the AGN as the FWHM has a much stronger impact in the BH mass determination than the luminosity. Therefore, we explore whether the CARS sample is representative also across this important parameter space.
Usually the Fe ii band flux is measured from the blue side of Hβ between 4434-4685Å, but is not covered by MUSE. We measure the Fe ii flux by fitting the wavelength region between 5070-5520Åwith the Fe ii templates from Kovačević et al. (2010). A first-order polynomial continuum sampled at both sides of the Fe ii blend is subtracted beforehand and the emission lines [N i] λ5199, [Fe xiv] λ5302 and He i 5411 are masked to avoid a potential contamination. During the fit, the individual groups defined in the Fe ii templates are scaled in flux per group, while all groups are forced to share the same line-of-sight velocity and line width. The total Fe ii flux was derived by integrating the resulting model between 5100Å and 5500Å. The measurement was repeated 100 times after randomly modulating the input spectrum with the error spectrum. The Fe ii flux values and errors are reported in Table 6 which correspond to the mean and standard deviation of the 100 measurements. This needs to be converted to the Fe ii(4434-4685Å) flux in order to compare with literature data. The mean ratio between the bands was measured to be Fe ii(4434−4685Å) Fe ii(5100−5550Å) = 1.44±0.55 by Kovačević et al. (2010). We scale our flux accordingly and present the FWHM of Hβ against R Fe = Fe ii(4434−4685Å) Hβ for the CARS sample in the upper panel of Fig. 11. A classification of the AGN into a population A (Hβ FWHM < 4000 km s −1 ) and B (Hβ FWHM > 4000 km s −1 ) has been introduced by Sulentic et al. (2000) based on radio properties and further refined by Sulentic et al. (2002) into several subcategories.
We compare the distribution of CARS with AGN from the SPIDERS (SPectroscopic IDentification of eROSITA Sources Dwelly et al. 2017) survey, which is an SDSS-IV wide-area spectroscopic follow-up program of X-ray selected source from ROSAT and XMM-Newton. The modeling of the spectra and resulting catalog of parameters are reported in Coffey et al. (2019). About ∼ 170 AGN from their catalog are at z < 0.06 consistent with our CARS sample selection. According to Fig. 11, the CARS sample is lacking extreme iron emitters at R Fe > 1.5. Nevertheless, the trend of increasing Eddington ratio as FWHM is decreasing and R Fe is increasing can be recovered. The bolometric luminosity against BH mass is presented in the lower panel of Fig. 11 for CARS and SPIDERS. Again, the distributions are comparable, but CARS contains a few more higher luminosity AGN due to the larger effective area and volume covered by HES compared to the X-ray surveys. Both samples contain BH masses, AGN luminosities and Eddington ratios in the range of 6.0 < log(M BH /[M ]) < 8.3, 43 < log(L bol /[erg s −1 ]) < 45.5 and −1.5 < log(λ) < 0.5, respectively. Hence, the CARS subsample of local AGN is representative for the population in all important AGN parameters independent of the optical versus Xray selection of type 1 AGN.

The ENLR size-luminosity relation
A relation between the ENLR size and the AGN luminosity has been observed in many studies based on imaging (e.g., Bennert et al. 2002;Fischer et al. 2018;Sun et al. 2018), long-slit spectroscopy (e.g., Greene et al. 2012;Hainline et al. 2013Hainline et al. , 2014Sun et al. 2017) and IFU observations (e.g., Liu et al. 2013a;Husemann et al. 2014;Kang & Woo 2018). While HST narrowband imaging offers the highest angular resolution, the limited sensitivity of HST for low surface brightness features can lead to an underestimation of the ENLR extent (e.g Greene et al. 2011;Husemann et al. 2013b). This leads to the primary issue of defining an objective ENLR size criterion based on the [O iii] light distribution. Like expensive multi narrow-band imaging, spectroscopic studies also have the advantage that the ionization mechanism of the gas can be verified to be AGN ionization as required for the ENLR. In Fig. 12 We assume that the ENLR sizes follow a log-Normal probability distribution function at a given luminosity: .
( 6) The intrinsic scatter of the relations is quite large with σ ENLR of 0.42±0.06 dex, 0.47±0.07 dex, 0.51±0.08 dex, for R ENLR,max , R ENLR,16 , and R ENLR,15 , respectively. In particular, we find numerous AGN host galaxies which display rather large ENLRs exceeding 10 kpc below L AGN < 10 45 erg s −1 which are similar to the sizes of significantly more luminous QSOs in the literature even for R ENLR,15 . One important aspect, which has not been discussed in the literature so far, is the implications of limited area and spatial resolution for IFU observations. This leads to the fact that the maximum and minimum ENLR size is significantly changing with target distance and redshift as shown in Fig. 13. These ranges also vary substantially between IFU instruments and observing site characteristics. The WFM of MUSE has by far the largest dynamical range for recovering the ENLR by area and likely also by depth compared to all other current IFU instruments. This might introduce tighter correlations and biased slopes as more luminous AGN are naturally targeted at higher redshifts than the lower redshift counter parts.
The systematic offset toward smaller ENLR size is evident in the detected ENLRs for unobscured QSOs observed with the PMAS IFU instrument (Husemann et al. 2013b) shown in the left panel of Fig. 12. This is expected for the small PMAS FoV compared to the big MUSE FoV when trying to recover R ENLR,max . Using the intrinsic high surface brightness limits for the ENLR size definition makes the ENLR size generally smaller so that the resolution limit of the observations becomes more important than the recoverable maximum size. This leads to the issue that the sizes will cluster close to the resolution limit and an increased fraction of ENLRs become actually unresolved as seen in our CARS data and the MaNGA data of Chen et al. (2019b).
Combining different samples from various instruments at different redshifts therefore inevitably introduces ENLR sizeluminosity relations with different slopes α depending on the details of target selection and analysis approaches. Slopes ranging from α = 0.22 ± 0.04 (Greene et al. 2012), α = 0.25 ± 0.02 (Liu et al. 2013a), α ∼ 0.3-0.4 (Hainline et al. 2013;Chen et al. 2019b), to α ∼ 0.5 (Bennert et al. 2002;Husemann et al. 2014) are reported in the literature. The slopes solely inferred from the CARS data are consistent with those reported by Greene et al. (2012) and Liu et al. (2013a) and are therefore on the shallower side of previous estimates. Nevertheless, the scatter in the observed relation is significant and measured slope variations might be entirely attributed to the observationally induced biases as discussed above. A slope of α = 0.5 is reminiscent of the BLR size-luminosity relation, but would require a constant ionization parameter U that demands a constant density with radius. This is not observed for the ENLR on kiloparsec scales (e.g., Bennert et al. 2006;Kakkad et al. 2018) and more detailed photoionization calculations are required to predict the shallower slopes inferred for most studies (Dempsey & Zakamska 2018). We cannot study the radial variations of U as our snapshot MUSE observations are not deep enough to map the electron density given the too low S/N of the [S ii] doublet on kpc scales. However, the photoionization calculations do not take into account variable ionizing flux from AGN on 10 5 yr time scales (Schawinski et al. 2015) and the various geometrical intersections of the ionizing radiation field with the gas distribution of the galaxies. The CARS survey is least biased with regard to R ENLR given the narrow redshift range and large dynamic range offered by MUSE (see Fig. 13). Therefore, the CARS survey is one of the best data set to explore the origin of the significant scatter in ENLR sizeluminosity relation and search for additional factors or more fundamental parameters controlling the ENLR size.
5.2. BH mass as a more fundamental driver for the ENLR size Husemann et al. (2008) already reported that the AGN luminosity cannot be the only parameter setting the ENLR size or its radial surface brightness distribution on kpc scales. They found that kpc-size ENLRs were preferentially present around AGN with Hβ FWHM larger than 4000 km/s and low Eddington ratios at a given AGN luminosity. Furthermore, large ENLR were preferentially detected around radio-loud AGN rather than radioquiet AGN (e.g., Stockton & MacKenty 1987). We can perform a more extensive study with CARS by looking at a variety of host galaxy and AGN parameters that may be important to control the ENLR size and represent more fundamental correlations. In Fig. 14, we show the Pearson correlation matrix between the three ENLR sizes, redshift, AGN parameters, and various host galaxy parameters as determined by Smirnova-Pinchukova et al. (2021). We find that the ENLR sizes are strongly correlated with each other (r > 0.5 nd p < 2 × 10 −4 in all cases) as naturally expected. Similar internal correlations are found between the AGN parameters and the host galaxy parameters. We find only weak correlations between the ENLR sizes and redshift or AGN luminosity with r < 0.45 and p > 0.02 in all cases. Strikingly, we find the strongest correlation between the maximum ENLR size R ENLR,max and the BH mass M BH with a correlation coefficient of r ENLR,max = 0.71 and p = 5 × 10 −6 . Because single-epoch BH mass estimates are computed based on AGN luminosity and the FWHM of the Hβ line (e.g., Kaspi et al. 2000;Vestergaard 2002;Peterson et al. 2004), it explains that 1) the ENLR shows a weak correlation with AGN luminos-lo g r E N L R , m a x lo g r E N L R , 1 0 1 6 lo g r E N L R , 1 0 1 5  ity and 2) the detection of the ENLR also depends on the FWHM of Hβ at fixed luminosity (Husemann et al. 2008). We plot R ENLR,max as a function of M BH in Fig. 15 for the CARS sample excluding targets observed with VIMOS and PMAS to avoid observational biases due to their smaller FoV. Following the same Bayesian approach as for the ENLR sizeluminosity relation we infer a scaling relation of R ENLR,max pc = 10 (3.81±0.06) M BH 10 7 M (0.50±0.10) with an intrinsic scatter of σ ENLR = 0.32 ± 0.04 dex. This intrinsic scatter is significantly smaller compared to a scatter of 0.42 ± 0.06 dex for the linear relation with L bol . Hence, M BH can be considered a more fundamental driver for the maximum ENLR size. As M BH is an inferred quantity derived from AGN luminosity and the FWHM of broad lines, we also infer a bivariate relation based directly on the actual observables which are the width and luminosity of the broad Hβ component: with an intrinsic scatter of σ ENLR = 0.31 ± 0.04 dex similar to the one using the M BH . As FWHM and Hβ luminosity are directly observed quantities, they are independent of the various prescriptions used to compute BH masses. It is also possible that the apparent correlation with BH mass is emerging as a result of the correlation with AGN luminosity and FWHM which depends on the actual physical cause of such relations.
It is important to note that the relations of ENLR size as a function of BH mass and the BLR parameters are only applicable to confirmed AGN with radiative-efficient accretion and type 1 AGN host galaxies, respectively. One cannot predict an ENLR size around non-AGN galaxies or radio AGN. The CARS sample is specifically selected to host BH with ongoing radiativeefficient accretion starting roughly at 1 % Eddington except for the changing-look AGN in the sample. It is also possible that the absolute normalization of the R ENLR -M BH relation will have a dependence on additional AGN or host galaxy parameters, but these dependencies are not detectable with the still relatively small sample size of CARS. While the relation may suggest that the BH mass could be estimated from the ENLR size for obscured type 2 AGN, the individual AGN lifetimes as discussed below introduce significant systematic uncertainties for individual galaxies. Further studies are needed to quantify such systematics and test additional applications.

ENLR sizes of changing-look AGN
The time variability of AGN luminosity becomes dramatically visible in changing-look AGN, which change their luminosity by orders of magnitude within a couple of years. They highlight that BH accretion is not necessarily preserved over long time periods and current observations provide only a current snapshot of AGN properties. However, the BH mass is a nearly constant quantity with time even for changing-look AGN and smoothly increases with time. Luckily, the CARS sample contains 2 well-known changing-look AGN, namely HE 0212−0059 aka Mrk 590 (Osterbrock 1977;Denney et al. 2014) and HE 0203−0019 aka Mrk 1018 (Cohen et al. 1986;McElroy et al. 2016;Krumpe et al. 2017), which we use here to test the impact on the ENLR size scaling.
We include both changing-look AGN as stripes in Fig. 12 representing the range in AGN luminosity from the observed maximum to 1/100 and 1/10 of the luminosity as observed for Mrk 590 (Denney et al. 2014) and Mrk 1018(McElroy et al. 2016Krumpe et al. 2017), respectively, during their minimum state. It is clear that both sources can be significantly offset from the mean ENLR size-luminosity relation depending on the actual phase (minimum/maximum) of accretion rates. Such drastic changes might happen in any AGN and would only be noticed with a continuous monitoring program over decades. Therefore, we are not aware of the long-term variations and current accre- tion rate status of all other AGN in our sample, which might play a role in the scatter of the ENLR size-luminosity relation. As shown in Fig. 15, both changing-look AGN Mrk 1018 and Mrk 590 are in agreement with the ENLR size-BH mass relation. For Mrk 590 we use the broad Hβ parameters as measured from the MUSE IFU data during a recent rebrightening phase, and for Mrk 1018 we used the BH mass estimated during its bright phase (Bennert et al. 2015) as the broad Hβ line remained too faint after the last changing-look event. The position of both changing-look AGN supports the notion that the M BH is a more fundamental parameter in predicting the ENLR size and insensitive to short term AGN variability.

Is the ENLR size an AGN episode lifetime indicator?
A new key result of this study is that the ENLR size is more strongly linked with the BH mass than with AGN luminosity. A major question is whether this relation has a physical origin directly linked to the BH mass or if it is a secondary correlation given that the BH mass is known to be linked with host galaxy properties. Indeed, we also observe a correlation between BH mass and the host galaxy size even for our predominately disk-dominated systems (see Fig. 14). Such a link between BH mass and host size (R host ) has been reported in the literature (e.g., van den Bosch 2016) even for disk-dominated galaxies. Although R host and R ENLR are correlated, R ENLR shows a much larger spread than R host . The ENLR sizes are often smaller than the host galaxies, which suggests that the ENLRs are likely to be ionization-bounded and the availability of gas clouds is not a major limiting factor for the overall ENLR size within the limited BH mass and AGN luminosity range probed by the CARS sample. Also we do not see the break-down of the ENLR size at ∼10 kpc that was previously reported (Hainline et al. 2013(Hainline et al. , 2014Dempsey & Zakamska 2018). While the complex and unknown gas distribution and ionization conditions will certainly have an impact on the ENLR size and vary from object to object, they should be marginalized statistically by the CARS sample size and comparable host galaxy properties.
Another mechanism predicted to scale with the BH mass is the condensation radius of gas cooling out of the turbulent hot halo of galaxies and groups via the Chaotic Cold Accretion (CCA) scenario (e.g., Gaspari et al. 2018Gaspari et al. , 2019. Indeed, more massive BHs are hosted in more extended and more massive halos (e.g., Krumpe et al. 2015), which implies a larger rain of warm and cold gas. Cool gas with ≤ 10 4 K is a requirement to detect AGN-ionized gas clouds via optical emission lines, so this mechanism can set an upper bound for R ENLR . Furthermore, the projected ENLR size and FWHM of Hβ could both be similarly affected by the inclination of the AGN central engine and thereby introduce a correlation between the two quantities. A sin i factor can impact the observed Hβ FWHM if the BLR clouds have a disk-like distribution as reported from velocity-resolved reverberation mapping (e.g., Grier et al. 2013;Pancoast et al. 2014;Grier et al. 2017;Williams et al. 2018). However, the inclination distribution of a random type 1 AGN sample has a welldescribed form when a simple cone geometry is assumed with a fixed half-opening angle θ ENLR (see Appendix A). The fraction of unobscured X-ray selected AGN is ∼50% at low redshift for AGN luminosities matching CARS sample (Merloni et al. 2014), which corresponds to a half-opening angle of θ ENLR = 60 • for the AGN ionization cones. The average inclination of a randomly selected unobscured AGN sample would then be i ENLR = 39 • . The corresponding 1σ confidence interval of the inclination distribution is [22.4 • , 54.5 • ] where inclinations toward 0 • become more and more unlikely. Such a narrow inclination range is consistent with the BLR inclination distribution directly determined from forward modeling of velocity-resolved reverberation mapping data of a random type 1 AGN sample (Williams et al. 2018). The inclination of the central engine alone might therefore lead to a log(R ENLR ) ∝ log(FWHM Hβ ) relation and could potentially explain the link between a kpc-scale region and the kinematics of clouds at <pc distance. However, the projection of the ENLR size should be invariant to the small inclination range of type 1 AGN due to the much larger ionization cone angle, so that a correlation between R ENLR and FWHM Hβ would not be introduced by the central engine inclination. Our large ionization cone angle assumption is supported by the fact that we do not find a systematic difference in R ENLR size between our type 1 AGN sample of CARS and the type 2 AGN sample of Chen et al. (2019b) (see middle panel of Fig. 12). Based on that, we assume that R ENLR is insensitive to the central AGN engine inclination and does not introduce a direct correlation with Hβ FWHM. This argument would break down if the ionizing radiation is not isotropic, but would lead to other major implications because an intrinsically isotropic AGN radiation field is a major assumption in many applications.
It is instructive to consider another possible physical parameter that can link the properties of the central AGN with the radiation field on kpc scales, which is AGN variability. A common implicit assumption is a constant ionizing flux of the AGN over time which may break down for a given AGN luminosity at various timescales (e.g., Keel et al. 2017). The distance of an AGN-ionized gas cloud to the AGN engine directly translates into a light travel time. In this way, the ENLR size can be directly converted into a lower limit for the AGN lifetime if the BH is still accreting or an upper limit if the BH is not accreting anymore. Light echos of shutdown AGN (e.g., Lintott et al. 2009;Keel et al. 2012;Schweizer et al. 2013;Keel et al. 2017) suggest an AGN life time for a single accretion episode to be ∼ 10 5 yrs, which is consistent with the fraction of young AGN without strong NLR emission (Schawinski et al. 2015). In addition, the ionization of neutral hydrogen or helium around highredshift QSOs, the so-called proximity zones, can be used to determine current QSO "on time" (t on ) for an individual source which is in the range of 10 5 to 10 7 yr (e.g., Eilers et al. 2017;Worseck et al. 2021). Inferring the AGN lifetime (t AGN ) from the ENLR or proximity zone size requires a statistical analysis to account for additional effects, such as neutral gas fraction, or random sampling of t on for a specific AGN from the AGN lifetime distribution. One issue here is that AGN lifetime is not a well defined quantity as AGN brightness is known to vary on different timescales for different luminosities and accretion levels.
In the picture we develop here, the best representation of AGN lifetime would be the stability time-scale of a radiative-efficient accretion disk with L bol /L Edd 0.01.
Such a statistical framework has recently been developed and applied for the He ii proximity zones around luminous z ∼ 3 QSOs . They assume a log-normal distribution of AGN lifetimes with a mean t AGN and dispersion σ AGN as well as a light-bulb light curve (top-hat function in time) for the AGN luminosity. Taking into account also random fluctuations in the initial He ii fraction they create probability density functions (PDFs) for the proximity zone size across a grid of t AGN and σ AGN values. Those PDFs are determined by applying a simple kernel density estimation for a large simulated population. In this way, the PDFs are effectively marginalized over all important stochastic relationships. Using the estimated PDFs, the likelihood function for the proximity zone measurements can be computed as a function of t AGN and σ AGN and the posterior distribution function can be sampled with a MCMC sampler.
In particular, the light-bulb approximation is a strong simplification of AGN variability and more complex light curves will be more realistic. In order to be consistent with the definition of Khrykin et al. (2021) for t AGN and considering our sample size, we follow their simple light-bulb assumption to avoid further complexities. It should be noted that the absolute normalization of t AGN is dependent on this assumption. Instead of the variable He ii fraction we need to incorporate some average inclination of our type 1 AGN sample with respect to our line-of-sight to account for the ionization cone projection effect. Here, we adopt an average inclination i ENLR = 39 • as discussed above. Details on the PDF generation and the construction of the likelihood function are provided in with a dispersion of σ AGN = 0.09 +0.08 −0.05 . The relation is shown in Fig. 15 and represents an upper envelope to the actual ENLR sizes while retaining the slope. Strikingly, our extrapolated relation to M BH > 10 9 M is consistent with the inferred AGN lifetime of log(t AGN /[yr]) = 6.22 +0.22 −0.25 for high-mass BHs at z ∼ 3 ) based on their He ii proximity zones. Our BH mass-dependent model for the duration of a single AGN lifetime episode can simultaneously predict the 10 5 yr lifetime of local shutdown AGN as well as the 10 6 yr lifetimes of more massive BH at higher redshifts. It is unclear if this is just a coincidence or a signature of an underlying physical mechanisms that is independent of redshift and mainly dependent on the gravitational potential of the central BH. It is interesting to note that the CCA framework mentioned above predicts a similar relation between AGN lifetime and BH mass. The CCA-driven AGN activity roughly scales with the hot gas cooling time, t cool ∝ T X , and M BH ∝ T 2 X , hence t AGN,CCA ∝ M 0.5 BH . At low radiative AGN accretion efficiencies, powerful radio jets can emerge from an AGN from which a much more precise AGN lifetime can be estimated for individual systems (e.g., Biava et al. 2021). The implied AGN lifetimes are typically on the order of tens of Myr given the observed jet lengths of several hundreds to thousands of kpc (e.g., Hardcastle et al. 1998) moving close to the speed of light (>0.1c), which likely depends on the ram pressure balance with the ambient medium. While this is much longer than the t AGN ∼ 10 5 yr lifetime for the radiative-efficient accretion in the CARS AGN, which might partially related to the very different accretion mode of the central engine, the discrepancy in time scales is actually not that large when considering the BH mass dependence. The fraction of galaxies hosting powerful radio-jets seems to be a strong function of stellar mass (e.g., Best et al. 2005) which is presumably tightly correlated with BH mass for bulge-dominated galaxies (Marconi & Hunt 2003;Häring & Rix 2004) reaching beyond M BH ∼ 10 10 M (e.g., Magorrian et al. 1998;McConnell et al. 2011;Mehrgan et al. 2019). The extrapolation of our relation to such high BH masses would predict radiative-efficient AGN accretion lifetimes of several tens of Myr. In addition, Gigahertz Peak Spectrum (GPS) radio sources exhibit exceptionally compact jets with estimated expansion times of 10 4 yr and are therefore interpreted as young radio source (e.g., Vink et al. 2006;Czerny et al. 2009). Hence, there is also a big range in observed jet length that can only be explained in a time evolution picture.
Interpreting the ENLR size to be induced predominantly by long-term AGN variability, as proposed by Schawinski et al. (2015), is therefore a possible scenario. Indeed, variability studies of AGN provide direct evidence that short-term and longterm X-ray AGN variability decreases with BH mass (e.g., Bian & Zhao 2003;O'Neill et al. 2005;Lanzuisi et al. 2014) and UV/optical variations on months to years timescales have a steeper slope with increasing BH (e.g., Caplar et al. 2017;Suberlak et al. 2021). However, none of those studies can provide evidence for a BH mass dependent AGN variability on a timescale of several thousand years. Of course, the overall AGN duty cycle must be much longer, on Gyr timescales, in order to explain the necessary BH mass growth. AGN therefore need to flicker on and off many times during their evolution in concert with the evolution of the host galaxies (e.g., Hickox et al. 2014;Schawinski et al. 2015;Gaspari et al. 2017). Based on our morphological classification we know that the majority of CARS galaxies are isolated disk-dominated galaxies for which we do not expect any significant impact on the flicker cycles by external factors such as mergers or interactions. In addition, the ENLR are clearly associated with a single galaxy and its central BH in the vast majority of cases. This ensures that we cannot systematically mix-up the life-cycles of two independent BHs. This might be an unknown source of confusion for the proximity zones at z ∼ 3 which exceed the sizes of galaxy groups that can host several independent BH being potentially active simultaneously (e.g., Husemann et al. 2019a). Furthermore, the time a BH spends nearly in a quiescent state (t off ) is crucial to understand the full AGN flicker cycle, which we cannot constrain from our observations. If t off is significantly shorter than our measured t AGN , we would already measure a superposition of several cycles, which we cannot rule out based on our observations and is a limitation of our model and data interpretation.

Conclusion
In this paper we describe the full optical IFU data set for CARS with detailed information on the data reduction, QSO subtraction scheme and modeling of the stellar and emission lines. We focus on the characterization of the AGN-ionized region, which is referred to as the NLR or ENLR ranging from subkiloparsec to tens of kiloparsec in size. The main advantage of the type 1 AGN sample of CARS is the possibility of determining primary AGN parameters, in particular BH mass and Eddington ratios, which are difficult to determine for type 2 AGN more commonly targeted for these types of studies. As the CARS sample is a representative subsample of the luminous AGN population at low redshift, we are able to determine the following key results: -We measure ENLR sizes ranging from a few 100 pc to a few tens of kpc for three different ENLR size definitions, which are all not strongly correlated with the AGN luminosity. -Comparing the ENLR sizes with all accessible AGN parameters, redshift and host galaxy properties, we find the strongest correlation between the maximum ENLR size log(R ENLR,max ) and BH mass log(M BH ). -Interpreting the maximum ENLR size as time scale indicator, we recover a scaling relation for the lifetime of a single AGN episode of t AGN ∝ M 0.45±0.07 BH taking into account projection effects and random AGN lifetime sampling applying a generative model to the data. Alternative interpretations are possible if secondary correlation are driving the observed BH mass dependence.
-Our BH mass dependent scaling relation for the AGN episode lifetime is consistent with the ∼3 Myr lifetime estimate of high-mass BH at z ∼ 3 based on the He ii proximity zone sizes ).
The apparent BH mass dependence of the ENLR size is the key result of our study. Nevertheless, we cannot fully rule out that the apparent correlation is a derivative of an even more fundamental hidden relation or through unknown selection effects considering the complexity of the ISM in galaxies, radiative transfer, and AGN variability. However, if the ENLR sizes are correctly interpreted as a proxy for the lifetime of a single AGN episode, it may provide possible explanations for various phenomena. First the famous Eigenvector 1 parameter space (e.g., Boroson & Green 1992;Sulentic et al. 2000) could be interpreted as a time difference in the AGN phase. AGN with more massive BHs would be statistically observed at much later times in their episodic phase than lighter BHs due to the difference in AGN lifetime despite similar luminosity. In this picture, the socalled narrow-line Seyfert 1 (NLS1) galaxies would be younger AGN as previously proposed (e.g., Mathur 2000;Grupe 2004). This may explain their extreme [O iii] blueshifts (e.g., Bian et al. 2005) as well as the relative faintness of [O iii] with respect to the broad Hβ component, because the accretion-disk winds and associated NLR would be more compact leading to higher velocity outflows and higher gas densities prone to collisional de-excitation of forbidden lines. If strong Fe ii emission is predominantly seen in AGN with narrower broad lines (smaller BH mass) as seen in Fig. 11, the Fe ii line strength might be linked to such time evolution effects and decay with time.
The relatively short timescale of t AGN of less than a Myr across a large range of BH masses can provide an explanation for the lack of total star formation suppression observed for moderately luminous AGN samples (e.g., Harrison et al. 2012;Shimizu et al. 2017;Rosario et al. 2018;Scholtz et al. 2020) and for our CARS sample (Smirnova-Pinchukova et al. 2021). For a large range in BH mass, the AGN lifetime would simply be too short to already have an impact on the total SFR of a galaxy with a characteristic time scale of several Myrs. For systems with more massive BHs, the AGN lifetime might be significantly longer and the time integrated energy input is much larger. Therefore, BH mass seems to be an important parameter in addition to the AGN luminosity in controlling the AGN feeding and feedback cycle of AGN, as also suggested by hydrodynamical simulations (e.g. including CCA; Gaspari et al. 2020

Appendix A: Bayesian model for AGN lifetime
In order to statistically infer the average AGN lifetime t AGN from an observed AGN sample, we follow the basic prescription developed by Khrykin et al. (2021). They use the He ii proximity zone sizes of about 20 z ∼ 3 luminous AGN to determine t AGN . Instead of the proximity zones, we use the size of the ENLR as an alternative proxy for the time an AGN has been active during the current episode, the "on-time" (t on ) where t on ∈ [0; t AGN ]. Only t on can be measured and we need to use a predictive model to infer t AGN from the whole sample as t AGN cannot be directly inferred for a given AGN by the nature of the problem. Following Khrykin et al. (2021) we start by assuming a lognormal distribution function for the AGN lifetime as where µ = log(t AGN /Myr]) and σ = σ log t AGN . In order to apply Bayesian inference we need to determine the likelihood function L AGN (R ENLR,i |µ, σ) for each AGN in the sample. We start with the construction of a parameter grid for µ and σ in the range of µ = [−5, 1.5] with a step size of ∆µ = 0.1 and σ = [0.01, 1.0] with ∆σ = 0.05. For each of those grid points we perform the following sequence of steps to determine the likelihood function for each parameter combination: (1) A sample of 1000 AGN lifetimes log(t AGN /[Myr]) is randomly drawn from a log-normal distribution function.
(2) For each of those 1000 values of t AGN we draw 10 values of t on adopting a uniform distribution with p(t on ) = 1 t AGN between 0 and t AGN . This takes into account the random sampling of AGN observations during their overall lifetime.
(3) We convert the drawn t on values into an intrinsic ENLR size R ENLR with the light speed.
(4) Adopting an intrinsic half-opening angle θ ENLR = 60 • for the ENLR ionization cones, the probability density function (PDF) p for the cone inclination i ENLR is given by the ratio of the surface area of a spherical ring dA ring (i) = 2πr 2 sin(i)di to the total area of the cone surface A cone = 2πr 2 (1 − cos(θ)): This PDF implies a mean inclination angle of i ENLR = 39 • , which we apply to convert intrinsic R ENLR sizes to projected R ENLR,proj sizes. We basically assume that the large opening angle of the cones are removing any dependence on object-to-object variations and apply only a global projection.
(5) The continuous PDF of p(R ENLR,proj ) is determined by applying a kernel density estimation (KDE) to the 10 000 drawn values of R ENLR,proj . This way the impact of the random sampling of t on is effeciently marginalized in the combined PDF.
We provide an example of all steps in Fig. A.1 for a single parameter set of µ and σ. The likelihood L AGN (R ENLR |µ, σ) can be obtained by evaluating the PDF of p(R ENLR,proj ) for a specific parameter set of µ and σ. Here, we use a simple linear interpolation of the PDF between our precomputed grid points. Since our data imply that µ is a function of log(M BH ), we assume a linear relation µ = m log(M BH )+b, which simply change the likelihood function to L AGN (R ENLR , log(M BH )|m, b, σ).
The joint likelihood function for our entire sample can then be computed as where N AGN = 32 is the number of AGN observed with MUSE excluding non-AGN and CLAGN sources as well as HE 0021−1810 due to its low data quality. As described in the main text we use the MCMC to sample L joint to infer the posterior distribution function for m, b and σ from our observed data.

Appendix B: Stellar and emission line maps
Article number, page 29 of 39 A&A proofs: manuscript no. CARS_presentation_v3