In order to analyze simultaneously all components of the integrated light from selected objects, we used the ULySS ¹ (Koleva et al., 2009) full spectrum fitting package, which fits spectroscopic observations against a bounded linear combination of nonlinear model components, convolved with a parametric line-of-sight velocity distribution.

The ULySS code was originally developed for studying the history of stellar populations and for analyzing stellar atmospheres. One of the main advantages of ULySS is its flexibility, allowing users to define arbitrary non-linear components, as explained in detail by Bon et al. (2014), Bon et al. (2016). In this study, we modified ULySS by introducing additional components designed to simultaneously fit the nebular continuum, stellar population templates, AGN continuum and a comprehensive set of emission lines. The parameters of the emission lines could be either tied or left free, depending on the specific fitting conditions.

For the purpose of studying spectra in NELG, the model of emission line galaxy used in the fit represents the linear combination of power law continuum, stellar population model convolved with a line-of-sight broadening function, and a sum of Gaussians, that represent emission lines, as described in Bon et al. (2014). The model is multiplied by a Legendre polynomial that supposes to absorb the errors in the flux calibration, Galactic extinction or any other cause that affects the shape of the continuum. The model is generated at the same resolution and with the same sampling as the observation. Spectra were fitted in the 3,700–6,800 Å wavelength range. The fit is performed in the pixel space.

We used Vazdekis stellar population models from the library of single stellar populations, computed with the MILES library and Salpeter IMF (Sánchez-Blázquez et al., 2006; Vazdekis et al., 2010).

By fitting a spectrum with the ULySS, we reconstruct the SSP-equivalent age and metallicity, mean stellar velocity and velocity dispersion, spectral index of the featureless continuum, shifts, widths and intensities of emission lines, as well as fractions of all model components.

For this study we fitted all Balmer lines and strong forbidden lines [O II] λ λ 3727,3729, [O III] λ λ 4959,5007, [N II] λ λ 6548,6585, and [S II] λ λ 6716,6732. We also included coronal lines like [Fe VII] λ 6,087 and [Fe X] λ 6,374, He II λ 4,686, and several other less prominent features when present.

In this work we analyzed only narrow emission lines, narrower than 600 km s^-1, so the profiles of emission lines were fitted with a single Gaussian function. Figure 1 represent examples of the fit in total fitted wavelength range, and zoomed H β and H α domains in Sy2 galaxy. In the upper panels of the graphs, the black line represents the observed spectrum, the blue line represents the best fit model, and the red line the multiplicative polynomial, while the green, light red, and violet lines represent components of the best fit model: violet–stellar population, red–emission lines, and green–AGN continuum. The bottom panel represents residuals of the best fit (black line). The solid green lines mark the 1 σ deviation, and the dashed line is the zero axis.

The common way of fitting [OIII] λ λ 4959,5007 and [NII] λ λ 6548,6583 lines is to tie their widths and to fix their intensities to ≈ 3. Here both width and intensity of Gaussians used to fit emission lines are free parameters in the fit, so we were able to analyze flux ratios of [OIII] λ λ 4959,5007 and [NII] λ λ 6548,6583 lines in four subclasses of NELG, and search for their differences.

Emission lines fluxes as well as their uncertainties were measured using IDL procedure INT_TABULATED from the best fit model of emission lines and residual of the fit, respectively.

Before measuring the flux ratios we removed from the set those spectra where: (a) relative error of line dispersion of any analyzed forbidden emission line is higher than 10%; (b) the difference between the width of two emission lines in the doublet [OIII] λ λ 4959,5007 lines is higher than 20%; the same criteria is used for two [NII] λ λ 6548,6583 lines; (c) relative error of the line flux is higher than 20%; and (d) overlap between [NII] 6,548 A ̊ and H α lines was higher than 1Å (to be as consistent as with the data set constrains of Dojčinović et al. (2023), where they selected spectra without overlapped emission lines).

In our analysis, each emission line was modeled using a single Gaussian profile. Even though the two lines of a forbidden doublet (e.g., [N II] λ λ 6548,6584 or [O III] λ λ 4959,5007 ) are physically expected to originate from the same ionized region and thus share similar kinematic properties, we chose not to tie their widths during the fitting process. Instead, we allowed the line widths to be fit independently and subsequently applied a post-fit consistency check: spectra in which the fitted widths of the two doublet components differed by more than 20% were excluded from the analysis. This approach was motivated by the relatively low spectral resolution of the SDSS data ( R ∼ 1800 , corresponding to ∼ 69–75 km s^-1), which can introduce fitting uncertainties, particularly when working with narrow and closely spaced lines. Enforcing tied widths under such conditions can sometimes lead to artificially poor fits or biased line fluxes due to local residuals or noise structures. Our adopted strategy allowed for greater fitting flexibility while ensuring that only physically plausible measurements-where the kinematics of the doublet components remain reasonably consistent-were retained for statistical analysis. We emphasize that the goal of this procedure was to retain empirical reliability without imposing potentially misleading constraints in low-resolution data.

