B2 0003+38A: A Classical Flat-spectrum Radio Quasar Hosted by a Rotation-dominated Galaxy with a Peculiar Massive Outflow

B2 0003+38A is classified as an FSRQ due to its BELs in the optical spectrum and the flat radio spectrum with α = −0.3 (Stickel & Kuehr 1994; Healey et al. 2007; Massaro et al. 2009). Resolved by very long baseline interferometry (VLBI) mapping at 2.29 GHz (Morabito et al. 1982), it shows a prominent core that emits more than 95% of the 2.3 GHz flux, along with a weak jet pointing to the southeast (Fey & Charlot 2000), as shown in Figure 1. Kinematic analyses show quasi-stationary knots at the jet base and relativistic motions downstream (Hervet et al. 2016). Aditya & Kanekar (2018) detected H i 21 cm absorption toward this source, finding the velocity-integrated H i 21 cm optical depths to be 1.943 ± 0.057 km s⁻¹, and the H i column density to be 3.54 ± 0.11 × 10²⁰ cm⁻².

Zoom In Zoom Out Reset image size

The quasar B2 0003+38A was observed by Keck II ESI long-slit spectrography in its echelle mode on 2015 September 9. The spectrum was taken in the optical band (wavelength coverage λ ∼ 3995–10198 Å in the observer's frame), covering 3250–8297 Å in the rest frame for our target at z ∼ 0.2. A 125 wide slit was employed, resulting in an instrumental dispersion σ of 22 km s⁻¹. The position angle of the slit is 41° (north to east, Figure 1). The slit width allows for both the blazar and its jet to be covered. Two exposures were taken with an integration time of 10 minutes each. Arc lamps were used during the observing campaign for wavelength calibration, and the standard star BD+28 4211 was observed about one hour before the target with identical settings for the purpose of flux calibration.

We use a customized routine based on the IDL package XIDL ⁵ to obtain a two-dimensional (hereafter 2D) spectrum. Our data reduction consists of four steps: (1) after bias corrections, flat-fielding, and removing cosmic rays, we stack the two exposures; (2) the nucleus of the target quasar is traced along 10 echelle orders, and for each order, we resample the 2D spectrum so that a resultant pixel corresponds to 10 km s⁻¹ in the wavelength direction and 015 in the spatial direction; (3) sky subtraction and flux calibration are applied on these resampled 2D spectra; (4) we combine these resampled 2D spectra from the 10 orders.

We extracted the 1D spectrum from the 2D spectrum as a result of data reduction by applying an aperture with a size of 3'' × 125 centered on the quasar nucleus, as shown in Figure 2 (black lines). AGN features, including BELs, narrow emission lines (NELs) and Fe ii bumps, are seen therein. However, the existence of high-order Balmer absorption lines indicate the continuum contribution from starlight. To disentangle different contributions of the spectrum, we construct a model consisting of four components: a stellar, a power-law, a BEL, and an Fe ii component.

Zoom In Zoom Out Reset image size

• 1.  
  The stellar component. We assume a single simple stellar population (SSP), utilizing the library of Bruzual & Charlot (2003). When fit to the data, the template is shifted to a redshift z_⋆, broadened by convolving to a Gaussian parameterized by a velocity dispersion parameter σ_⋆, and multiplied by the extinction curve of the Small Magellanic Cloud (SMC) parameterized by E(B − V)_⋆ (Pei 1992).
• 2.  
  The power-law component. We assume a formulation of f_λ (λ) = C λ^α , where α is the exponent and C is the normalization.
• 3.  
  The BEL component. For this component, we mainly account for the Balmer series, He i λ5876 and He ii λ4686. We assume that the BEL profiles can be represented by a linear combination of two Gaussians, whose mean and standard dispersion are fixed but the flux is allowed to vary. Furthermore, we fix the flux ratios of higher-order Balmer BELs to the Hγ BEL to those under the "Case B" situation (Storey & Hummer 1995), assuming a typical circumstance for broad-line regions with T_e = 15,000 K and n_e = 10⁹ cm⁻³.
• 4.  
  The Fe ii component. We employ the Fe ii template constructed by Véron-Cetty et al. (2004). After experiments with various combinations of different Fe ii lines, we conclude that only narrow Fe ii lines are necessary for fitting the observed Fe ii complex. We further assume that these different narrow Fe ii lines can be represented by single Gaussians with a fixed redshift and a fixed width. The flux ratios between different Fe ii lines are taken from the template given in Véron-Cetty et al. (2004). In addition, we assume the SMC extinction curve for the dust attenuation and reddening of the Fe ii component parameterized by E(B − V)_{Fe II} .

We then fit the spectrum of the nuclear region, where the contamination from NELs is minimal, resulting in a minimized reduced chi-square of 1.99. This fit can be further improved if the regions affected by telluric absorptions are dismissed (reduced χ² = 1.66, Figure 2). The individual components of the best-fit spectra are depicted by colored lines in the same figure, and the best-fit parameters are tabulated in Table 1.

3.1.1. A Narrow Absorption Line System

The nuclear spectrum of B2 0003+38A shows deep and narrow Na i λλ5890,5896 absorption lines (Figure 3), corresponding to a redshift of ∼0.22883, which is about 70 km s⁻¹ blueshifted relative to the stellar redshift (0.22911). Unlikely to originate from the stellar populations, these are likely absorption lines from absorbers that happen to lie in our line of sight toward the quasar nucleus. Narrow Ca ii λ3934 absorption is detected at the same redshift as that of Na i absorption, implying the same absorption system, though we cannot affirm the existence of the Ca ii λ3969 absorption line due to the influence of the [Ne iii] λ3869 and H emission lines. Assuming that the starlight is not affected by the absorber mentioned above, we consider the case of partial coverage:

$$\begin{eqnarray}&&f(\lambda )={f}_{\mathrm{qso}}(\lambda )\times (1-{C}_{f}+{C}_{f}\times {e}^{-\tau (\lambda )})+{f}_{\mathrm{stellar}},\end{eqnarray} \tag{ 1 }$$

where f_qso is the flux density of the emission from the quasar nucleus (including the power-law continuum, the BEL, and the Fe ii components), f_stellar is the flux density of stellar emission, C_f is the covering factor, and τ(λ) is the optical depth of absorption. The τ(λ) profile that we adopt for each absorption line is Gaussian with a fixed velocity dispersion. We perform this nuclear spectrum fitting in the vicinity of the narrow absorption lines, and the resultant best-fit model and the corresponding parameters are given in Figure 3 and Table 2, respectively. We find the redshift of these absorption lines to be z = 0.228827 ± 0.000002, highly consistent with that of the H i 21 cm absorption (z ∼ 0.2288; Aditya & Kanekar 2018). In addition, the H i 21 cm absorption line is spectrally unresolved, in line with the narrowness of the Na i and Ca ii absorption lines (σ ∼ 23 km s⁻¹). These consistencies strongly imply a relation between the absorption of Na i/Ca ii and that of H i. Aditya & Kanekar (2018) measured the H i column density to be (3.54 ± 0.11) × (T_s /100 K) × 10²⁰ cm⁻², where T_s is the spin temperature. The H i column density N(H) and the Na i column density N(Na I) are related by the following equation (e.g., Rupke et al. 2005):

$$\begin{eqnarray}&&N({\rm{H}})=N(\mathrm{Na}\,{\rm\small{I}}){\left(1-y\right)}^{-1}{10}^{-(a+b)},\end{eqnarray} \tag{ 2 }$$

where y is the ionization fraction, a is the Na abundance, and b describes the depletion onto dust. If we adopt y = 0.9 following Rupke et al. (2005), a solar abundance (so that a = −5.69), and the canonical Galactic depletion value of b = − 0.95, then we see that the measured N(Na I) corresponds to N(H) ∼ 1.1 × 10²¹ cm⁻², a value close to the result of radio spectral analysis.

Zoom In Zoom Out Reset image size

Table 2. Parameters for Modeling the Narrow Absorption Lines

Parameter	Value	Unit
z	0.228827 ± 0.000002
σ	23.4 ± 0.6	km s−1
Cf	0.78
NNa I	26 ± 2	1012 cm−2
NCa II	3.4 ± 0.7	1012 cm−2

Download table as:  ASCIITypeset image

BELs in the optical spectrum are unambiguously emitted from the quasar nucleus. We measure the Balmer decrement of the BELs to be 8.07 ± 0.05, significantly higher than the theoretical "Case" B value of 2.7 (Gaskell 2017), indicating heavy dust reddening toward the quasar nucleus. Assuming an SMC extinction curve results in an estimation that E_B−V ∼ 1.2. The power-law component contributes 60%–70% of the total continuum flux in the wavelength range of 4000–7000 Å. The quasar continuum is remarkably red with a power-law index α of 1.32 ± 0.02, probably a result of the synchrotron emission from the radio jet, and/or the reddened thermal emission from the quasar nucleus.

In view of the heavy dust extinction in optical bands, we use the infrared continuum to estimate the bolometric luminosity of the quasar. As the first step, we obtain the 5 μm monochromatic luminosity using the infrared SED constructed from the AllWISE photometry of the quasar, finding that ν L_ν (5 μm) = 8 × 10⁴⁴ erg s⁻¹. We then follow Richards et al. (2006) to apply a bolometric correction factor of 8, reaching a bolometric luminosity of 6 × 10⁴⁵ erg s⁻¹, though this value may have been overestimated due to the possible contribution of synchrotron emission from the jet at 5 μm.

To investigate the gas kinematics of the host galaxy through the 2D profiles of NELs, we remove the emission from the quasar nucleus and the stars, each of which is represented by double Gaussians in the spatial regions where the NELs are negligible. The two top panels of Figure 4 show the rest-frame 2D spectra in the neighborhood of the [O i], [S ii] λ λ 6717, 6731 doublet, the Hα, N ii doublet, and the Hβ, [O iii] λ λ 5007, 4969 doublet (note that we employ the stellar redshift, z = 0.22911). The NEL structure extends to a size of ∼22, greater than the FWHM of the point-spread function (PSF) 071 (see Appendix A), and thus is spatially resolved.

Zoom In Zoom Out Reset image size

These 2D NEL spectra reveal three components. In particular, the 2D spectra of [O i]- and [S ii]-emitting NELs show a velocity gradient across the spatial extent, indicating a rotation-dominated disk. The velocity gradient can be seen in the 2D NEL spectra produced using both the two independent exposures (see Appendix B), and thus is reliable. The 2D spectra of Hα and Hβ NELs (Figure 4) imply the existence of an extended region located ∼4–20 kpc southwest of the quasar nucleus, whose position and velocity is reminiscent of EELRs around quasars (e.g., Stockton & MacKenty 1987; Fu & Stockton 2009). The 2D spectra of the [O iii] NEL is more complicated: besides the two components mentioned above, an additional blueshifted component centered on the quasar nucleus is evident. Such a profile is frequently seen in the 2D spectra of quasars' [O iii] emission and is conventionally considered to be suggestive of the existence of outflowing gas. Therefore, the three components of the 2D NEL spectra include a rotating disk, an EELR, and an outflow. In addition, we see double-peak profiles in [O iii] Hα and Hβ lines in the southwest of the nucleus, indicating multiple components.

To delineate the gas kinematics, we perform a two-step spectral fit to decompose these three components using the Python package MPFIT. In the first step, for each spatial element (spaxel), we fit a double Gaussian to the profile of Hβ, [O iii], Hα, and N ii. The resultant best-fit Gaussian models for a single spaxel are demonstrated in Figure 5. Our reconstructed 2D spectra of the red Gaussian show the EELR only, while the blue Gaussian is a superposition of a rotating disk and an outflow (Figure 4, the third and fourth rows). In the second step, for each individual spaxel, we perform a single-Gaussian fit to the [O i] emission line and the [S ii] λ λ 6717, 6731 doublet. The best-fit spectra of [O i] and the [S ii] doublet in seven spaxels are shown in Figure 6 (left panels), from northeast to southwest. Due to the higher complication of their profiles, we fit the Hα and Hβ profiles to a single or double Gaussian, for which the decision is made in a way similar to that in Liu et al. (2014). It turns out that the five spectra closest to the nucleus demand a double Gaussian. These fits allow us to successfully decompose the contribution from the rotating disk and a nearly spherical outflow structure, and the best-fit Hα and Hβ spectra of seven spaxels are shown in Figure 6 (right panels), from northeast to southwest.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We scrutinize the rotating disk by measuring the velocity relative to the redshift of host galaxy and the velocity dispersions of best-fit spectra of the rotation components in the NELs of [O i], [S ii], Hβ, Hα and N ii. In Figure 7, we plot them as a function of the distance from the quasar nucleus (positive is to the southwest and negative to the northeast). In addition, we only show those spaxels with a high signal-to-noise ratio (S/N of [O i] ≳ 3), corresponding to ∼5 kpc from the nucleus. We see that the overall velocity distribution and dispersion of these NELs are in line with typical galactic disks (Rix et al. 1992; Proshina et al. 2020). We find that the velocity profiles measured from different NELs are consistent, and thus use the median of the measurement results from these different NELs, following the method employed by Courteau (1997), Weiner et al. (2006), and Drew et al. (2018). In detail, assuming an axis ratio of b/a = 1, the distribution of velocity is formulated as:

$$\begin{eqnarray}\begin{array}{rcl}{V}_{\mathrm{obs}}(r) & = & \left(\displaystyle \frac{2}{\pi }{V}_{{\rm{a}}}\arctan \left(\displaystyle \frac{{r}_{\mathrm{obs}}}{{r}_{{\rm{t}}}\cos (i)}\right)+c\displaystyle \frac{{r}_{\mathrm{obs}}}{\cos (i)}\right)\\ & & \times \sin (i)+{V}_{\mathrm{rel}},\end{array}\end{eqnarray} \tag{ 3 }$$

where V_a is the asymptotic velocity, r_t is the knee radius, r_obs is the projection distance from the center of the galaxy, c is the outer-galaxy slope, V_rel is the velocity of the ionized gas relative to the stars along the line of sight, and i is the inclination angle of the galaxy. Here V_a and r_t are dimensional scaling parameters, whereas c characterizes the shape of the rotation curve. Meanwhile, if a Gaussian fit is used, the profile of the velocity dispersion is given by:

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)=A\exp \left(-\displaystyle \frac{{\left(r-m\right)}^{2}}{2{\omega }^{2}}+{\sigma }_{0}\right),\end{eqnarray} \tag{ 4 }$$

