© The Authors 2023

1. Introduction

Narrow-line Seyfert 1 (NLSy1) galaxies and quasars with relatively narrow permitted lines are generally believed to have relatively high Eddington ratios (e.g., Mathur 2000; Wang & Netzer 2003; Grupe 2004). The NLSy1 galaxies are contextualized as a sub-population within a larger population of sources accreting at relatively high rates (Eddington ratio, λ_Edd > 0.2 for Population A, following Sulentic et al. 2000a). However, uncertainty as to how high these ratios can be remains. This is an important issue from the theoretical point of view, as we have a basic understanding of the accretion when the Eddington ratio is moderate, but the models are unreliable at the extremely high accretion rates. Moreover, strong outflows might, in principle, prevent too high accretion rates from happening in active galactic nuclei (AGN). Therefore, the study of sources considered as highly super Eddington is of extreme importance.

The determination of the Eddington ratio requires both the measurement of the bolometric luminosity of a source and its black hole mass. Currently, the reverberation mapping method is considered the most reliable method for black hole mass measurement (Peterson 1993; Peterson et al. 2004; Cackett et al. 2021). It measures the profiles of the broad emission lines coming from the broad-line region (BLR) and the delay between the variable continuum and the line response, which constrains the velocity and the BLR size. The method is based on the assumption that the BLR is predominantly in Keplerian motion, and the results generally agree with the results based on the black hole mass – bulge velocity relation (see e.g., Ferrarese & Merritt 2000; Gültekin et al. 2009), which supports the underlying assumption. However, some objects with an Hβ full width at half maximum (FWHM) of about 2000 km s⁻¹ show Hβ lags much shorter than the well-known radius-luminosity relation (Bentz et al. 2013) when we compared objects with the same luminosity (Du et al. 2014, 2015, 2016, 2018; Wang et al. 2014a). This implies a surprisingly small black hole mass and, for a measured bolometric luminosity, a very high Eddington ratio.

Such high Eddington objects challenge the theoretical models of the accretion process, and they are important for understanding the rapid growth of the black hole mass at high redshifts. Theoretically, these objects should be modeled with slim accretion disks (Abramowicz et al. 1988). The innermost part of such a disk is geometrically puffed up, and a significant fraction of the radiation is expected to be released in a collimating funnel. In such a picture, the irradiation of the outer parts of the disk, where the BLR is located, is geometrically less efficient (Wang et al. 2014b). Thus, the BLR may be located elsewhere, such as in the biconical flow (Corbett et al. 2000). If so, the BLR we observe would not be in Keplerian motion, and the black hole mass measurement may be highly biased. In addition, in extreme cases of a source seen along the symmetry axis, the lines can also be much narrower due to a purely geometrical factor (e.g., Baldi et al. 2016). A comparison of black hole masses estimated using the reverberation mapping measurements with the masses estimated from the stellar dispersion may not give the final agreement. Due to an evolutionary effect, NLSy1 galaxies are sometimes argued to be outliers from the black hole mass-stellar dispersion relation. Their black holes are too small for their bulge masses, and they accrete vigorously, increasing mass and moving toward the standard relation appropriate for mature AGN (e.g., Mathur 2000; Mathur et al. 2001). Robinson et al. (2011), with the spectropolarimetric measurements for 16 NLSy1, concluded that NLSy1 galaxies represent an extreme realization of one or more physical parameters, such as black hole mass and/or accretion rate, rather than simply being preferentially oriented close to face-on viewing angles. Thus, NLSy1 galaxies are not simply a sub-population that is preferentially viewed close to the AGN symmetry axis, and hence, the orientation of the accretion disk cannot be the main parameter governing the broad-line widths in these sources.

Using spectropolarimetry, it is possible to obtain new and independent insight into the geometry and full velocity field of the line-emitting material as well as into the location of the fully ionized scattering medium. Spectropolarimetry revealed the nature of type 2 AGN as active galaxies with the BLR hidden by the dusty-molecular torus (Antonucci & Miller 1985). Moreover, Smith et al. (2004, 2005) later showed the advantage of using this technique to study type 1 AGN as well. The wavelength-dependent polarization angle also indicates whether the scattering takes place in the polar or the equatorial region and measures the viewing angle of the system. In low viewing angle sources, emission lines in polarized light are much broader than in the unpolarized spectrum. Even simply comparing the kinematic line width in total and polarized light can verify the statement about the actual value of the black hole mass and the Eddington ratio of a source (Baldi et al. 2016). A new independent method of black hole mass measurement has been introduced, and it is based on the change of the polarization angle across the broad emission line. This change depends on the rotational velocity in the BLR, and it allows for direct measurement of the black hole mass (Afanasiev & Popović 2015). The method was proposed by Afanasiev et al. (2014, 2015) and applied to mostly Hα-line based observations (e.g., Afanasiev & Popović 2015; Afanasiev et al. 2015, 2019) using the spectropolarimetry data from the Very Large Telescope (VLT)/FORS2. Meanwhile, Jiang et al. (2021) has shown that for Fairall 9, the mass measurements from Hα and Hβ lines are consistent.

The first spectropolarimetric black hole mass measurement with the MgII line using the rotation of the polarization angle was recently performed by Savić et al. (2020). Their result is in good agreement with other black hole mass estimation techniques (see Afanasiev & Popović 2015; Afanasiev et al. 2019).

We present the spectropolarimetric measurements and spectral decomposition of the Hα spectral region of three NLSy1 galaxies, Mrk 1044, SDSS J080101.41+184840.7, and IRAS 04416+1215¹, performed with the VLT/FORS2 instrument. We also show the spectral decomposition of the archival spectra for these three sources in Hβ and Hα regions. All three sources are part of the super-Eddington accreting massive black holes (SEAMBH) Project sample (Wang et al. 2014a), selected for reverberation monitoring as candidates for super-Eddington accreting sources. In the following paragraphs, we briefly describe the archival information of the three chosen sources.

The source Mrk 1044 (z = 0.016451 ± 0.000037; RA 02:30:05.5, Dec −08:59:53 from NED²) is one of the nearest and brightest Seyfert galaxies (V mag = 14.5, from NED; Véron-Cetty & Véron 2006). The line width of the Balmer lines was already measured by Rafanelli & Schulz (1983). The authors reported an FWHM of 3600 km s⁻¹ for the broad components (BCs) of both the Hα and Hβ lines. They stressed the comparable intensity of the narrow components (NCs), with an FWHM width of 1000 km s⁻¹, and the very low value of the [OIII]/Hβ ratio (below 0.6). The source was later reclassified as an NLSy1 galaxy (Osterbrock & Pogge 1985). This source was also observed in X-rays, and the viewing angle obtained by Mallick et al. (2018) is i = $47.2_{-2.5}^{+1.0}$, from the joint fitting of Swift, XMM-Newton, and NuSTAR X-ray spectra.

The source SDSS J080101.41+184840.7 (z = 0.13954 ± 0.00001; RA 08:01:01.41, Dec +18:48:40.78 from NED) with V mag = 16.88, from Véron-Cetty & Véron (2010), is the second object in our sample. Liu et al. (2021) modeled XMM-Newton spectrum of SDSS J080101.41+184840.7 and obtained a steep photon index (2.33 ± 0.06) using a power-law model modified by Galactic absorption. However, in the model used by Liu et al. (2021), the viewing angle was not one of the model parameters.

Tortosa et al. (2022) recently reported a detailed analysis of the broadband observations of IRAS 04416+1215 (r mag = 16.24, z = 0.089). This source shows a narrow Hβ line (FWHM = 1670 km s⁻¹; Moran et al. 1996) and a very broad [OIII] lines (FWHM = 1150 km s⁻¹; Véron-Cetty et al. 2001). In their work, Tortosa et al. (2022) showed that the best-fitting model in the X-ray band is composed of a soft excess, three ionized outflows, neutral absorption, and a reflection component. Prominent soft X-ray excess supports the claim for a super-Eddington accretion rate in this source. Through a dedicated monitoring campaign using the Lijiang 2.4 m telescope under the SEAMBH Project, Du et al. (2016) reported the measurements of the Hβ FWHMs for the three sources (see Table 5 in their paper): Mrk 1044 (1178 ± 22 km s⁻¹), SDSS J080101.41+184840.7 (1930 ± 18 km s⁻¹), and IRAS 04416+1215 (1522 ± 44 km s⁻¹).

This paper is organized as follows: Details on the performed observations and data reduction are reported in Sect. 2. In Sect. 3, we present the analysis of spectropolarimetric data obtained using the VLT and the archival spectroscopic data for each source. In Sect. 4, we show the spectral decomposition for natural and polarized light for our objects and spectral modeling. Next, we perform polarization radiative transfer modeling using STOKES and compare it with our observed results in Sect. 5. We discuss the implications of these results, outline the physical parameters responsible for the trends observed in our sample using principal component analysis, and summarize our findings from this study in Sect. 6. In this work, we use ΛCDM cosmology (Ω_Λ = 0.7, Ω_M = 0.3, H₀ = 70 km s⁻¹ Mpc⁻¹).

2. Observations and data reduction

Our spectropolarimetric observations were performed with the FORS2 instrument mounted on the UT1 telescope of the 8.2 m ESO VLT. The observations were taken using the GRISM-300I in combination with the blocking filter OG590. The spectral range of our data is 6000−10 000 Å. The 0.7-arcsec wide slit was oriented along the parallactic angle, and the multi-object spectropolarimetry observations were performed with a 2048 × 2048 pixel CCD with a spatial resolution of 0.126 arcsec pixel⁻¹.

The spectropolarimetric observation of Mrk 1044 was accomplished on 2016 October 27. The source SDSSJ080101.41+184840.7 was observed on 2016 December 26, 2016 December 30, and 2017 January 2. Finally, the IRAS 04416+1215 was observed on 2016 November 6. For all observational runs, we used the Patat & Romaniello (2006) prescription of observing with four quarter-wave plate angles (0, 22.5, 45, 67.5 deg.). For Mrk 1044, two exposures were taken at each angle (4 × 2 spectra), with each exposure lasting 206 s. For IRAS 04416+1215, three exposures lasting 455 s were taken at each angle. For SDSS J080101.41+184840.7, three exposures of 250 s each were taken on 2017 January 2, and three exposures of 280 s each were taken during the observations in December 2016.

Together with the sources, two standard stars were observed: one star that is highly polarized and one star that is unpolarized (Ve 6−23 and HD 62499, respectively). Observations were performed in the service mode (program ID: 098.B-0426(B), PI: B. Czerny), and observing conditions were good, with an atmospheric seeing ∼0.8. Since only one unpolarized star was observed and it was not located close to any of our sources in the sky, we decided not to subtract the effect of the interstellar medium from our measurements. This is certainly justified for Mrk 1044 and SDSS J080101.41+184840.7, although not quite for IRAS 04416+1215 (see Appendix A).

To reduce data, we used ESO REFLEX³ (Hook et al. 2008). The workflow combines the bias frames into a master bias that is then subtracted from the science image. The same procedure is performed with the lamp flats, and the science image is flat-fielded. The workflow then removes the cosmic ray events and performs the wavelength calibration using the standard He-Ar arc lamps. As final products, among others, we obtained the extracted 1D spectra with the Stokes parameters (U and Q), total linear polarization (L), and flux in ADU/s (I) as a function of wavelength with the associated uncertainties.