We used this method for the statistical analysis of the spectra in the central kpc of the selected NELG. The narrow emission lines can be used to classify the dominant ionization energy source in emission line galaxies. Consequently we made the well known H α /[NII]6548,6583 A ̊ vs. H β /[OIII]4959,5007 A ̊ BPT diagnostic diagram (Baldwin et al., 1981) to separate four classes of objects - Seyfert 2, LINERs, HII regions and so-called” composite objects “(Kewley et al., 2001; Kauffmann et al., 2003). By analyzing the BPT diagram, we distinguished 226 Sy2s, 89 LINERs, 3,117 starburst galaxies and 1,078 composite objects.

Therefore, from 229 Sy2s, 89 LINERs, 3,117 starburst galaxies and 1,078 composite objects, these additional constraints have limited our samples to 37 Sy2, 48 LINERs, 262 composite spectra and 1,227 starbursts for [OIII] ratio analysis, while for the case of [NII] ratio analysis our analyzed samples counted 37 Sy2s, 47 LINERs, 202 composites and 1,513 starbursts for flux ratio analysis.

In order to remove extreme outliers we performed sigma-clipping procedure within standard 3 σ threshold. This process primarily affected the starburst sample, where approximately ∼ 2.5 % of objects were removed. In contrast, for other galaxy subclasses, the number of spectra remained largely unchanged after sigma clipping; either no objects were excluded, or at most one or two spectra were clipped. The final number of spectra, after sigma clipping, used for calculation of the arithmetic average and their errors are reported in the Table 1. After applying all our selection criteria and performing σ -clipping to remove outliers, we found that the broadest emission line in our sample exhibits a velocity dispersion of σ = 149 km s^-1, which corresponds to a full width at half maximum (FWHM) of 352 km s^-1, assuming a Gaussian profile. This value sets the upper limit for the line widths considered in our analysis and ensures consistency across all fitted components.

The observed emission line fluxes were corrected for interstellar reddening using the Balmer decrement, specifically the H α /H β ratio.

To accurately measure the intrinsic fluxes of emission lines in our sample of galaxies, we corrected for interstellar reddening by applying an extinction correction based on the observed Balmer decrement, specifically the H α / H β ratio. This approach allows us to estimate the degree of dust extinction, which primarily affects shorter wavelengths and thereby alters the observed flux ratios of emission lines. The reddening correction was applied using the Cardelli et al. (1989) extinction.

The hydrogen Balmer lines are commonly employed to estimate the level of extinction affecting the observed emission lines of a photoionized plasma. This is because the intrinsic Balmer decrement values are largely insensitive to variations in gas temperature and density (if electron density is low, ∼ 100   c m − 3 ). This principle holds true for H II regions that are ionized by hot stars. Deviations from this intrinsic value in the observed ratio indicate the presence of dust extinction. However, for the narrow-line regions (NLRs) of Seyfert 2 galaxies, where the ionizing spectrum is significantly harder, the intrinsic H α /H β Balmer decrement is generally estimated to be in the range of 3.0–3.1. This is slightly higher than the recombination Case B value (e.g., Gaskell and Ferland, 1984; Veilleux and Osterbrock, 1987). Therefore, for HII regions we adopt intrinsic Balmer decrement value of 2.86 as expected for standard Case B recombination at T ≈ 10,000   K and n e ≈ 100   c m − 3 , while for other cases LINER, Composites and Seyfert galaxies, where collisional excitation, harder radiation, or dust effects may occur, the intrinsic value is assumed to be 3.1.

To quantify the extinction, we applied the reddening law, which expresses the total extinction k ( λ ) as a function of wavelength: k λ = a x + b x R V , (1) where x = 1 λ (in units of μ m − 1 ) is the inverse wavelength, while R V presents the ratio of total-to-selective extinction, and determines the shape of the reddening law, which varies with the physical and chemical properties of the dust along the line of sight. We varied the R V = ( 2.1 , 3.1 , 4.1 ) , to examine how different reddening assumptions affect the results for each sample (see, Table 1).