where Σ(r) is the value of the fitted Gaussian at each radius, m is the central position (if the ionized gas is isotropically distributed around the center, then m is 0), ω is the width of the Gaussian, and σ₀ is the isotroptic component of the velocity dispersion. Here we correct the observed velocity dispersion for the intrinsic instrumental dispersion using $\sigma =\sqrt{{\sigma }_{\mathrm{obs}}^{2}-{\sigma }_{\mathrm{inst}}^{2}}$, where σ is reported in Figure 7, and σ_inst is the combined instrumental- and spectral-seeing dispersion that we measure to be ∼23 km s⁻¹. The resultant best-fit velocity dispersion σ is 126 ± 10 km s⁻¹.

Zoom In Zoom Out Reset image size

We fit the velocity and velocity dispersion profiles to the median values of the measurements results from different NELs ([O i], [S ii], Hβ, Hα, and N ii), utilizing a Monte Carlo (MC) simulation that we run 50 times to calculate the errors of parameters. Errors on the median are calculated as $\sigma ={\sigma }_{\mathrm{NMAD}}/\sqrt{n-1}$, where σ_NMAD is the normalized median of the absolute deviations (Hoaglin et al. 1983) and n is the number of samples, here n = 5. The best-fit profiles are shown in Figure 7, along with the corresponding parameters listed in Table 3. The velocity profile reveals a roughly symmetric gas motion pattern. Gas in the quasar's nuclear region has a velocity of ∼−75 km s⁻¹, and is blue-/redshifted on the southwest/northeast side, respectively. The velocity profile flattens beyond a radius of 4 kpc in both of these two directions, with a velocity of −150 to −180 km s⁻¹ (southwest) and ∼40 km s⁻¹ (northeast) within the distance range of 4–6 kpc. The rotating component is blueshifted relative to the systemic (stellar) velocity with V_rel = − 71 km s⁻¹. Since the narrow absorption line system is blueshifted with a similar relative velocity (−85 km s⁻¹; see Section 3.1.1), this kinetic consistency is suggestive of the same origin of the two.