A comparison between the natural light spectrum and the one from the SDSS catalog showed that the AGN continuum form was different in both spectra. The same effect was observed in the reduced standard stars obtained with ESO REFLEX, where an expected blackbody continuum was not recovered. This difference suggests that one correction was not applied by the pipeline. To correct the continuum form in the natural light spectrum, we used the IRAF routines “standard” and “sensfunc”, and the standard star HD 64299. After the correction, a good agreement was obtained between the natural light spectrum and the one from the SDSS catalog.

3. Data analysis

3.1. Standard star

In Table 1, we present the measurements for the standard stars HD 64299 and Ve 6−23, used in our work. Cikota et al. (2017) provide measurements for Ve 6−23⁴ between 7.20% and 7.03% for the I band. In our case, we obtained 7.12%, which is comparable to the average value of Ve 6−23 polarization (7.13%) from Cikota et al. (2017). Sosa et al. (2019) obtained a 6.23 ± 0.03% polarization level for this star, but observations were performed with a different instrument⁵. The polarization angle for Ve 6−23 is 170, which is consistent with the angle Cikota et al. (2017) calculated (171.95 ± 0.01). In the case of HD 64299, Sosa et al. (2019) performed the measurements for this star in the I band and obtained P = 0.16 ± 0.01% (Table 3 therein). Our measurements for this star, P = 0.018%, differ from Sosa et al. (2019), which may be connected with different CasPol instrumental polarization in the I band. Sosa et al. (2019) report the inconsistent measurements in the I band between CasPol literature. However, results obtained by CasPol in B and V bands, seem to be consistent with the ones from the literature, including FORS2, (for details see Sosa et al. 2019). The instrumental polarization is small, and the standard measurements were reproduced within errors, in particular the angle that depends on the half-wave plate (HWP) achromaticity.

3.2. Spectropolarimetry of NLSy1

In this work, we present three types of polarization measurements: the mean polarization of each source, the polarization of the continuum, and the polarization of the line corrected for continuum polarization. By a mean polarization (P_mean), we refer to the polarization calculated from the mean Q and U of the whole wavelength range of the spectra. As a continuum polarization (P_cont), we consider the polarization calculated from the mean Q and U at two windows surrounding Hα from its blue and red sides. Each window has a width of 100 Å approximately, and they are separated from the center of Hα by at least 300 Å. The continuum polarization angle (χ_cont) was also calculated from the mean Q and U of the same windows. For the line polarization measurement (P_line), we considered the polarization of the region with a width of approximately 50 Å surrounding the center of the emission line.

In addition, we corrected for the contamination of the continuum underneath the emission line (from the Stokes parameters, which we mark as ${\tilde{U}}_{\mathrm{line}}$ and ${\tilde{Q}}_{\mathrm{line}}$). We estimated the Stokes parameters, Q and U (normalized to the intensity, I), using the two windows (width = 100 Å) surrounding Hα as described above. Then, we calculated the wavelength-averaged vector of Q_cont ⋅ I_cont and U_cont ⋅ I_cont for the continuum in these two windows. Finally, we calculated the mean value of Q_cont ⋅ I_cont and U_cont ⋅ I_cont for the two windows and subtracted this mean contribution from each point of the Hα line vector, as shown in formulas below:

$$\begin{array}{cc}& Q_{\mathrm{line}}=({\tilde{Q}}_{\mathrm{line}}·I_{\mathrm{line}}-Q_{\mathrm{cont}}·I_{\mathrm{cont}})/I_{\mathrm{line}}\hfill \end{array}$$(1)

$$\begin{array}{cc}& U_{\mathrm{line}}=({\tilde{U}}_{\mathrm{line}}·I_{\mathrm{line}}-U_{\mathrm{cont}}·I_{\mathrm{cont}})/I_{\mathrm{line}}.\hfill \end{array}$$(2)

The further procedure of calculating the fraction of linearly polarized radiation is the same as for the P_mean and P_cont.

3.3. Archival spectra of NLSy1

We decided to explore the physical parameters of our sample in the optical band starting from archival reduced data. For Mrk 1044, we used spectrum from Jones et al. (2009), and for SDSS J080101.41+184840.7 and IRAS 04416+1215, we used data from SDSS (Blanton et al. 2017). For these spectra, we performed spectral decomposition of the Hβ and Hα regions. Then, we used values obtained for the Hα region from archival spectra as initial values for the decomposition of non-polarized spectra of Hα from FORS2/PMOS.

To analyze all spectra we used the specfit task from IRAF (Kriss 1994), including various components in two ranges: 4400−5200 Å for the Hβ region and 6200−6800 Å for the Hα region. In the Hβ range, we followed the method presented in Negrete et al. (2018). In both regions, we modeled the accretion disk emission and the starlight contribution to the spectrum assuming a power-law shape without performing a more complicated procedure for the starlight component, since the ranges are relatively narrow. To model the continuum, we used the continuum windows around 4430, 4760, and 5100 Å in case of Hβ (e.g. Francis et al. 1991), and in case of Hα, we used 5700−6000 Å (e.g. Kuraszkiewicz et al. 2002). In both regions, we considered the Fe II template from Marziani et al. (2009).

The [OIII]λλ4959,5007 doublet lines in the Hβ region were modeled with two Gaussians that represent narrow and semi-broad components. We assumed the theoretical ratio of 1:3 between the strengths of the components (Dimitrijević et al. 2007). The [N II]λλ6548,6584 and [S II]λλ6716,6731 doublet lines in the Hα region were modeled with Gaussian profiles, with the assumption of the theoretical ratio 1:3 between their strengths and the FWHM of lines fixed to the [OIII] narrow component (NC). For each Balmer line, we decided to use the Lorentzian shape for the broad component (BC) since it is considered to be more suitable for the NLSy1 than the Gaussian (e.g. Laor et al. 1997). We modeled the Balmer emission lines taking into account two components: one Lorentzian component and one narrow Gaussian component (to reproduce the peak of the line). Moreover, we kept the FWHM of the narrow Gaussian component fixed to the FWHM of the NCs for all forbidden lines and relaxed it only for the last iteration of the fitting. To estimate errors, we used the maximum and minimum continuum levels, and for each case, we performed the fit of the spectrum. We assumed that the distribution of errors follows the triangular distribution (d’Agostini 2003). For each line measurement, we calculated the variance using the following formula for the triangular distribution:

$$\begin{array}{c}\hfill {\sigma }^2(X)=\frac{{\mathrm{\Delta }}^2{x}_{+}+{\mathrm{\Delta }}^2{x}_{-}+\mathrm{\Delta }{x}_{+}+\mathrm{\Delta }{x}_{-}}{18}\end{array}$$(3)

where Δx₊ and Δx₋ are the differences between measurements of the maximum and best continuum and between the best and minimum continuum, respectively. This formulation can be explained as the linear decrease in either side of the maximum of the distribution (which is the best fit) to the values obtained for maximum and minimum contributions of the continuum. We find this analytical method easy and sufficiently precise. We propagated uncertainties using standard formulas of error propagation for the reported values in this work.

4. Results

4.1. NLSy1 spectra decomposition in natural light

We fit each spectrum individually, and for each, we followed the same steps to keep the same accuracy. As a first step, we decomposed the Hβ and Hα emission line profiles in natural light from archival spectra. Then, we fit Hα in natural light obtained from our FORS2/VLT observations. To remain consistent, we kept the same values for the Hα components obtained from archival spectra as input and then fit the Hα in natural light from VLT. We present our measurements of FWHM of Balmer lines’ components in Table 2. Each line was modeled with three components: broad (Lorentzian), narrow (Gaussian), and blue (Gaussian). In all cases, the FWHM of the Hα BC is slightly narrower than the Hβ one, whereas the blue components are broader. We present the spectral decompositions together with those done for archival spectra in Fig. 1. We did not see any significant differences between Hα fits from archival spectra and the VLT spectra. For further analysis, we used the archival spectra of Hβ (since we did not perform observations for this line) and the VLT Hα spectra. The Balmer lines are slightly asymmetric, which may indicate the presence of an outflowing component. The presence of a blue, outflowing component in Balmer lines is visible for all three objects. In the case of IRAS 04416+1215, it is visible for both lines and Mrk 1044, it is visible for Hβ. However, the Hβ blue component in Mrk 1044 is almost negligible. For SDSS J080101.41+184840.7, the blue components are the weakest but are still needed to obtain the lowest residual values. All objects from our sample show prominent Fe II emission around the Hβ region, but the [OIII] doublet and Hβ NC differ considerably between objects. The Fe II emission around the Hα region is not significant, and the lower contamination in this region is a general property observed in other sources (e.g., in I Zw 1; Véron-Cetty et al. 2004 Fig. 7 therein) and in CLOUDY photoionization models (Sarkar et al. 2021).

Fig. 1. Natural light spectra for Mrk 1044 (top row), SDSS J080101.41+184840.7 (middle row), and IRAS 04416+1215 (bottom row). Top row: normalized 6dF spectrum of Mrk 1044 from Jones et al. (2009) with the fitted Hβ (left) and Hα (middle) part. The right plot shows the normalized unpolarized spectrum obtained with FORS2/VLT. The middle and lower rows show as follows for SDSS J080101.41+184840.7 and IRAS 04416+1215, respectively: fitted Hβ and Hα part from SDSS normalized archival spectrum and normalized unpolarized spectrum obtained with FORS2/VLT. Data are marked in gray, and the model is marked in black. The fit components are given in the following: for Hβ (left panel), the Hβ BC and [O III] semi-broad component doublet are red. The NCs of the Hβ and [OIII] doublet are black. The Hβ blue component is marked in blue. The power-law component is marked in blue, and the Fe II pseudo continuum is green. For Hα (middle and right panel), the Hα NC and [NII] and [SII] doublets are marked in black. The Hα BC is marked as red. The Hα blue component is marked in blue, the Fe II pseudo continuum is marked in green, and the power-law component is marked in blue. The lower panels of each plot correspond to the residuals, in radial velocity units km s−1. All spectra are normalized to Hα intensity.

In the case of Mrk 1044, the Hβ and Hα NC intensities are similar to BCs of those lines. We reproduced relatively well the profile of the [OIII] doublet using just one component, NC. The noticeable asymmetry in [OIII]5007 Å was nicely reproduced by fitting the Fe II template. In the case of SDSS J080101.41+184840.7, the Hβ and Hα NC of those lines is around one-third of the BC. The intensity of the semi-broad component of [OIII] is comparable to the NC of this line. In IRAS 04416+1215, the NC of Hβ and Hα play the least significant role in the modeling of the full profile. The most prominent, in comparison to the other two sources, asymmetry in the red side of Hα is caused by the [NII] doublet, whose intensity is the highest for this source. The semi-broad component of [OIII] is dominant, and it is the only source from our sample with that feature.

4.2. Polarization measurements

We present the spectropolarimetric measurements for binned data in Table 3. For each source, P_mean is less than 1%. We calculated P_mean (Col. 2) from the mean of the Stokes parameters Q and U of the whole wavelength range of the spectra using the following formula:

$$\begin{array}{c}\hfill P_{\mathrm{mean}}=\sqrt{Q_{\mathrm{mean}}^2+U_{\mathrm{mean}}^2}.\end{array}$$(4)