The ratio R V = 3.1 is the standard value for interstellar dust in the diffuse Milky Way ISM. This value is commonly used for star-forming regions as extinction law for diffuse clouds Cardelli et al. (1989). For Seyfert 2, composite, and LINERs R V may differ, reflecting an extinction law of different steepness. These other values are motivated by observations of AGN and similar environments, where the extinction curve often deviates from the standard Milky Way curve due to differences in the dust composition and density (Maiolino et al., 2001) stop.

The selection of different R V = ( 2.1 , 3.1 , 4.1 ) for AGN-related objects is supported by studies suggesting that the grain size distribution in AGN and LINER environments may differ from that in the Milky Way. Specifically, these regions may contain a larger fraction of small grains, or an increased presence of large grains, leading to variations in R V . Observational studies of reddening curves in AGN environments often indicate steeper extinction laws, providing justification for considering values of R V > 3.1 .

The functions a ( x ) and b ( x ) were computed using the series expansions for the optical wavelength range (see Cardelli et al., 1989)

The reddening correction was applied to the observed fluxes F ( λ ) of lines by using the relation: R corr = R obs × 1 0 − E B − V k λ , (2) where R is the ratio of emission lines and E ( B − V ) is the color excess determined from the Balmer decrement (see,for e.g., Shivaei et al., 2020): E B − V = 2.5 × log H α / H β obs H α / H β int × 1 k H β − k H α mag . (3) where: ( H α / H β )_obs is the observed flux ratio of H α to H β , ( H α / H β )_int is the intrinsic flux ratio of H α to H β , k ( H α ) and k ( H β ) are the extinction coefficients for H α and H β , respectively and E ( B − V ) is the color excess due to reddening.

The forbidden [O III] lines at 4,958.911 Å and 5,006.843 Å are emitted from transitions within the O 2 + ion, corresponding to the transitions 2 s 2 2 p 2 1 D 2 → 2 s 2 2 p 2 3 P 1 and 2 s 2 2 p 2 1 D 2 → 2 s 2 2 p 2 3 P 2 respectively, and their ratio under the assumption of a low-density environment is theoretically expected to be 2.89 (Galavis et al., 1997), or 2.98 by including relativistic corrections to the magnetic dipole operator (Storey and Zeippen, 2000).

Observationally, some deviations have been noted: Iye et al. (1987) obtained 3.17 ± 0.04, and Leisy and Dennefeld (1996) reported 3.00 ± 0.08 (Popović et al. (2005)), 2.92 ± 0.08, Dimitrijević et al. (2007) 2.993 ± 0.014 and Laker et al. (2022) between 3.017 and 3.022.

We present our results of [OIII] λ λ 4,959, 5,007 flux ratios for every analyzed sub-sample of NELG in the top panel of Figures 2,3, and Table 1. Table 1 reports the sample size, the [OIII] λ 5,007/[OIII] λ 4,959 and [NII] λ 6,583/[NII] λ 6,548 intensity ratio averages and their associated errors, for three extinction curves with values or R V = 2.1 , 3.1 , 4.1 . The sample averages were computed applying a sigma-clipping algorithm iteratively, continuing the process until the values converged. The errors are the error on the average (i.e., the sample standard deviation divided by the number of degrees of freedom, sample size N − 1 ).

Figure 2 presents the comparison between our measurements of the average values calculated for three extinction laws ( R V = 2.1 , 3.1 , 4.1 ) and the theoretical values of Storey and Zeippen (2000) and Galavis et al. (1997). The statistical significance of the results was estimated by computing a z estimator defined as z = | < [ O I I I ] λ 5007 / [ O I I I ] λ 5007 obs > − [ O I I I ] λ 5007 / [ O I I I ] λ 5007 theor | / σ mean which follow a Student’s t distribution. We find that [OIII] line ratios correspond to theoretical prediction of Storey and Zeippen (2000) for samples Sy, LINER and composite for any of the considered values of R V , while for HII regions this solution correspond to the cases 3.1 ≤ R V ≤ 4.1 . The case R V = 4.1 is marginally different from 2.98 (probability P ∼ 1 − 6.1 0 − 3 ), but definitely not consistent with 2.89.

In the Figure 3 we present the line ratio distributions for each sample (Seyfert 2, LINER, composite, and H II regions), assuming R V = 3.1 for reddening corrections. The left panel displays the results for [O III], while the right panel shows the results for [N II].

We applied the Shapiro-Wilk test to assess Gaussianity of the distributions presented in the Figure 3. The test confirmed Gaussianity only for the distribution of line ratios in the composite sample, across all tested combinations: both line ratio distributions and all R V values. For all other samples, the distributions deviate from a Gaussian profile. In our plots, the Gaussian overlay is included primarily for illustrative purposes, serving as a reference rather than an indication of an exact fit to the data.