In the central region, we note that a velocity difference of ∼50 km s⁻¹ exists at a > 3σ significant level between [S ii] and other NELs. This is likely due to the outflow, though the possibility of influence from the sky emission line in the blue outskirts of [S ii] cannot be fully ruled out, but if [S ii] is excluded from these analyses, our results change minimally.

The EELR extends to the southwest of the nucleus (see Figure 4), which is most evident in the [O iii] emission, where two compact knots and a diffuse region are seen. For the line spectra at each position, we measure the integrated flux, median velocity, and line width (W₈₀, as detailed in Liu et al. 2013) of the Hα, Hβ, and [O iii] emission lines by fitting their profiles to Gaussian models, and conclude that the three NELs depict consistent velocity profiles (Figure 8).

Zoom In Zoom Out Reset image size

The integrated flux profile of all three NELs unambiguously show a series of three peaks at galactocentric distances of about 7, 13, and 18 kpc. Hence, we divide the EELR into three annuli, with radius ranges of 2–8, 8–16, and 16–25 kpc, respectively, and for each of which we extract a 1D emission-line spectra (Figure 9). The best-fit integrated fluxes, median velocities, and line widths are listed in Table 4, where it can be seen that the emission lines from 2–16 kpc are relatively narrow and are redshifted relative to the stellar redshift, while at 16–25 kpc they are broader and blueshifted.