The polarization of the continuum region (Col. 4) and line region (Col. 6) is calculated similarly, including different wavelength ranges of the spectra. Additionally, the polarization of the lines includes the “continuum polarization correction”, which we explained in Sect. 3.2. The “continuum polarization correction” does not change the physical interpretation of the results. For the SDSS J080101.41+184840.7, we obtained a P_line of 0.20% with the correction, and a P_line of 0.42% without this correction. Both values are lower than the P_mean and P_cont for this source, but by including this correction, the fraction of linearly polarized radiation was reduced.

The polarization position angle (which is the angle of maximum polarization) was calculated using the following formulas (Eq. (8) in Bagnulo et al. 2009):

$$\begin{array}{c}\hfill \begin{array}{cc}& \chi =\frac{1}{2}{tan}^{-1}\left(\frac{U}{Q}\right)(\mathrm{if}Q>0\mathrm{and}U\ge 0);\hfill \\ \hfill & \chi =\frac{1}{2}{tan}^{-1}\left(\frac{U}{Q}\right)+{180}^{°}(\mathrm{if}Q>0\mathrm{and}U<0);\hfill \\ \hfill & \chi =\frac{1}{2}{tan}^{-1}\left(\frac{U}{Q}\right)+{90}^{°}(\mathrm{if}Q<0).\hfill \end{array}\end{array}$$(5)

We used the mean Q and U from the whole region (Col. 3), the continuum region (Col. 5), and the line region (Col. 7). We show the intensity of the Hα, the polarization percentage, the Stokes parameters, and the polarization angle for each object in Fig. 2. We show the windows for continuum and the line polarization for objects from our sample in red. We used the same windows for each object. Our results are not sensitive to the exact width of the continuum. We tested a few cases between 100 Å and 400 Å and did not notice a significant change in the obtained values.

Fig. 2. Data of Mrk 1044, SDSS J080101.41+184840.7, and IRAS 04416+1215. Plots, from top to bottom: total flux, Stokes parameters U and Q, flux of U and Q, and the polarization angle. We mark in red the windows in which we calculated polarization. For continuum polarization, we considered two surrounding Hα from its blue and red side, each window has a width of ≈100 Å, and they are at a distance of at least 300 Å from the center of Hα. For the polarization of the line, we considered the red window mark on the line. Each central window has width of ≈50 Å. The vertical blue dotted line marks the central wavelength for Hα. The data were rebinned with a 2 Å binning.

Overall, the level of polarization in our three sources is very low, less than 1%, which tends to be positively biased (see Simmons & Stewart 1985). The measurement of the change of the polarization angle across the broad emission line is hard to measure for sources with low polarization. Therefore, we concentrated on the estimates of the viewing angle to improve the black hole mass measurement.

In addition, we noticed that the polarized flux profile has a peak-like appearance, with a width not much larger than the one measured in natural light, at least for Mrk 1044 and SJ080101.41+184840.7 (the case of IRAS 04416+1215 is less clear due to its pattern in degree of polarization being dominated by noise). This can be an indication of a relatively high viewing angle (as discussed below in Sect. 4.4) or as evidence of a polar scatterer (Sect. 5). The evidence of a polar scatterer is the strongest in the case of Mrk 1044. In the continuum not corrected for polarization, the degree of polarization of the line (0.33%) is higher than the one of the adjacent continuum (0.22%), an increase by a factor of 1.45. With the correction, the continuum and line polarization degrees are on a similar level. Moreover, χ changes by approximately 90°, which is a typical signature of polar scattering as well.

4.3. NLSy1 spectra decomposition in polarized light

In Fig. 3, we show the spectral decomposition of VLT Hα emission line profiles in the polarized light. The expected equatorial placed scatterer was modeled as a single BC broader than the Hα in natural light. The case in which just a single component is enough to model the line profile is visible in the bottom panel of the figure (IRAS 04416+1215). The cases of Mrk 1044 and SDSS J080101.41+184840.7 are different. Modeling the full line profile using just a BC was not possible. Thus, we decided to add the NCs marked in black. The NCs are most probably associated with a non-equatorial scatterer, and we investigated Mrk 1044 and SDSS J080101.41+184840.7 using the STOKES radiative transfer code in this regard (see Sect. 5.2). The residuals for Mrk 1044 and IRAS 04416+1215 are higher than for SDSS J080101.41+184840.7. Due to low S/N in our data, we were not able to obtain better quality fits; however, there is no specific pattern in the residuals, which supports the reliability of our performed fits.

Fig. 3. Fitted Hα part from the polarized spectrum obtained with FORS/VLT. The data are marked in gray, and the model is marked in black. For each fit, we considered two components of Hα: a BC (red) and an NC (black). The blue line corresponds to the power-law continuum. The lower panels of each plot correspond to the residuals, in radial velocity in units of km s−1. With an arrow, we mark the presence of an “absorption” in the fit of SDSS J080101. The “absorption” is meant to account for the low S/N of the data and to ease the fit of a symmetric Gaussian constrained by the blue side of the polarized flux profile.

4.4. Viewing angles

To estimate the viewing angles, we used Eq. (8) from Collin et al. (2006):

$$\begin{array}{c}\hfill \mathrm{\Delta }{V}_{\mathrm{obs}}\approx V_{\mathrm{Kep}}{[{(H/R)}^2+{sin}^{2}i]}^{1/2}.\end{array}$$(6)

We assumed that the FWHM of the line in polarized light measures the true Keplerian velocity at the distance R (i.e., V_Kep), while the observed value of ΔV_obs is measured as the FWHM of the line in non-polarized natural light, and i is the viewing angle of the system toward us. With H/R, we denote the aspect ratio of the disk at any radius. From the observer’s line of sight, the “full” Keplerian velocities from the rotating disk are not visible, as the equatorial scattering region and the rotating disk are on the same plane (see Fig. 9 Smith et al. 2005). However, the equatorial scattering region is well positioned for observations of unprojected velocities from the rotating disk.

The ratio of the thickness of the BLR region to its radius is indicated by H/R, and it represents the random motion of BLR clouds. The physical justification for an H/R higher than 0.1 comes from Collin et al. (2006). As the authors underlined, the illumination required for the BLR implies a large opening angle of the BLR. The most typical H/R ratio found for the BLR is 1:3 (e.g., Shaw et al. 2012), and we adopt this value in further considerations. The results are reported in Table 4.

4.5. Estimation of M_BH

The reverberation mapping technique allowed us to determine the black hole mass via a virial relationship:

$$\begin{array}{c}\hfill M_{\mathrm{BH}}=f\frac{R_{\mathrm{in}}FWH{M}^2}{G},\end{array}$$(7)

where G is the gravitational constant, and f is the (dimensionless) virial factor that depends on the structure, kinematic, and viewing angle of the BLR. The emissivity-averaged radius of the BLR (also known by the size of the BLR) is denoted as R_in, and FWHM is the full width at half maximum of the emission line profile (for details, see Panda et al. 2019b). The radius of the BLR is determined from the observed time-lag between the continuum and the broad Hβ emission line, and the FWHM is determined from the spectral modeling of the Hβ. The virial factor is more puzzling since it depends on more than one property, and mostly, it is assumed as a fixed factor. With the spectropolarimetry technique, it is possible to estimate the black hole mass using the FWHM of the polarized line in the formula.

The analytical formula for the virial factor in the case of the disk-like structure of thickness H (assumed ≪R_in) is given by:

$$\begin{array}{c}\hfill f={[4(sin{(i)}^2+{(H/R_{\mathrm{in}})}^2)]}^{-1}\end{array}$$(8)

If we use a polarized line to estimate the black hole mass, we assume that we are observing the line-emitting medium locating the observer on the disk-plane. Thus, sin(i) = 1 (i.e., i = 90°).

Including the f factor with polarized line assumptions, we used the following equation to compute the black hole mass for our sources,

$$\begin{array}{c}\hfill M_{\mathrm{BH}}=\frac{R_{\mathrm{in}}FWH{M}^2}{4G}·\end{array}$$(9)

Here, FWHM is the full width at half maximum of the BC of the polarized Hα line obtained in this work.

In Table 5, we present a comparison between the estimation of the black hole mass for our sample from Du et al. (2016) using the reverberation mapping technique (Col. 2) and using Eq. (7) (Col. 3). We note that in Eq. (9) we still used values obtained with reverberation mapping (i.e., the inner radius of the BLR). In this work, the inner radius of the BLR is assumed to be equivalent to the emissivity-weighted radius of the BLR, which we estimated using reverberation mapping studies.

After using our method, the mass estimations for IRAS 04416+1215 and SDSS J080101 are both within the (one sigma) error bars in comparison to the reverberation mapping method (Du et al. 2016). The difference in the mass estimations is only significant for Mrk 1044, as we obtained a lower value from our measurement. However, the Mrk 1044 spectra quality is the best, and we did not have the same problem estimating the viewing angle as we had for IRAS 04416+1215.

4.6. Determination of R_FeII

The strength of the optical Fe II emission, that is, the ratio of the Fe II optical flux measured between 4434−4684 Å range to the broad Hβ flux (also known as R_FeII), seems to be an indirect tracer of the Eddington ratio (e.g., Shen & Ho 2014; Marziani et al. 2018; Panda et al. 2018, 2019a,c; Martínez-Aldama et al. 2021). Since we had already modeled the optical spectra around the Hβ region, we decided to also explore our sample in the context of R_FeII. We show equivalent widths of Fe II, Hβ, and R_FeII in Table 6. The source IRAS 04416+1215 has the highest R_FeII value in our sample and may be called an extreme case, as R_FeII > 1.3, according to Śniegowska et al. (2018), and R_FeII > 1, according to Sulentic & Marziani (2015). From the latter criterion, Mrk 1044 is a borderline high-accreting source. We expected relatively high values of R_FeII (≈1) for all sources in our sample, since a high Eddington ratio was one of our criteria to create the sample. The two other sources, Mrk 1044 and SDSS J080101 are not extreme Fe II emitters, but they do have significantly above average R_FeII values of 0.64. The median value is 0.38 in the Shen et al. (2011) catalog if only the measurements with errors below 20% are considered (Śniegowska et al. 2018).

5. STOKES modeling

Our analysis of the polarized and unpolarized spectra shows that our sources, selected as highly accretion objects, indeed have intrinsically narrow lines and small black hole masses. Therefore, high accretion rate values suggested by Wang et al. (2014a) are not an artifact of the special orientation of these sources. Such highly accreting sources are expected to have very intense outflows. For example, Mrk 1044 shows clear signatures of ultra-fast outflows (UFOs; Krongold et al. 2021). Such strong, high density, highly ionized wind with column density of ∼10²³ cm⁻² and a velocity of 0.15 light speed suggests the presence of a copious amount of material surrounding the black hole, which could lead to efficient polarization while our observations imply very low polarization in these sources. To address the potential discrepancy between strong outflow and low polarization, we performed simulations of the polarization properties of our sources. With those simulations, we intended to independently constrain the amount of scattering material surrounding the AGN.

We used the radiative transfer code STOKES (Goosmann & Gaskell 2007; Marin 2018; Savić et al. 2018) to understand the polarization in the Hα emitting BLR of these three objects. The STOKES code is a 3D radiative transfer code that incorporates the Monte Carlo approach in order to trace every single photon emitted from a source, of which a considerable number are scattered by the intervening media. The code includes physical processes such as the electron or dust scattering until the photons get absorbed or eventually escape and reach the distant observer. We used the latest version (v1.2) of the code⁶.

5.1. Parameters of the model