The Figure 3 present the distribution of [OIII] λ λ 4959, 5007 flux ratios in a form of histograms.

The [N II] emission lines at 6,548 Å and 6,583 Å arise from transitions within the N + ion. Two transitions are commonly observed in astrophysical spectra: (1) 2 s 2 2 p 2 1 D 2 → 2 s 2 2 p 2 3 P 1 ([N II] λ 6548.05 Å) and (2) 2 s 2 2 p 2 1 D 2 → 2 s 2 2 p 2 3 P 2 ([N II] λ 6583.45 Å). These radiative transitions are forbidden in electric dipole radiation due to parity violations and inter-combination restrictions ( Δ S = 0 ) . However, the [N II] λ 6548 and λ 6583 lines are often strong in optical spectra from low-density, photoionized regions due to the vast emitting volumes.

The critical density for collisional de-excitation of the ¹D₂ level is 6.6 × 1 0 4 cm − 3 (Osterbrock and Ferland, 2006). In environments with higher electron densities, this metastable level is more likely depopulated by electron impacts rather than by radiative transitions, causing [N II] λ λ 6,548, 6,583 Å lines to be absent in spectra. These lines are often observed adjacent to the H α line, making them useful for diagnostics in star-forming regions as well as in AGN environments.

The theoretical intensity ratio of [N II] λ 6583 to [N II] λ 6548 is 2.94 (e.g., Condon, 1934; Galavis et al., 1997) or 3.05 with including relativistic corrections to the magnetic dipole operator (Storey and Zeippen, 2000).

Figure 3 shows the distribution of [NII] λ λ 6548,6583 flux ratios for four subclasses of NELG spectra. Results of sigma-clipped average value of the distributions for each type of galaxy are presented in Table 1.

The result for starburst spectra are in a good agreement with theoretical expectation of Storey and Zeippen (2000), especially for the case of R V = 4.1 . Results for LINERs, composites and Sy2 are significantly lower and hence closer to theoretical values of Galavis et al. (1997), for all versions of analyzed extinction laws. The z estimators confirms a significant difference with Storey and Zeippen (2000) at a confidence level ≳ 3 σ in all extinction cases and for all classes save the HII nuclei. However, due to the large dispersion, Sy2 and LINER classes are only marginally in disagreement with [Storey and Zeippen (2000)]. The composite class, due to the large sample size, is in significant disagreement with both 2.94 and 3.05 (see Figure 2).

The reason could be physical and environmental conditions that can result with slight changes in [NII] line ratio as mentioned above, for example, if the density fluctuations reach values close to critical density, due to shocks (see, Rich et al., 2011). The [NII] lines are more sensitive to physical conditions then [OIII]. The highest discrepancy (3.6%) between our result and theoretical values is in the case of LINERs, where shocks are typically expected. We can conclude that our measurements are in overall agreement with theoretical predictions.

In addition to the [OIII] line ratio analysis, we also investigated the [NII] λ 6584 / λ 6548 ratios across different galaxy classes and extinction laws. Unlike the case of [OIII], where all subsamples agree closely with the theoretical value of 2.98 (given by Storey and Zeippen, 2000), the results on the [NII] ratios are less straightforward to assess. Our measured [NII] ratios vary depending on both galaxy type and the adopted extinction law, spanning the interval between two commonly cited theoretical predictions in the literature: 2.94 (Galavis et al., 1997) and 3.05 (Storey and Zeippen, 2000).

To quantify the agreement between our data and theoretical expectations, we performed a two-tailed z -test for each galaxy class and extinction curve ( R V = 2.1 , 3.1, 4.1), using both theoretical values as references. The results reveal that the H ii galaxies are in best agreement with the Storey and Zeippen (2000) prediction, particularly under R V = 4.1 where p = 0.74 . In contrast, the LINERs, Seyfert 2s, and Composite galaxies show statistically significant disagreement with the Storey and Zeippen (2000), while being more consistent with the Galavis et al. (1997) value.

These deviations are plausible. The [NII] lines are more sensitive to gas density and ionization structure than [OIII], and are particularly sensitive to shock excitation and other processes in AGN-driven environments. The largest deviation is found in LINERs, where shocks are expected to play a significant role (Rich et al., 2011). Overall, the measured [NII] ratios are consistent with theoretical predictions when the environmental and excitation conditions specific to each galaxy type are considered.