Zoom In Zoom Out Reset image size

The intensity ratios of the emission lines facilitate our analysis on the physical conditions of the EELR. Specifically, we use the [O iii]/Hβ and Hα/Hβ ratios to quantify the degree of ionization and the dust attenuation, respectively. To obtain higher S/N ratios, we bin the spatial pixels within 0.45'' before measuring the median values of these line ratios, for which the uncertainties are, again, calculated using the MC method (Figure 10). We find that [O iii]/Hβ is larger than 10 at all galactocentric radii under our consideration (top panel therein), implying a high-ionization state in general. Meanwhile, considering the theoretical ratio Hα/Hβ ∼ 2.86 for "case B" (n_e = 100 cm⁻³, T_e = 10,000 K; Storey & Hummer 1995), we plot the result in Figure 10 (blue dashed line, bottom panel). The Hα/Hβ ratio at galactocentric radii of 2–16 kpc are roughly constant, corresponding to A_V ∼ 1. However beyond 16 kpc, Hβ is too weak for accurate measurement of dust attenuation.

Zoom In Zoom Out Reset image size

The electron density is achievable from [O ii] I(3729)/I(3726) and [S ii] I(6717)/I(6731) ratios. For this purpose, we stack the spectra taken from locations 2–16 kpc away from the center, and fit the [O ii] and [S ii] doublet emission lines by fixing the kinematics of the two lines (Figure 11). Assuming an electron temperature of 10,000 K and at a 68% confidence level, we estimate the electron density to be n_e = 137${}_{-30}^{+36}$ cm⁻³. This result is generally lower than that of typical narrow-line regions (n_e ∼ 1000 cm⁻³; Greene et al. 2011), but close to that of the EELR (Fu & Stockton 2009).