In our modeling, we assumed the continuum source to be point-like located at the center, and emitting isotropic unpolarized radiation. The generated flux thus generated takes the form of a power-law spectrum F_cont ∝ ν^−α with α = 2. This assumption of the value for the spectral index is justified based on prior works (Savić et al. 2018; Jiang et al. 2021) where the focus was to reproduce the polarization properties around a specific emission line (in our case Hα) and the continuum around that line. The continuum source is surrounded by the line-emitting BLR of cylindrical geometry that is specified by three spatial parameters, namely (a) the distance between the continuum source to the center of the BLR (R_mid); (b) the height of extension of the BLR above (or below) the mid-plane joining the source and the BLR (a); and (c) the half-width of the BLR (b). These three spatial terms can be estimated using the information of the BLR’s inner (R_in) and outer (R_out) radii, which we extracted from previous studies. Thus, for the inner radius of the BLR, we used the reverberation mapped estimates compiled by Du et al. (2016) for the three objects considered in this work, while we estimated the outer radius of the BLR using the analytical relation by Netzer & Laor (1993) for the dust sublimation radius, which has the following form:

$$\begin{array}{c}\hfill R_{\mathrm{out}}^{\mathrm{BLR}}=0.2L_{\mathrm{bol},46}^{0.5}\end{array}$$(10)

where L_bol, 46 is the bolometric luminosity in units of 10⁴⁶ erg s⁻¹. The L_bol for our objects were estimated by scaling the observed monochromatic luminosity at 5100 Å compiled in Du et al. (2016) by a luminosity-dependent bolometric correction factor (Netzer 2019). The assumption of a slightly smaller inner radius than the reverberation mapped estimate (or the emissivity-weighted radius) retrieves nearly identical results for the Stokes parameters, suggesting that, for these sources, assuming the emissivity-weighted radius as a proxy for the inner radius is justified. This is due to the fact that the vertical extension of the BLR scales with the radial location of the onset of the BLR which is estimated based on a fixed half-opening angle of 15°. We tested two scenarios where the R_in was assumed to be smaller by a factor of two and three with respect to the emissivity-weighted radius. The R_mid (which is computed as an average of the R_in and R_out) changes from ∼26 ld to 22.5 ld and 23.37 ld, for the two alternate cases, respectively. This corresponds to an increase by ∼50 km s⁻¹ in the average velocity estimated at R_mid. The results from the STOKES modeling show very little change in the total and polarized spectra, while we noticed marginal variations in the polarization angle due to the change in the location of the inner radius. The remaining parameter to complete the information for the BLR is the velocity distribution for the Hα line. This requires the value of the average velocity (v_avg) along the x − y direction (or ϕ direction; here we assumed the z-axis to be aligned with the BH spin axis). We estimated the average velocity by assuming a Keplerian velocity distribution at a given radius (here, the radius considered is the mean value between the BLR’s inner and outer radii) from the source with a predetermined BH mass, that is, $v_{\mathrm{avg}}=\sqrt{GM/R}$, where G is the gravitational constant.

For the scattering region that represents the torus, we assumed a geometry of a flared disk characterized by inner and outer radii as well as a half-opening angle. The light becomes scattered predominantly due to free electrons (i.e., Thomson scattering), where the content of this medium is regulated by the number density parameter. In addition to these assumptions, the velocity distribution was again set similarly to the scenario for the BLR but assuming the appropriate inner and outer radii consistent for the torus. The inner radius for the scattering region was set based on the R–L relation from the infrared reverberation mapping (Kishimoto et al. 2007; Koshida et al. 2014). We adopted the prescription of Savić et al. (2018) for the outer radius of this region, that is, the location of the outer radius of the torus was set to be at the location where the outer radius of the BLR subtends a half angle of 25°. The combined information from the inner and outer radius of the scattering region plus the electron number density allowed us to estimate the optical depth of this region. We obtained a grid of electron number density solutions to arrive at a preferred solution that fits best the spectral information gathered from the observed data. Table 7 lists the various input parameters used in our modeling for the three objects.

In Fig. 4, we show the schematic illustrations of the model setup for our equatorial scatterer, which includes the continuum source, the BLR, and the equatorial scattering region. We utilized the equatorial scatterer option, as in Jiang et al. (2021), but we also considered the option of the polar scatterer (see e.g., Smith et al. 2005). To model the polar scatterers, we assumed the scatterers to be located at the same distance from the continuum source and analogous in size (R_out − R_in) to the equatorial scattering region. We also kept the structure of the polar scatterers identical to their equatorial counterparts. We modulated the net density of the polar scatterer until we recovered a good fit for the polarized spectrum and the polarization fraction. This setup is illustrated in Fig. 5.

Fig. 4. Schematic view of the models with the equatorial scatterer. The continuum source is marked in gray, the BLR is marked in blue, and the equatorial scattering region is marked in red. We marked the half-opening angle of the scattering structure with a dashed curve between the gray and black dashed lines.

Fig. 5. Schematic view of the models with two scatterers: polar and equatorial. The various components are colored identically to Fig. 4. The polar scattering region is marked in purple.

5.2. Comparison with the observed polarization properties

To reconstruct the observed spectrum and the corresponding polarized emission of the Hα emission line and the continuum around the line, we generated a simulated spectrum using STOKES for each of the three objects. Specifically, we considered the spectrum in natural light (F_λ), the polarized spectrum (p × F_λ), the polarization fraction (p) and the polarization angle (χ). We ran a simulation of these four entities using the parameterization as listed in Table 7, and we compared our results to those from observations.

In addition to the above setup, we also included a polar scattering region that helps recover the peak distribution in the polarized spectra and fractional polarization for Mrk 1044 and SDSS J080101.41+184840.7. These two objects are the ones where we recovered a viewing angle value that is slightly on the higher side (i.e., ∼54° for Mrk 1044 and ∼31° for the SDSS J080101.41+184840.7). However, for IRAS 04416+1215, which is the source for which we recovered the smallest viewing angle (i.e., ∼4°), we could reconstruct the polarized emission using the model with only the equatorial scattering region.

The polarization fraction (p%) for the three sources (Mrk 1044, SDSS J080101.41+184840.7, and IRAS 04416+1215) including the polar scattering region are 0.40, 0.36, and 0.3, respectively. For IRAS 04416+1215, we found the p% for the case with only the equatorial scatterer to be similar (i.e., 0.28). We note that the parameterization for the polar scattering region is unrestrained, as we lacked observational pieces of evidence to break the degeneracy between the optical depths, the density, the location, and the extension of these media. In addition, the half-opening angle for the polar scatterers was assumed to be identical to their counterparts in the equatorial region. With the assumption that the density of the medium remains rather unaffected there is a coupling between the opening angle and the optical depth of the scattering medium, i.e., an increase/decrease in the opening angle will lead to a decrease/increase in the optical depth along a given line of sight. The parameters used to model this polar scattering region for the two objects are tabulated in Table 8.

We illustrate the performance of the recovery of the spectrum in natural and polarized light for our three objects, Mrk 1044, SDSS J080101.41+184840.7, and IRAS 04416+1215, in Figs. 6–8, respectively. The polarization fractions on these plots (third panel from the top in Figs. 6–8) are notably higher than the mean values reported in Table 3. We note that the mean values of the computed polarization fractions depend on the wavelength range employed and on the order of the executed steps. For example, for Mrk 1044, we obtained 0,67% when computing the mean value for P using P value for each wavelength point for the full wavelength range (6101−10 295 Å). Whereas, when we computed P value using the means of the respective Q and U values (as explained in Sect. 4.2), we obtained lower values (as presented in Table 3). We note that the difference between averaging P and averaging (Q, U) before computing P is due to the noise that introduces a bias on P. In the absence of noise, both the averages would be identical.

Fig. 6. STOKES modeling for Mrk 1044 with the equatorial plus polar scattering region and comparison with observational estimates. The case shown has a viewing angle of 54°. From top to bottom: spectrum in natural light (Fλ), polarized spectrum (p × Fλ), polarization fraction (p), and polarization angle (χ). The vertical red dashed line marks the central wavelength for Hα.

Fig. 7. STOKES modeling for SDSS J080101.41+184840.7 with the equatorial plus polar scattering region and comparison with observational estimates. The case shown has a viewing angle of 31°. The panels depict the same parameters as described in Fig. 6.

Fig. 8. STOKES modeling for IRAS 04416+1215 with only the equatorial scattering region and comparison with observational estimates. The case shown has a viewing angle of 4°. The panels depict the same parameters as described in Fig. 6.

The inclusion of the polar scatterers provided an additional flux close to the line center that the equatorial scatterers are not sufficient in providing – at least in the case of Mrk 1044 and SDSS J080101.41+184840.7. Nevertheless, the equatorial scatterers are necessary to recover wider line profiles. For comparison purposes, we show the results of our STOKES modeling for Mrk 1044 and SDSS J080101.41+184840.7 when excluding a polar scattering region in Figs. B.1 and B.2. Noticeably, the models without the polar scatterers are unable to recover the “peaky” profiles obtained in the observed spectrum, both in natural and polarized light. In addition, the polarization fraction is substantially underpredicted in these models, especially at the position marking the line center for Hα. We show the case with both the polar and equatorial scattering regions for IRAS 04416+1215 in Fig. B.3. As shown in this figure, the polarization level becomes significantly lower when the polar scatterers are included in the STOKES modeling. We note that the inclination angles used to compare the models to the observations in each of these plots were set as per the estimated value from the observations. For the cases of the SDSS J080101.41+184840.7 and IRAS 04416+1215, the wings in the full profile are much wider than predicted in the models. Larger viewing angles can account for these deficits. We discuss this aspect in the Sect. 6.3.