Zoom In Zoom Out Reset image size

As we mentioned in the Introduction, there are abundant works studying the host galaxies of blazars, particularly to distinguish whether the host galaxies are disk or elliptical galaxies. Most of them focused on imaging BL Lac objects, while none of them successfully found a disk-galaxy-hosted FSRQ. Because blazars are rare and can only be found in the distant universe, it is hard to spatially decompose the nuclei and the host galaxies well, even when using the Hubble Space Telescope's highest spatial resolution in the optical and NIR bands. However, we could use a spatially resolved spectrum to constrain the host galaxy properties. In disk galaxies, young stars and interstellar gas and dust rotate in disks around bulging nuclei, while in elliptical galaxies, old stars randomly swarm and gas and dust are lacking. Thus, the properties and kinematics of stars and gases can be used to distinguish between disk and elliptical galaxies.

In Section 3.4, we find that the kinematics of gas in the host galaxy of B2 0003+38A are dominated by rotation, and the curvatures of velocity and velocity dispersion are similar to those of disk galaxies (Ho & Martin 2020). In addition, there are other hints that support B2 0003+38A being hosted by a disk galaxy. The velocity of the ionized gas (rotating component) relative to the systemic (stellar) velocity with V_rel = −71 km s⁻¹ is similar to the velocity of the narrow absorption line system (−85 km s⁻¹), suggesting the same origin for these two. Furthermore, emission lines from host galaxy can be detected significantly, indicating a gas-rich host galaxy. The spectrum lacks the features of old stars, such as TiO molecular bands, indicating that old stars are not dominant. Meanwhile, the spectrum is red, suggesting that the starlight is dust reddened. Young stars and rich interstellar dust are also characteristic of disk galaxies. This is consistent with the analysis of a nuclear spectrum, which prefers a young stellar population (∼450 Myr) with large dust reddening (E(B − V)_⋆ = 0.62), under the assumption of a single SSP and quasar spectrum. Therefore, we conclude that B2 0003+38A is mostly likely hosted by a disk galaxy.

We find the EELR to the southwest of the nucleus, and it extends to a projected distance of up to 25 kpc. EELRs with such sizes were found around a substantial fraction of radio-loud quasars.

Before discussing the origin of the EELR further, we estimate some parameters of the EELR. The gases in parts one and two in Figure 9 have redshifted velocities of v₀ ∼ 120 km s⁻¹, and their distance to the galaxy nucleus is R₀ ∼ 15 kpc (see Figure 9). We estimate a dynamical timescale to be t ∼ 1.2 × 10⁸ yr, i.e., the time taken by the gas from the nucleus to reach such a distance with an average velocity of v₀. The mass of gas in the two parts can be estimated using Hβ luminosity L_Hβ and electron density n_e (e.g Liu et al. 2013; Harrison et al. 2014). After correction for dust attenuation, the summed luminosity in the two parts is ∼2.4 × 10⁴⁰ erg s⁻¹. The total mass of this ionized extended gas can be estimated as:

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{\mathrm{gas}}}{2.82\times {10}^{9}\,{M}_{\odot }}=\left(\displaystyle \frac{{L}_{{\rm{H}}\beta }}{{10}^{43}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}\right)\left(\displaystyle \frac{{n}_{e}}{100\,{\mathrm{cm}}^{-3}}\right).\end{eqnarray} \tag{ 5 }$$