To disentangle the contribution from the polar scatterers to our modeling, we performed an additional test. Keeping the optical depth for the equatorial scattering region unchanged for our three sources (see values in the last column in Table 7), we prepared STOKES models for a range of τ_pol/τ_eq = [10⁻², 10] (i.e., the ratio of the optical depth in the polar scattering region to the optical depth in the equatorial scattering region⁷. The results are shown in Fig. 9. The left panels in this figure show the reduced-χ² distributions⁸. These reduced-χ² distributions show multiple minima. We mark the location of the minimum χ² and the second-to-minimum χ² for each case. The adjacent panels show the fit for the total flux (F_λ) corresponding to the two cases of τ_pol/τ_eq corresponding to the two minimum χ² results. The remaining parameters (p × F_λ, P and the χ) show similar behavior. The two minima are comparable, suggesting that the best fit can be obtained with a polar region whose optical depth is either low (∼0.01−0.032) or moderate (∼0.562−1.0, where 1.0 indicates the same optical depth as that of the equatorial scattering region). Therefore, in summary, the presence of a polar scatterer is quite certain since the representative solutions indicate a non-zero value for the τ_pol/τ_eq. We however note that obtaining a good fit for all the observables for our sources is not trivial and is affected by the noise in the observed spectra.

Fig. 9. Distribution of χ2 for a range of τpol/τeq for our three sources, Mrk 1044 (top rows), SDSS J080101 (middle rows), and IRAS 04416 (bottom rows). This range is the ratio of the optical depth in the polar region to the equatorial region. The location and the corresponding value of the ratio (τpol/τeq) of the minimum and second-to-minimum χ2 values are marked with dashed and dotted lines, respectively. Middle and right columns: fit for the total flux for each source corresponding to the minimum χ2 and the second-to-minimum χ2.

6. Discussion and conclusions

We performed VLT-FORS2 spectropolarimetric observations of three AGNs that belong to the NLSy1 class and were claimed to accrete at super-Eddington rates. Our goal was to check whether the Eddington rates are not overestimated by the specific orientation of the sources. Polarimetric measurements, allow a unique way to estimate the viewing angle for the observed sources, even for top-view sources, and can reveal the true range of the Keplerian velocities of the BLR.

The viewing angle of 54° for Mkn 1044 measured in this paper when comparing the natural light width and the width in the polarized light is comparable to i = $47.2_{-2.5}^{+1.0}$ obtained by Mallick et al. (2018) from the joint fitting of Swift, XMM-Newton, and NuSTAR X-ray spectra for Mrk 1044.

The viewing angle of SDSS J080101 using our method is 31°, and we did not find any independent measurements of the viewing angle for this source in the literature. As for IRAS 04416+1215, our determination of the viewing angle gave a very small value of 4°. Our formal (systematic) error on this value is small, but what we do not include here is a possible systematic error due to uncertainty in the H/R factor in Eq. (6). On the other hand, Tortosa et al. (2022) analyzed XMM-Newton and NuSTAR observations and were only able to get a highly model-dependent upper limit for the viewing angle for this source in the range of 25−45°. If the viewing angle is indeed higher than 4°, it would mean that the random motions in this source are less important (i.e., H/R in Eq. (6) is smaller than the one-third value adopted in the current paper). We discuss the issue of the viewing angle of IRAS 04416+1215 more in Sect. 6.3.

Our results confirmed that the sources are high accretors. The black hole mass values we obtained for each of the three sources from the polarized spectra are not systematically higher than the masses obtained earlier using reverberation mapping (see Table 5). The high accretor character of all three sources is also supported by our analysis of their X-ray spectra, which revealed features characteristic of high Eddington ratio sources: steep X-ray spectra and large soft X-ray excess (Tortosa et al. 2022; Mallick et al. 2018; Liu et al. 2021).

6.1. Polarization level

The overall level of polarization in the three sources is low, 0.5% or less, typical for Type 1 AGNs (see Smith et al. 2004; Robinson et al. 2011; Afanasiev et al. 2019).

The polarization in Mrk 1044 has been measured previously by Grupe et al. (1998), who gave a polarization continuum of 0.52% ± 0.05 and a polarization angle of 144°. This is roughly consistent with our results for the continuum (although our value is even lower, 0.23% ± 0.04). Their observations were performed with a 2.1 m telescope at the McDonald Observatory, which did not allow for specific study of the wavelength-dependent polarization. In general, a higher value for continuum polarization (Goodrich 1989) may be due to a wavelength range that is different from our measurements, as we do not have data for the rising blue side of the spectra.

Goodrich (1989) reports spectropolarimetric observations of Mrk 1044 and the general effect of degree of polarization rising toward the blue side. We did not notice this effect, but our observations have different wavelength ranges (for the Mrk 1044 observation in Goodrich (1989), it is 4436−7210 Å, whereas our measurements start from ∼6000 Å). The mean degree of polarization reported for Mrk 1044 in our work is 0.15% ± 0.01%, whereas in Goodrich (1989), the mean degree of polarization is 0.46% ± 0.05%, which confirms the low level of degree of polarization.

After analyzing measurements for a sample of 16 NLSy1, Robinson et al. (2011) concluded that the spectropolarimetric properties have a physical rather than a geometrical origin. NLSy1 galaxies tend to have a high Eddington ratio and high accretors that show more extreme physical conditions, such as higher metallicities (e.g., Shin et al. 2013; Sulentic et al. 2014; Panda et al. 2020) Moreover, Robinson et al. (2011) reported a peculiar, prominent red wing in Hα in the polarized flux of Mrk 1239, a characteristic of a polar scattering outflow. We did not notice redshifted asymmetries in the polarized Hα spectra of our sample.

The polarization measurements made by Lira et al. (2021) using ESO/FORS1 observations of NGC 3783, which has a low mass (2.4 − 3.5 × 10⁷ M_⊙; Onken & Peterson 2002), and Mrk 509, which has a relatively high mass (1.4 × 10⁷ M_⊙; Peterson et al. 2004), exhibit polarization levels of ∼0.7% and ∼0.9, respectively. These are low Eddington ratio sources that show relatively high levels of polarization, unlike the three sources studied in this work, even after correcting for contamination from the interstellar medium (ISM).

Older papers by Smith et al. (2004, 2005) contained many objects with polarization larger than 0.5%. The recent sample by Capetti et al. (2021) shows the mean level of polarization at 0.75% (the median is at 0.59%), but their sources do not preferentially belong to the high Eddington ratio class. What seems to be a difference between the sources of previous spectropolarimetry sources studied in previous surveys, such as the ones of Afanasiev & Popović (2015), Afanasiev et al. (2019), is the requirement of an additional polar scatterer. We hypothesize that the occurrence of a polar scatterer could be a genuine difference between sources radiating at high Eddington ratios and sources radiating at lower Eddington ratios (Sulentic et al. 2000b).

The low polarization in our sample requires the presence of both the equatorial and polar scatterers for the two objects Mrk 1044 and SDSS J080101.41+184840.7. In the third object, the polar scatter is not required, but it is not excluded. High Eddington ratio sources in general are expected to have more massive outflows that might lead to higher polarization rather than lower polarization, but this is not what we observed. The upper limits for the optical depth of the polar scatter are quite high, but the geometric location of the scatterer along the symmetry axis apparently does not imprint high polarization.

6.2. Comparison of sources’ properties with the sample from Capetti et al. (2021)

We decided to compare the spectral properties of our sample to a larger set of spectropolarimetric measurements which contains 25 AGN from the same instrument carried out recently by Capetti et al. (2021, 24 objects) and by Jiang et al. (2021, 1 object). In the paper of Capetti et al. (2021), the whole sample is a set of 25 objects. However, the BLR polarization measurements of source J145108 were affected by a cosmic ray, as the authors mentioned, and we did not include this source. For J140700, two observations were performed, and we chose an observation with a higher S/N = 413 (per resolution element at 6750 Å).

Additionally, the P_line from this work is the same quantity as P_BLR from Capetti et al. (2021). To illustrate this, we calculated the window from the definition postulated in Capetti et al. (2021) for Mrk 1044. The FWHM(Hα) from our spectral fitting is 1290 km s⁻¹, which gives 28 Å. After doubling this value, as in Capetti et al. (2021), we obtained 56 Å, which is comparable with the value of ≈50 Å in this work.

We decided to use the P_line term since, as we show for our objects, they may contain a narrow line component in this part of the polarized spectrum. In Fig. 10, we show three plots of the P_line/P_cont − λ_Edd plane⁹. Each plot contains a sample from Capetti et al. (2021), Fairall 9 measurements from Jiang et al. (2021), and measurements for our sources. For Fairall 9 and our sources, we labeled the points. We color coded the P_line/P_cont Eddington ratio (λ_Edd) to identify three quantities: the viewing angle, FWHM Hα, and R_FeII. To calculate P_line/P_cont from our sample, we used measurements from Table 3.

Fig. 10. Plane of Pline/Pcont − λEdd for sources from Capetti et al. (2021, marked as dots), Fairall 9 (labeled and marked as a star), and source from this work (with labels). Left panel: the color code is based on the viewing angle. We indicate Fairall 9 with the viewing angle measurement from Jiang et al. (2021) with a star and label it as Fairall (II). We indicate Fairall 9 with the viewing angle measurement based on the method from this work as a dot and label it as Fairall (I). The black point is the source J140336 from the sample of Capetti et al. (2021). We do not provide the viewing angle for this source since the FWHMunpolarized > FWHMpolarized (see Eq. (6)). Middle: the same sample color coded with the FWHM of a broad Hα. The sources from Capetti et al. (2021) are plotted with a broad Hα measurement from Capetti et al. (2021), and Fairall 9 is plotted with a broad Hα measurement from Jiang et al. (2021). In the case of our sources (the three sources observed with VLT), we used measurements from Table 4. Right panel: the same sample color coded with RFeII. The sources from Capetti et al. (2021) with RFeII were calculated from Shen et al. (2011) values. The black points are sources: J154743 (no EW(Fe II) measurement in Shen et al. 2011) and J084600 (no EW(Hβ) BC and no EW(Fe II) measurement in Shen et al. 2011).

In order to calculate the Eddington ratio of the sources from our sample, we used the following equations. For bolometric luminosity, we used the bolometric correction of Netzer (2019, Eq. (3) therein), which in our case (for L_5100 Å) is

$$\begin{array}{c}\hfill k_{\mathrm{BOL}}=40{\left(\frac{L_{5100\mathrm{\AA }}}{{10}^{42}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{-0.2}.\end{array}$$(11)

For Eddington luminosity, we used the masses that we calculated and present in Table 5.

For color coding, we used the viewing angle and FWHM(Hα) BC from Table 4 and the R_FeII from Table 6.

We plot in Fig. 10 the sample from Capetti et al. (2021) using the P_line/P_cont values from Table 3 in their paper (ratio of Cols. 5 and 2). For our sample, we used values from Table 3 (for Mrk 1044, P_line/P_cont = 1.59; for SDSS, J080101.41+184840.7 P_line/P_cont = 0.55; and for IRAS 04416+1215, P_line/P_cont = 0.19).

For Fairall 9, we used values from the analysis of Jiang et al. (2021): P_line/P_cont = 0.89 (P_line = 0.94, P_cont = 1.06). To calculate the Eddington ratio of the sources using the same method as our sample, we used M_BH and L_5100 Å from Table 1 from Capetti et al. (2021). For the color coding, we used the viewing angle we calculated from Eq. (6), wherein for the V_Kep we input values of inter-quartile width (between 25% and 75%) of the BLR in total intensity from Capetti et al. (2021; Table 3, Col. 8 in their paper), and for ΔV_obs, we input the FWHM of non-polarized Hα BC (Capetti et al. 2021; Table 1, Col. 7 in their paper). We also used the FWHM of polarized Hα BC (Table 3, Col. 9 from their paper). Finally, we calculated R_FeII from EW(Fe II) and EW(Hβ) BC from the catalog of Shen et al. (2011) for each source from the sample of Capetti et al. (2021).

For Fairall 9, Jiang et al. (2021) determined the viewing angle of the source as 50° based on the polarization level and polarization angle fit to the data. For our sources, we determined this viewing angle from a simple comparison of the FWHMs in natural and polarized light, following Eq. (6). To achieve consistency with our sample and to gain insight into the accuracy of the viewing angle determination, we used the data from Jiang et al. (2021) and adopted our method, based on Eq. (6). In this case, we obtained the value ∼27° for the viewing angle, which is smaller compared to the original determination by Jiang et al. (2021). We used both values in the plot, marking them as Fairall 9, our value, and Fairall (II), the original value. However, as Jiang et al. (2021) have commented, if the object has a more complex structure (i.e., a warped structure), a single value of the viewing angle may only be illusory. We computed the Eddington ratio with the method used for the rest of the sample. We used values for L_5100 Å and M_BH from Du et al. (2015). This value of the Eddington ratio is plotted in the middle and right panels of Fig. 10.

The outlying source on the P_line/P_cont − λ_Edd diagram with the highest P_line/P_cont ratio is J110538. This high ratio is caused by the extremely low continuum polarization in comparison to the rest of the (Capetti et al. 2021) sample. In the left panel, color coded with the viewing angle, we do not have the viewing angle calculation for J140336 since the FWHM of non-polarized Hα is broader than the polarized FWHM, thus making Eq. (6) unusable for this case. We mark this source as a black dot in this plot. For the right panel, we color-coded with R_FeII, and after cross match g with the catalog of Shen et al. (2011), we noticed no EW(Fe II) measurement for J154743 and no EW(Hβ) BC measurement for J084600. After visual inspection of those spectra, we saw no strong Fe II emission for J154743 nor for J084600, and we decided not to fit the Fe II template in order to avoid biasing the sample.

A quick look at those three panels does not reveal any obvious trend for the color-coded values. We suppose there is a weak trend with decreasing non-polarized FWHM of Hα along the x-axis (i.e., Eddington ratio). However, to confirm this trend (and others), we performed a principal component analysis for the full sample, which is shown in the next section. The increasing R_FeII along the Eddington ratio is a confirmation of the expected behavior from many previous works (Sulentic & Marziani 2015; Panda et al. 2019c).

6.3. Searching for the driver with principal component analysis

Principal component analysis (hereafter, PCA) is a “linear” dimensionality reduction technique where the dataset is projected onto a set of orthogonal axes that explain the maximum amount of variability in the dataset. The PCA technique works by initially finding the principal axis along which the variance in the multidimensional space (corresponding to all recorded properties) is maximized. This axis is known as “eigenvector 1”. Subsequent orthogonal eigenvectors, in order of decreasing variance along their respective directions, are found, until the entire parameter space is spanned (see, for example, Boroson & Green 1992; Kuraszkiewicz et al. 2009; Wildy et al. 2019; Martínez-Aldama et al. 2021). The PCA method is particularly useful when the variables within a dataset are highly correlated. Correlation indicates that there is redundancy in the data. Due to this redundancy, PCA can be used to reduce the original variables into a smaller number of new variables (principal components, or PCs), explaining most of the variance in the original variables. This approach allows us to determine correlated parameters, and in the context of our work, we utilized this technique to determine the physical parameter(s) that lead to the correlations shown in Fig. 10.

We used the sources shown in Table 9 for our PCA. We administered the PCA on the sample using only the direct observables, that is, the FWHM of the Hα emission line for both the polarized and the unpolarized profiles, the AGN monochromatic luminosity at 5100 Å (L_5100 Å), and the equivalent widths (EW) of the Hβ and the optical FeII blend within the 4434−4684 Å that is key to the optical plane of the quasar main sequence (Boroson & Green 1992; Sulentic et al. 2000b; Shen & Ho 2014; Marziani et al. 2018; Panda et al. 2018, 2019a,c). This choice of using only the direct observables allowed us to remove redundancies in the PCA that enter when the parameters that are derived (e.g., viewing angle, black hole mass, and Eddington ratio) using the direct observables are incorporated within the analysis (see Martínez-Aldama et al. 2021, for a detailed study on this issue). Thus, in our PCA, we used 25 sources with five properties.

Similar to Wildy et al. (2019), Martínez-Aldama et al. (2021), we used the prcomp module in the R statistical programming software. In addition to prcomp, we used the factoextra¹⁰ package to visualize the multivariate data at hand and especially to extract and visualize the eigenvalues and variances of the dimensions.

There is no well-accepted way to decide how many principal components are enough, we used the “scree plot” method to evaluate the number of principal components that best describe the variance in our dataset. A scree plot shows the variances (in percentages) against the numbered principal component and allows the number of significant principal components in the data to be visualized. The number of components is determined at the point beyond which the remaining eigenvalues are all relatively small and of comparable size (Peres-Neto et al. 2005; Jolliffe 2011).

The top panel in Fig. 11 shows the scree plot obtained from our PCA. It can be seen that the first four principal components explain 96.7% of the total variance in our dataset, with the first principal component (PC1) contributing 41.5% to this. The subsequent principal components (PC2, PC3, and PC4) account for 26%, 18.8%, and 10.4%, respectively. The first two principal components are sufficient to explain a major fraction of the variance (together they amount to 67.5%), yet we considered the first four principal components to obtain a more consistent picture for the sample using the PCA. The lower panels of Fig. 11 show the contributions of the five input properties to the first two principal components. The red dashed line on each of these sub-panels indicates the threshold (i.e., the expected average contribution). If the contribution of the variables were uniform, the expected value would be 1/length(variables) = 1/5 (i.e., 20%). For a given component, a variable with a contribution larger than this cutoff could be considered important and as contributing to the component. We thus note that the two FWHMs (polarized and then unpolarized) are the two dominant parameters contributing ∼75% to the PC1. The remaining parameters are below the threshold. Similarly, for the PC2, the EW(FeII) dominates, followed by EW(Hβ).

Fig. 11. Contribution of principal components. Top panel: scree plot using the PCA and showing the contribution (as a percentage) of the five principal components (PCs) to the overall variance in the dataset. Bottom panels: contributions of the original variables: the AGN luminosity at 5100 Å (L5100 Å, in units of erg s−1); the FWHM of the polarized and unpolarized Hα emission line (in units of km s−1); and the EW for the corresponding Hβ and optical Fe II emission within 4434−4684 Å (in units of Å) to the first two PCs. The red dashed line indicates the expected average contribution.

We show the 2D projection maps from our PCA in Fig. 12. We only show the projection map from the first two PCs (PC1 and PC2). The left panels in this figure show the projection maps with the corresponding PCs where their contribution to the total variance is also reported within parentheses. The 25 sources are shown along with the five input properties that are represented in the form of vectors. The magnitude and direction of these vectors show their importance to the respective PCs. The sources are colored as a function of the “squared cosine” (cos2) that shows the importance of a principal component for a given observation. It indicates the contribution of a component to the squared distance of the observation to the origin and corresponds to the square of the cosine of the angle from the right triangle made with the origin, the observation, and its projection on the component (see Abdi & Williams 2010, for an in-depth review on PCA). The panels on the right show the same projections but mark the source number and locate our three sources (Mrk 1044, SDSS J080101.41+184840.7, and IRAS 04416+1215) and Fairall 9. The PCA also serves as a clustering technique that can allow for the identification of groups of sources that have similar physical properties in multidimensional space. We derived the geometrical distance between every pair of observations on each projection map (PC1 to PC5) and retrieved two unambiguous mini-clusters, namely, (a) #10 (SDSS J113422.47+041127.7) and #17 (SDSS J154007.84+141137.0), and (b) #25 (IRAS 04416+1215) and (blue) #12 (SDSS J140621.89+222346.5). We set an upper limit for this distance (δr) at 0.57 to identify the mini-clusters and discard the imposters (see Fig. 12). These projection maps allowed us to scrutinize the real mini-clusters of sources from the imposters. What we mean by this is that if one focuses only on the dominant projection map, that is the PC1–PC2 plane, one may end up marking more mini-clusters. But the majority of these mini-clusters were found to be imposters – their relative positions vary significantly when traced on the remaining projections. The two pairs of sources have many identical/similar properties. For example, the first pair (#10 and #17) have identical luminosities, log L_5100 Å = 42.9 (in erg s⁻¹) and almost identical black hole masses, log M_BH = 8.4 and 8.46 (in M_⊙), in addition to similar EW(Hβ) and EW(FeII). For the second pair (#12 and #25) which includes one of our sources (i.e., IRAS 04416+1215), we found similar luminosities, log L_5100 Å = 44.2 and 44.47; black holes masses close to each other, log M_BH = 7.36 and 6.97; and similar EW(Hβ) and EW(FeII). The sources in each pair also have similar Eddington ratios as well as similar R_FeII values (see Table 9).

Fig. 12. Two-dimensional projections of the first two principal components (PCs) using the PCA for the sources shown in Table 9 after filtering out the sources marked with an asterisk. Only direct observables from the spectra were used in the PCA, namely, the AGN luminosity at 5100 Å (L5100 Å, in units of erg s−1), the FWHM of the polarized and unpolarized Hα emission line (in units of km s−1), and the EW for the corresponding Hβ and optical Fe II emission within 4434−4684 Å (in units of Å). The percentages reported in parenthesis mark the contribution of the respective PC to the overall variance in the data. The importance of the direct observables to each PC is shown using vectors. The sources are color coded according to their “squared cosine” values and are a function of the squared distance of the observation to the origin. The sources are numbered as in Table 9. The location of our three sources (Mrk 1044, SDSS J080101.41+184840.7, and, IRAS 04416+1215) and Fairall 9 are indicated. The two mini-clusters containing pairs of sources that were revealed from our clustering selection are shown with green shaded ovals. The clustering selection was based on the geometrical distance between two data points, δr ≲ 0.57, across all the projections.

We show the spectra of both pairs of objects in Fig. 13. They indeed show similar spectral properties, such as the shapes of the continuum, Balmer lines (Hα, Hβ, Hγ), [OIII] doublet, and Fe II pseudo-continuum emission. Also, they have a similar redshift, namely, z = 0.119 for J154007 and z = 0.108 for J113422. For the second pair, z = 0.097 for J140621 and z = 0.089 for IRAS 04416. The overall similarity of the spectra and the derived parameters thus suggests that an upper limit for the viewing angle for IRAS 04416+1215 can be ∼34°. We note that assuming this value of the inclination angle gives us a BH mass, log M_BH ∼ 6.397, which is smaller by a factor of three relative to the BH mass estimated assuming inclination angle of 4° (see Table 5). Such an analysis can be useful for larger samples to identify sources with similar spectral and physical properties. In Fig. B.4, we show the STOKES modeling and a comparison with the observed spectra for IRAS 04416+1215 assuming an inclination angle of 35° and including the polar scattering region. We note that the model with only the equatorial scattering region gives qualitatively similar results to Fig. B.4, although the polarized flux has a significantly lower flux level at the line center (highlighted with the red dashed line). The modeled spectrum in the natural light has a better agreement relative to the 4° case (see Fig. 8), although we need to assume an inclination closer to 54° to have the wings modeled with the lowest residuals. In our STOKES modeling, since we assume the FWHM that is estimated from the spectral fitting of the observed spectrum in the polarized light, which is quite noisy, the FWHM can be larger than the assumed value. Thus, we have a degeneracy between the estimation of the inclination angle and the FWHM of the polarized Hα spectrum. We need better-quality data to break this degeneracy.

Fig. 13. Spectra of the grouped objects due to our PCA analysis. In the left panel, we show (green) #10 (SDSS J113422.47+041127.7) and (blue) #17 (SDSS J154007.84+141137.0). In the right panel we show (green) #25 (IRAS 04416+1215) and (blue) #12 (SDSS J140621.89+222346.5). The vertical dotted line represents the Hβ and Hα rest-framed position.

Finally, we show the correlations of the important derived parameters for the first principal component to check the primary driver(s) of this sample. In Fig. 14, we show AGN luminosity at 5100 Å in addition to the derived parameters (i.e., black hole mass, Eddington ratio, R_FeII, and the viewing angle). We performed an ordinary least-square fit for each panel and report the Spearman’s rank correlation coefficient and the corresponding p-value for the panels only when the p-value was less than 0.001. This criterion was taken from our previous study (Martínez-Aldama et al. 2021), and it allowed us to identify true and robust correlations. We noticed that the first principal component (PC1) is primarily driven by a combination of black hole mass, Eddington ratio, and R_FeII. The remaining PCs do not show any significant correlations and hence are not shown. As noticed in earlier works, the parameter R_FeII is a primary driver of the eigenvector 1 in the optical plane of the quasar main sequence (Boroson & Green 1992; Sulentic et al. 2000b; Shen & Ho 2014; Marziani et al. 2018; Panda et al. 2018, 2019a,c) and is a direct tracer of the Eddington ratio. Smaller black hole masses and larger luminosities can push the Eddington ratios to higher values, and we observed this effect in our correlations. Also, such sources show stronger Fe II emissions. The correlation between the PC1 and black hole mass, and the anti-correlation between the PC1 and R_FeII highlight this trend appropriately. Our three sources are highlighted in green and occupy the region with the lowest black hole masses and highest Eddington ratios. Their L_5100 Å and R_FeII are comparable to some of the sources from Capetti et al. (2021), mostly those occupying the more luminous regions and with higher R_FeII values. In no situation did we observe the viewing angle (panels in the last column) to be significant in the PCA for this sample. A reason for this result could be the overall similarity of the two FWHMs (unpolarized and polarized) vectors as seen from the projection maps (see Fig. 12), and since the viewing angle was estimated using the ratio of these two FWHMs, the overall correlation is affected. The FWHMs individually have strong correlations with the PCs, especially with the PC1 (see Fig. 11), which dampens the correlation with the viewing angle.

Fig. 14. Correlations between the first principal component (PC1) and the physical properties of the sample. From left to right: correlation with the AGN luminosity at 5100 Å (L5100 Å, in units of erg s−1), the black hole mass (MBH, in units of M⊙), the Eddington ratio, the ratio RFeII, and the viewing angle (i). Our three sources (Mrk 1044, SDSS J080101.41+184840.7, and, IRAS 04416+1215) are marked in green, while the other sources are marked in red. The Spearman’s rank correlation coefficients (ρ) and the p-values are reported for the correlations whenever the p-value < 0.001. The ordinary least-square fit for each panel is shown using a red solid line. The black dotted lines mark the confidence intervals at 95% for the 1000 realizations (dark gray lines) of the bootstrap analysis. The corresponding prediction intervals are shown in the background using a light gray color.

6.4. Summary

In this work, we present new VLT-FORS2 spectropolarimetric measurements for three objects selected as high Eddington ratio candidates. Our findings from this work are as follows:

• The viewing angles recovered for the three sources indicate that they cover a large range in the viewing angle values, from an almost face-on orientation (IRAS 04416+1215) through an intermediate case (SDSS J080101.41+184840.7) to a highly inclined orientation (Mrk 1044). We find it important to highlight that the viewing angle estimation for IRAS 04416+1215 should be taken with caution.

• Despite the large differences in the viewing angles, we confirm the small values of the black hole mass in these sources and their high Eddington nature.

• We were successful in recovering the observed Hα line profile both in the natural and polarized light using the STOKES modeling. We recovered the polarization fractions of the order of 0.2−0.5% for the three sources, although the recovery of the phase angle is sub-optimal, mainly due to the noise in the observed data.

• Our principal component analysis shows that the sample of the 25 sources, including our sources, is mainly driven by the black hole mass and Eddington ratio. We reaffirm the connection of the strength of the optical Fe II emission with the Eddington ratio, but the dependence on the viewing angle is ordinary and resembles more of a secondary effect.

Acknowledgments

We are very grateful to the Referee for very helpful and constructive remarks which help to improve considerably the content of the manuscript. M.S. acknowledges support from Polish Funding Agency National Science Centre, project 2021/41/N/ST9/02280 (PRELUDIUM 20), the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement 950533) and from the Israel Science Foundation, (grant 1849/19). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. [951549]). M.S. is grateful to S. Bagnulo for the discussion about data reduction and all tutors and organizers of Torun Summer School 2019: Polarimetry as a diagnostic tool in astronomy for providing a review of methods and techniques developed for the polarimetry data. We also thank organizers and participants of the conference AGN and Polarimetry conference: Towards a panchromatic understanding of the polarization of Active Galactic Nuclei for discussion which helps to improve this work. The project was partially supported by the Polish Funding Agency National Science Centre, project 2017/26/A/ST9/00756 (MAESTRO 9), and MNiSW grant DIR/WK/2018/12. S.P. acknowledges the financial support of the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) Fellowships 164753/2020-6 and 300936/2023-0. Đ.S. acknowledges support by the Ministry of Education, Science and Technological Development of the Republic of Serbia through contract no. 451-03-9/2021-14/200002 and by the Science Fund of the Republic of Serbia, PROMIS 6060916, BOWIE. M.L.M.A. acknowledges financial support from Millenium Nucleus NCN19-058 (TITANs). L.Č.P. is supported by the Ministry of Education, Science and Technological Development of Serbia, the number of the contract is 451-03-68/2020-14/200002 and L.Č.P. acknowledges the support by the Chinese Academy of Sciences President’s International Fellowship Initiative (PIFI) for visiting scientist. Funding for the Sloan Digital Sky Survey (SDSS) has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the US Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS website is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium (ARC) for the Participating Institutions. The Participating Institutions are The University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, The Johns Hopkins University, Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, University of Pittsburgh, Princeton University, the United States Naval Observatory, and the University of Washington. Authors’ contribution: M.S. performed the data reduction and spectral fitting of the observational data. S.P. performed the STOKES modeling and the PCA analysis. M.S., S.P., and B.C. primarily wrote the manuscript. Đ.S. helped with the STOKES modeling. M.L.M.A. helped with the spectral fitting. P.M. helped with text writing and computations. J.M.W., P.D., and L.C.P. helped with the discussions. C.S.S. helped with the data analysis of Planck-derived products. Software: MATPLOTLIB (Hunter 2007); NUMPY (Oliphant 2015); SCIPY (Virtanen et al. 2020); SKLEARN (Pedregosa et al. 2011); STATSMODELS (Seabold & Perktold 2010); STOKES (Marin 2018); HEALPIX (Górski et al. 2005); ASTROPY (Astropy Collaboration 2013).

References

Appendix A: Contamination of the polarized emission of studied sources with respect to the interstellar extinction and host galaxy

As clearly discussed and modeled by Marin (2018), in the case of observations of Seyfert 1 galaxies the contamination due to the ISM and dilution of the observed flux by the host can considerably affect the polarization measurements. Dust scattering in our Galaxy along the line of sight from the source imprints its own polarization signature that is predominantly set by the total dust column. The maximum of the effect is seen at ∼5450 Å (Serkowski et al. 1975). The ISM polarization in our Galaxy is primarily due to dust dichroic extinction, that is, differential linear extinction for the two waves polarized along and perpendicular to the direction of alignment of the light propagation vector before reaching the ISM.

The best way to address the role of the ISM is to measure the observed polarization for several unpolarized stars located close to the selected object. Such an approach was adopted by Jiang et al. (2021) for Fairall 9, which was also observed with VLT/FORS2, and it showed that for this source the effect is unimportant. Unfortunately, for the sources discussed in the present paper, such observations were not performed, and the two comparison stars (polarized and unpolarized) were located far from any of the objects. The other way to estimate the possible role of interstellar polarization is through the measurement of the Galactic extinction in the direction of the source. We give these values, parameterized by A_V (taken from NED Database) in Table A.1.

We observed that for two of the sources, Mrk 1044 and SDSS J080101.41+184840.7, the extinction is not considerably higher than for Fairall 9, and the effect of the ISM can be neglected. However, this is not the case for IRAS 04416+1215, as the Stokes parameter Q changes sign and affects the measurements the most.

Therefore, we consulted the Planck polarization maps, as recommended by Capetti et al. (2021) and Pelgrims (2019). We used the intensity (I) and polarization (Q and U) thermal dust maps at the 353 GHz frequency channel from the 2018 Planck data release¹¹ described in Planck Collaboration IV (2020). These dust maps were produced using the Generalized Needlet Internal Linear Combination (GNILC) component separation method, which uses spectral information as well as the angular power spectrum to disentangle specific diffuse foreground components. The Stokes parameters I, Q, and U from the Planck datasets are in units of the thermodynamic temperature K_CMB¹², in milli-Kelvins. We converted these Stokes parameters to MJy sr⁻¹, assuming the CMB temperature at z = 0 to be 2.73 K. The values of the polarization level and the polarization angle of the foreground thermal dust contribution at the position of the discussed sources are given in Table A.2, and their positions are indicated in Figures A.1 and A.2. We estimated the polarization contribution from the ISM for sources from our sample and for Fairall 9. To correct the Stokes parameters computed from optical observation, we first estimated the contribution of the ISM from submillimeters to V-band using a prescription from Pelgrims (2019), Capetti et al. (2021) and values from Table A.1 using the following equation:

$$\begin{array}{c}\hfill Q_{\mathit{CORR}}=-0.238\ast A_V\ast Q_{\mathit{GNILC}}\end{array}$$(A.1)

Fig. A.1. Degree of polarization (P in %) computed using the GNILC 353GHz thermal dust map from Planck Collaboration IV (2020) at HEALPIX resolution, Nside = 2048. The location of our three sources (black dots) and Fairall 9 (black star) are indicated based on their galactic coordinates reported in Table A.1. The values for P for these sources are reported in Table A.2.

Fig. A.2. Same as Figure A.1 but for the polarization angle (χ) computed using Equation 5. The values for χ for these sources are reported in Table A.2.

which has a similar expression for U_CORR. Next, we corrected our observed values of the Stokes parameters and computed P_corr and χ_corr, which we report in Table A.2. For our three sources, Mrk 1044 and SDSS J080101.41+184840.7, and Fairall 9, there is no significant change in the χ_corr and P_corr. For Mrk 1044, the degree of polarization before correction is 0.15 ± 0.02%, and after correction it is 0.19%. For SDSS J080101.41+184840.7, the degree of polarization before correction is 0.52%, and after correction it is 0.49%. Hence, the ISM correction does not change our original inferences. For IRAS 04416+1215, however, we observed a noticeable change after the correction: 0.17% before the correction and 0.64% after the correction. Compounded with other degeneracies, such as problems with the estimation of the viewing angle (see Section 6.3 for details), this result makes it difficult to draw any clear conclusion for IRAS 04416+1215 with the available data.

Robinson et al. (2011) argues that variations of P and χ across the Balmer lines indicate that sources are also intrinsically polarized. We did not see such signatures in our sources, but it may be due to short exposures.

In the case of the host galaxy’s contamination of our results, we compared the shape of the degree of polarization along the wavelength from our measurements to models from Marin (2018; Figure 3 therein). Models from Marin (2018) suggest a small decrease in the degree of polarization along the wavelength in the case of contamination from the ISM or host galaxy. For our sources, however, we noticed a rather small increase in the degree of polarization along the wavelength. This small contribution to the measured flux comes from the fact that in our VLT/FORS2 observations, we used a very narrow slit of 0.7″, so the spectrum contains much less starlight than the typical measurements available in the NED Database¹³ that reflect a red spectrum contaminated by the host.

We performed an additional test using the available acquisition image for Mrk 1044. The exposure time was only ∼1s (in OG590 filter) but the quality of the image was good. We analyzed the image using the method described in Bentz et al. (2013) for the starlight subtraction using Hubble Space Telescope images. We fitted the background and the active galaxy profile assuming a Gaussian profile for the point spread function and the Sersic profile for the host, with a fixed value of the profile index n. After fitting the parameters we integrated the profiles for the long slit used in the data acquisition (0.7 arcsec). We assumed Sersic profiles for two cases of Sersic index n: n = 2 and n = 4. The contamination weakly depended on the value of n and was about 0.1852 in the case of n = 2, and 0.1858 in the case of n = 4. The numbers represent the fraction of the host galaxy flux to the AGN flux (which is dimensionless). Hence, the contribution from the host galaxy is not dominant.

Appendix B: Supplementary plots

Fig. B.1. STOKES modeling and comparison with observational estimates for Mrk 1044 with only the equatorial scattering region included. The case shown has a viewing angle of 54°. From top to bottom, the panels show the spectrum in natural light (Fλ), the polarized spectrum (p × Fλ), the polarization fraction (p) and the polarization angle (χ). The vertical red dashed line marks the central wavelength for Hα.

Fig. B.2. STOKES modeling and comparison with observational estimates for SDSS J080101.41+184840.7 with only the equatorial scattering region included. The case shown has a viewing angle of 31°. The panels are similar to Figure B.1.

Fig. B.3. STOKES modeling and comparison with observational estimates for IRAS 04416+1215 with the equatorial plus polar scattering region included. The case shown has a viewing angle of 4°. The panels are similar to Figure B.1.

Fig. B.4. STOKES modeling and comparison with observational estimates for IRAS 04416+1215 with the equatorial plus polar scattering regions included. The case shown has a viewing angle of 35°. The panels are similar to Figure B.1.

Fig. B.5. Corollary plots for Figure 10.

All Tables

All Figures

Fig. 1. Natural light spectra for Mrk 1044 (top row), SDSS J080101.41+184840.7 (middle row), and IRAS 04416+1215 (bottom row). Top row: normalized 6dF spectrum of Mrk 1044 from Jones et al. (2009) with the fitted Hβ (left) and Hα (middle) part. The right plot shows the normalized unpolarized spectrum obtained with FORS2/VLT. The middle and lower rows show as follows for SDSS J080101.41+184840.7 and IRAS 04416+1215, respectively: fitted Hβ and Hα part from SDSS normalized archival spectrum and normalized unpolarized spectrum obtained with FORS2/VLT. Data are marked in gray, and the model is marked in black. The fit components are given in the following: for Hβ (left panel), the Hβ BC and [O III] semi-broad component doublet are red. The NCs of the Hβ and [OIII] doublet are black. The Hβ blue component is marked in blue. The power-law component is marked in blue, and the Fe II pseudo continuum is green. For Hα (middle and right panel), the Hα NC and [NII] and [SII] doublets are marked in black. The Hα BC is marked as red. The Hα blue component is marked in blue, the Fe II pseudo continuum is marked in green, and the power-law component is marked in blue. The lower panels of each plot correspond to the residuals, in radial velocity units km s−1. All spectra are normalized to Hα intensity.

In the text

Fig. 2. Data of Mrk 1044, SDSS J080101.41+184840.7, and IRAS 04416+1215. Plots, from top to bottom: total flux, Stokes parameters U and Q, flux of U and Q, and the polarization angle. We mark in red the windows in which we calculated polarization. For continuum polarization, we considered two surrounding Hα from its blue and red side, each window has a width of ≈100 Å, and they are at a distance of at least 300 Å from the center of Hα. For the polarization of the line, we considered the red window mark on the line. Each central window has width of ≈50 Å. The vertical blue dotted line marks the central wavelength for Hα. The data were rebinned with a 2 Å binning.

In the text

Fig. 3. Fitted Hα part from the polarized spectrum obtained with FORS/VLT. The data are marked in gray, and the model is marked in black. For each fit, we considered two components of Hα: a BC (red) and an NC (black). The blue line corresponds to the power-law continuum. The lower panels of each plot correspond to the residuals, in radial velocity in units of km s−1. With an arrow, we mark the presence of an “absorption” in the fit of SDSS J080101. The “absorption” is meant to account for the low S/N of the data and to ease the fit of a symmetric Gaussian constrained by the blue side of the polarized flux profile.

In the text

	Fig. 4. Schematic view of the models with the equatorial scatterer. The continuum source is marked in gray, the BLR is marked in blue, and the equatorial scattering region is marked in red. We marked the half-opening angle of the scattering structure with a dashed curve between the gray and black dashed lines.
In the text

	Fig. 5. Schematic view of the models with two scatterers: polar and equatorial. The various components are colored identically to Fig. 4. The polar scattering region is marked in purple.
In the text

	Fig. 6. STOKES modeling for Mrk 1044 with the equatorial plus polar scattering region and comparison with observational estimates. The case shown has a viewing angle of 54°. From top to bottom: spectrum in natural light (Fλ), polarized spectrum (p × Fλ), polarization fraction (p), and polarization angle (χ). The vertical red dashed line marks the central wavelength for Hα.
In the text

	Fig. 7. STOKES modeling for SDSS J080101.41+184840.7 with the equatorial plus polar scattering region and comparison with observational estimates. The case shown has a viewing angle of 31°. The panels depict the same parameters as described in Fig. 6.
In the text

	Fig. 8. STOKES modeling for IRAS 04416+1215 with only the equatorial scattering region and comparison with observational estimates. The case shown has a viewing angle of 4°. The panels depict the same parameters as described in Fig. 6.
In the text

Fig. 9. Distribution of χ2 for a range of τpol/τeq for our three sources, Mrk 1044 (top rows), SDSS J080101 (middle rows), and IRAS 04416 (bottom rows). This range is the ratio of the optical depth in the polar region to the equatorial region. The location and the corresponding value of the ratio (τpol/τeq) of the minimum and second-to-minimum χ2 values are marked with dashed and dotted lines, respectively. Middle and right columns: fit for the total flux for each source corresponding to the minimum χ2 and the second-to-minimum χ2.

In the text

Fig. 10. Plane of Pline/Pcont − λEdd for sources from Capetti et al. (2021, marked as dots), Fairall 9 (labeled and marked as a star), and source from this work (with labels). Left panel: the color code is based on the viewing angle. We indicate Fairall 9 with the viewing angle measurement from Jiang et al. (2021) with a star and label it as Fairall (II). We indicate Fairall 9 with the viewing angle measurement based on the method from this work as a dot and label it as Fairall (I). The black point is the source J140336 from the sample of Capetti et al. (2021). We do not provide the viewing angle for this source since the FWHMunpolarized > FWHMpolarized (see Eq. (6)). Middle: the same sample color coded with the FWHM of a broad Hα. The sources from Capetti et al. (2021) are plotted with a broad Hα measurement from Capetti et al. (2021), and Fairall 9 is plotted with a broad Hα measurement from Jiang et al. (2021). In the case of our sources (the three sources observed with VLT), we used measurements from Table 4. Right panel: the same sample color coded with RFeII. The sources from Capetti et al. (2021) with RFeII were calculated from Shen et al. (2011) values. The black points are sources: J154743 (no EW(Fe II) measurement in Shen et al. 2011) and J084600 (no EW(Hβ) BC and no EW(Fe II) measurement in Shen et al. 2011).

In the text

Fig. 11. Contribution of principal components. Top panel: scree plot using the PCA and showing the contribution (as a percentage) of the five principal components (PCs) to the overall variance in the dataset. Bottom panels: contributions of the original variables: the AGN luminosity at 5100 Å (L5100 Å, in units of erg s−1); the FWHM of the polarized and unpolarized Hα emission line (in units of km s−1); and the EW for the corresponding Hβ and optical Fe II emission within 4434−4684 Å (in units of Å) to the first two PCs. The red dashed line indicates the expected average contribution.

In the text

Fig. 12. Two-dimensional projections of the first two principal components (PCs) using the PCA for the sources shown in Table 9 after filtering out the sources marked with an asterisk. Only direct observables from the spectra were used in the PCA, namely, the AGN luminosity at 5100 Å (L5100 Å, in units of erg s−1), the FWHM of the polarized and unpolarized Hα emission line (in units of km s−1), and the EW for the corresponding Hβ and optical Fe II emission within 4434−4684 Å (in units of Å). The percentages reported in parenthesis mark the contribution of the respective PC to the overall variance in the data. The importance of the direct observables to each PC is shown using vectors. The sources are color coded according to their “squared cosine” values and are a function of the squared distance of the observation to the origin. The sources are numbered as in Table 9. The location of our three sources (Mrk 1044, SDSS J080101.41+184840.7, and, IRAS 04416+1215) and Fairall 9 are indicated. The two mini-clusters containing pairs of sources that were revealed from our clustering selection are shown with green shaded ovals. The clustering selection was based on the geometrical distance between two data points, δr ≲ 0.57, across all the projections.

In the text

	Fig. 13. Spectra of the grouped objects due to our PCA analysis. In the left panel, we show (green) #10 (SDSS J113422.47+041127.7) and (blue) #17 (SDSS J154007.84+141137.0). In the right panel we show (green) #25 (IRAS 04416+1215) and (blue) #12 (SDSS J140621.89+222346.5). The vertical dotted line represents the Hβ and Hα rest-framed position.
In the text

Fig. 14. Correlations between the first principal component (PC1) and the physical properties of the sample. From left to right: correlation with the AGN luminosity at 5100 Å (L5100 Å, in units of erg s−1), the black hole mass (MBH, in units of M⊙), the Eddington ratio, the ratio RFeII, and the viewing angle (i). Our three sources (Mrk 1044, SDSS J080101.41+184840.7, and, IRAS 04416+1215) are marked in green, while the other sources are marked in red. The Spearman’s rank correlation coefficients (ρ) and the p-values are reported for the correlations whenever the p-value < 0.001. The ordinary least-square fit for each panel is shown using a red solid line. The black dotted lines mark the confidence intervals at 95% for the 1000 realizations (dark gray lines) of the bootstrap analysis. The corresponding prediction intervals are shown in the background using a light gray color.

In the text

	Fig. A.1. Degree of polarization (P in %) computed using the GNILC 353GHz thermal dust map from Planck Collaboration IV (2020) at HEALPIX resolution, Nside = 2048. The location of our three sources (black dots) and Fairall 9 (black star) are indicated based on their galactic coordinates reported in Table A.1. The values for P for these sources are reported in Table A.2.
In the text

	Fig. A.2. Same as Figure A.1 but for the polarization angle (χ) computed using Equation 5. The values for χ for these sources are reported in Table A.2.
In the text

Fig. B.1. STOKES modeling and comparison with observational estimates for Mrk 1044 with only the equatorial scattering region included. The case shown has a viewing angle of 54°. From top to bottom, the panels show the spectrum in natural light (Fλ), the polarized spectrum (p × Fλ), the polarization fraction (p) and the polarization angle (χ). The vertical red dashed line marks the central wavelength for Hα.

In the text

	Fig. B.2. STOKES modeling and comparison with observational estimates for SDSS J080101.41+184840.7 with only the equatorial scattering region included. The case shown has a viewing angle of 31°. The panels are similar to Figure B.1.
In the text

	Fig. B.3. STOKES modeling and comparison with observational estimates for IRAS 04416+1215 with the equatorial plus polar scattering region included. The case shown has a viewing angle of 4°. The panels are similar to Figure B.1.
In the text

	Fig. B.4. STOKES modeling and comparison with observational estimates for IRAS 04416+1215 with the equatorial plus polar scattering regions included. The case shown has a viewing angle of 35°. The panels are similar to Figure B.1.
In the text

	Fig. B.5. Corollary plots for Figure 10.
In the text