We measure an electron density of n_e ∼ 137 cm⁻³ using [O ii] and [S ii] doublets, if the electron temperature is 10,000 K. We find M_gas ∼ 9.3 × 10⁶ M_⊙. Combined with the average velocity of 120 km s⁻¹, the total kinetic energy of the gas can be estimated as:

$$\begin{eqnarray}&&{E}_{\mathrm{kin}}=\displaystyle \frac{1}{2}{M}_{\mathrm{gas}}{v}_{\mathrm{gas}}^{2}=1.33\times {10}^{54}\,\mathrm{erg}.\end{eqnarray} \tag{ 6 }$$

If assuming the lifetime of the structure is the dynamical timescale, we can also estimate a mass rate $\dot{M}$ is 7.8 × 10⁻² M_⊙ yr⁻¹, and a kinetic energy rate ${\dot{E}}_{\mathrm{ink}}\sim 3.5\times {10}^{38}$ erg s⁻¹.

However, this electron density value may not be the bulk density of the EELR gas. The EELR might be illuminated by the quasar. If so, the ionization parameter U of the EELR can be estimated as:

$$\begin{eqnarray}&&U=\displaystyle \frac{{Q}_{{\rm{H}}}}{4\pi {r}^{2}{{cn}}_{{\rm{H}}}},\end{eqnarray} \tag{ 7 }$$

where Q_H is the rate of the hydrogen ionization photons from the quasar, r is the distance from the quasar nucleus to the EELR, and n_H is the hydrogen number density. Assuming that the intrinsic SED of the quasar is that given in Mathews & Ferland (1987), we estimate a Q_H of 1.2 × 10⁵⁶ s⁻¹ using the bolometric luminosity estimated in Section 3.2. Assuming a distance of 10 kpc, and assuming n_e = 1.2n_H for highly ionized plasma, the inferred U is about 10^−2.5. However, the observed [O iii]/Hβ values of 12–16 indicate a higher U value of ∼10⁻¹. This may be because the [O iii] emitting gas has a lower density than the previously estimated value. Using the mixed-medium model (Robinson et al. 2000), Stockton et al. (2002) divided the EELR of 4C 37.43 into two components, and found one component with a density of several hundred cm⁻³, and another with a density near 1 cm⁻³. Most of the [O iii], Hβ, and Hα come from regions with low density; [O ii] and [S ii] come from the regions with a density several hundred times higher. If this is the case for the EELR in B2 0003+38A, the actual value of M_gas is two orders of magnitude higher.

There are various possibilities for the origin of an EELR. In the case of an isolated galaxy, the gas can be an inflow, i.e., a cold accretion flow from the intergalactic medium, an outflow triggered by an AGN or a starburst, or the recycling of gas (of the outflow). In the case of a galaxy with interactions, the gas can also belong to tidal features. First, inflows are generally isotropic, which excludes the possible of the EELR being inflow driven. Second, EELRs are common in radio quasars, and these EELRs are generally related with outflows driven by radio jets (e.g., Fu & Stockton 2009), so the EELR in B2 0003+38A may also be the same. Third, verifying the possibility of the recycling of gas requires a panoramic image of the circumgalactic gas. The long-slit spectrum is not sufficient, and future integral field unit (IFU) observation data are needed. We did not see a companion galaxy around B2 0003+38A, and neither did we see any asymmetry from the brightness profile of the starlight. As for the possibility of galaxy interactions, we do not find any direct evidence of them, including any signature of asymmetry from the brightness profile of the starlight, or any companion galaxy around it. However, we note that the current data quality is not enough to rule out the possibility. The EELR may also originate in tidal features. If so, starlight accompanied with the EELR should be seen. This cannot be tested using the existing data. Higher-quality future observations, such as IFU, are required to fully investigate the origin of EELRs.