Deuterium fractionation of the starless core L 1498

We conducted radio and sub-millimeter (submm) observations toward L 1498 using the Institut de Radio-Astronomie Millimétrique (IRAM) 30-m telescope, the James Clerk Maxwell Telescope (JCMT), and the Green Bank Telescope (GBT). Observational parameters are summarized in Table 1 and the pointings of spectral line observations are shown in Fig. 1.

The continuum maps of L 1498 at near-infrared (NIR), mid-infrared (MIR), and sub-millimeter (submm) wavelengths are shown in Fig. 2. With the benefit of the deep NIR and MIR observations, the cloudshine phenomenon (Foster & Goodman, 2006) is detected at J, H, and $K_{s}$ bands, while the coreshine phenomenon (Pagani et al., 2010; Steinacker et al., 2010; Lefèvre et al., 2014) is detected at IRAC1 and IRAC2 bands. This cloud/coreshine detection indicates the presence of dust grain growth (Steinacker et al., 2015). We analyze the dust extinction in the above deep J, H, $K_{s}$, and IRAC2 images to derive the visual extinction map with a beam size of 50^′′ shown in Fig. 2g. We attempt to derive the total column density all over the cloud using the $A_{\rm V}$ map with the benefit that the dust extinction is dependent on the dust density but independent of the dust temperature (Pagani et al., 2004, 2015; Lefèvre et al., 2016). However, because of the high extinction toward the core center, the lack of sufficient $K_{s}$-band stars prevents us from deriving the actual $A_{\rm V}$ peak value but a lower limit of $A_{\rm V}$ at 25 mag toward the core center. This $A_{\rm V}$ map is then provided as a constraint on the outer density profile derived by the N₂H⁺ line emission (see Sect. 4.4).

The 850 $\mu$m dust continuum emission observed by JCMT (Fig. 2h) is optically thin and sensitive to the cold dust in the core. With a much smaller 850 $\mu$m beam size of 14^′′ compared to that of the $A_{\rm V}$ map, it clearly reveals the northwest–southeast-elongated core shape with a concave edge at the southern side. Compared with the N₂H⁺ (1–0) integrated line intensity map (Fig. 2i) with a beam size of 26^′′ obtained from Tafalla et al. (2004), we can see that both the 850 $\mu$m continuum emission and the N₂H⁺ (1–0) emission peak at the same position and have similar emission distributions, which implies that both trace the same region. Therefore, we use the 850 $\mu$m continuum map and the N₂H⁺ (1–0) integrated intensity map to determine the core center of L 1498 to be (RA, Dec)_J2000 = (4^h10^m52$\aas@@fstack{s}$97, +25°10^′ 18$\aas@@fstack{\prime\prime}$0). We can see that the chosen core center also coincides with the absorption center on the IRAC4 image (Fig. 2f), where the absorption feature is associated with the inner region of the core (Lefèvre et al., 2016). We denote the core center as a cross shown in each panel of Fig. 2.

Figure 1 shows the pointings of our observations. Our IRAM 30-m and GBT pointing observations are performed along a northeast–southwest cut across the core center, and one row in the JCMT pointings is overlaid on the same northeast–southwest cut. We used this cut (hereafter, the main cut) to analyze the physical and chemical properties of L 1498 (Sect. 4). Figure 3 shows the spectra of each emission line in the $T_{\rm A}^{*}$ scale along the main cut.

The L 1498 core is quiescent. Each hyperfine component of N₂H⁺ (1–0) and N₂D⁺ (1–0) is well separated. Along the main cut, the coverage of the N₂H⁺ (1–0), N₂H⁺ (3–2), DCO⁺ (1–0), and DCO⁺ (2–1) spectra span across the whole core. Including the N₂D⁺ and o-H₂D⁺ observations, we can see that the spectral intensities of these four cations peak toward the central region, suggesting that these cations are tracers for the core region. Meanwhile, N₂H⁺ (4–3) shows no detection beyond the 3$\sigma$ significance because of the low temperature inside the core and low central density compared to the critical density for this line ($\sim 2\times 10^{7}$ cm^-3 at 10 K). On the other hand, the o-H₂D⁺ (1₁₀–1₁₁) line is detected in a smaller elongated region with an extent of $\sim$30^′′ along the main cut and $\sim$60^′′ perpendicular to the main cut (also see Fig. 10). This indicates that o-H₂D⁺ traces the innermost region of the L 1498 core. We note that this is the first o-H₂D⁺ (1₁₀–1₁₁) detection toward L 1498. Caselli et al. (2008) conducted a survey of o-H₂D⁺ (1₁₀–1₁₁) in nearby dense cores with the CSO 10.4-m antenna. However, L 1498 is one of the two undetected targets out of ten starless cores in their sample, indicating that the deuterium fractionation in L 1498 is relatively low compared with the other starless cores of their study. Their non-detection is probably because their rms noise of $\sim$130 mK in the $T_{\rm MB}$ scale is larger than our rms noise of $\sim$50–85 mK in the same scale (see Table  1). In addition, we also have a few C¹⁸O (2–1) spectra observed along the main cut, of which the intensity shows a clear decrease in the central region although the minimum spectral intensity does not occur at the core center but at ($-$20^′′, $-$20^′′).

We aim to provide a 3D physical description of L 1498 and estimate the chemical timescale of the deuteration fractionation. To accomplish this goal, we utilized non-LTE radiative transfer modeling to analyze molecular line emissions, enabling us to evaluate the physical structure and molecular abundance profiles of the core. This abundance analysis is essential for estimating the chemical timescales of observed $X$(N₂D⁺)/$X$(N₂H⁺) ratios and other deuterated molecular abundances through time-dependent chemical modeling.

In the non-LTE radiative transfer modeling, we approximated the physical structure of L 1498 as two half onion-like models stuck together, comprised of multiple concentric homogeneous layers. The parameters in each layer are number density ($n_{\rm H_{2}}$), gas kinetic temperature ($T_{\rm kin}$), relative abundances with respect to H₂ of the observed species ($X_{\rm species}=n_{\rm species}/n_{\rm H_{2}}$), turbulent velocity ($V_{\rm turb}$⁶⁶6$V_{\rm turb}$ is defined as the 3D isotropic non-thermal velocity dispersion. If $T_{\rm kin}$ and $V_{\rm turb}$ are constant along the line of sight, the FWHM of the line profile is given by $\delta v_{\rm FWHM}^{2}=8\log(2)\cdot\left(\frac{k_{\rm B}T_{\rm kin}}{m}+\frac{1}{3}V_{\rm turb}^{2}\right)$, where $k_{\rm B}$ is the Boltzmann constant and $m$ is the mass of the observed species.), rotational velocity field ($V_{\rm rot}$), and radial velocity field ($V_{\rm rad}$). The layer width is chosen to be $10\sqrt{2}$ arcsec ($\approx 14.1$^′′ and 1980 au at the distance of 140 pc), the same as the diagonal spacing of our IRAM 30-m/GBT pointing observations along the main cut. We therefore can determine these parameters at each layer sequentially from the outermost to the innermost layer by sampling sightlines of progressively decreasing radius along the main cut across L 1498.

Despite observations indicating an infall motion within the L 1498 envelope (Lee et al., 2001; Lee & Myers, 2011; Magalhães et al., 2018), the L 1498 core remains quiescent, as evidenced by the single-peaked spectra observed in our study. Furthermore, the velocity gradient of N₂H⁺ within this region is characterized by small magnitudes and seemingly random orientations, with a reported mean magnitude of 0.5$\pm$0.1 km s^-1 pc^-1 and a direction of $9^{\circ}$$\pm$10^∘ (Caselli et al., 2002). Based on these observations, we have proceeded with the assumption that $V_{\rm rad}=0$ and $V_{\rm rot}=0$.

Then our procedure to determine these parameters at each layer is as follows. First, we used multi-transition N₂H⁺ spectra to determine the $n_{\rm H_{2}}$, $T_{\rm kin}$, and $X$(N₂H⁺) profiles, assuming a uniform $V_{\rm turb}$ value. Second, we determined the abundance profiles of N₂D⁺, DCO⁺, and o-H₂D⁺ by assuming that they share the same density, temperature, and kinematic properties of N₂H⁺. The details of the radiative transfer modeling are explained in Sect. 4.2. With the obtained abundance profiles of N₂H⁺, N₂D⁺, DCO⁺, and o-H₂D⁺, we employed a pseudo-time-dependent chemical model (Pagani et al., 2009b) to estimate the chemical timescales of each layer with our derived density and temperature profiles. In this adopted chemical model, the free parameters include the abundances of CO and N₂, the initial OPR of H₂ (OPR_intial(H₂)), the average grain radius ($a_{\rm gr}$), and the cosmic ray ionization rate ($\zeta$). These parameters are either determined or assumed with chosen values during the modeling process. Further details of the chemical modeling are explained in Sect. 4.3.

We basically followed the procedure in Lin et al. (2020), with the main difference being the construction of an asymmetric onion-like model for L 1498. Additionally, while Lin et al. (2020) independently determined $n_{\rm H_{2}}$ using the NIR/MIR dust extinction measurements, we simultaneously determined $n_{\rm H_{2}}$ along with $T_{\rm kin}$, $X$(N₂H⁺), and $V_{\rm turb}$ by reproducing our multi-transition N₂H⁺ spectra through radiative transfer modeling. This was due to the lack of sufficient $K_{s}$-band stars, which prevented us from deriving the actual $A_{\rm V}$ peak value but allowed us to establish a lower limit of $A_{\rm V}$ toward the core center (see Sect. 4.4). Thus our dust extinction map served primarily as a constraint on the $n_{\rm H_{2}}$ profile at outer radii, where $A_{\rm V}$ can be determined. We will further discuss the agreement of the dust extinction with our $n_{\rm H_{2}}$ in Sect. 5.1.

Figure 4a and 4b show two onion-like models of L 1498. One is 1D spherically dissymmetric and the other is 3D asymmetric (see Appendix D). With our spectral observations along the main cut, we assumed that both models comprise nine and five layers on the southwest and northeast sides, respectively. In order to reproduce our observed spectra, we adopted a 1D spherically symmetric non-LTE radiative transfer code (MC; Bernes, 1979; Pagani et al., 2007) for the 1D spherical case, and the Simulation Package for Astronomical Radiative Transfer/Xfer code (SPARX⁷⁷7https://sparx.tiara.sinica.edu.tw/) for the 3D asymmetric case. In both codes, we turned on the hyperfine-line-overlapping feature, and provided the updated hyperfine-line-resolved collisional rate coefficients (Lique et al., 2015; Pagani et al., 2012; Lin et al., 2020). We then obtained our modeled spectra calibrated in the $T_{\rm A}^{*}$ scale, whereas the $T_{\rm MB}$ scale can introduce an overcorrection for low main-beam efficiencies, $\eta_{\rm MB}$ (Bensch et al., 2001). Particularly in the 1.3–0.8 mm range observed with the IRAM 30-m telescope and also in the 4 mm range observed with the GBT, the error beams can pick up a substantial fraction of the signal of extended sources (see Sect. 2.1.2). This leads to overestimates of the temperature and/or density and also has an impact on the molecular abundance and the abundance ratio. Calibrating in the $T_{\rm A}^{*}$ scale by considering the telescope coupling contributions from the main beam and error beams allows for a better representation of extended sources without having to handle the corrections of the error beam pick-up in the $T_{\rm MB}$ scale (see Appendix A).

We search for the best-fit profiles of $n_{\rm H_{2}}$, $T_{\rm kin}$, abundances of our observed species, and $V_{\rm turb}$ with the 1D spherical model by our spectral fitting procedure. Initially, we ran MC to reproduce our N₂H⁺ (1–0, 3–2, and 4–3) spectra across the main cut. This allowed us to determine $n_{\rm H_{2}}$, $T_{\rm kin}$, $X$(N₂H⁺), and $V_{\rm turb}$ layer by layer, sequentially from the outermost to the innermost layer, for both the northeast and southwest sides. We found that a uniform $V_{\rm turb}$ of $0.09$ km s^-1 can reproduce the spectral line widths along the main cut, with all spectra consistent with a systematic velocity of 7.8 km s^-1. These initial $n_{\rm H_{2}}$, $T_{\rm kin}$, and $X$(N₂H⁺) profiles served as the initial guess for refinement. Subsequently, we used the simulated annealing (SA) algorithm provided by the Modeling and Analysis Generic Interface for eXternal numerical codes (MAGIX; Möller et al., 2013) to optimize the fit of the N₂H⁺ spectra, obtaining the best-fit profiles of $n_{\rm H_{2}}$, $T_{\rm kin}$, and $X$(N₂H⁺) from our initial profiles. The 1$\sigma$ uncertainty of each parameter at every layer was determined using the interval nested sampling error estimation method provided by MAGIX, with the other parameters fixed. Afterward, we determined the abundance profiles of N₂D⁺, DCO⁺, and o-H₂D⁺, along with their respective 1$\sigma$ uncertainties, by fitting the spectra of N₂D⁺ (1–0, and 2–1), DCO⁺ (1–0, 2–1, and 3–2), and o-H₂D⁺ (1₁₀–1₁₁). This fitting process assumed identical $n_{\rm H_{2}}$ and $T_{\rm kin}$ profiles, as well as a constant $V_{\rm turb}$, as for N₂H⁺. Because of our shorter spatial coverage in the N₂D⁺ observation, we limit $X$(N₂D⁺) in the outer layers with $X$(N₂D⁺)/$X$(N₂H⁺) lower than the observed minimum $X$(N₂D⁺)/$X$(N₂H⁺) of 0.07. The best-fit central $n_{\rm H_{2}}$ was found to be $1.6_{-0.3}^{+3.0}\times 10^{5}$ cm^-3, and the central $T_{\rm kin}$ is 7.5${}_{-0.5}^{+0.7}$ K.

Next, we transform the above best-fit profiles from the 1D spherical model to the 3D asymmetric model because carrying out a non-LTE calculation with a 3D asymmetric model repeatedly in the fitting is much more numerically expensive than the 1D case. Our 3D asymmetric model was designed to follow the emission morphologies of the 850 $\mu$m continuum and the N₂H⁺ (1–0) integrated intensity. We set the southwest sides in both models to have identical physical parameters and abundance profiles. For the northeast side, the spherical model was stretched along the sightline direction by a factor of $\sim$2 to match the thickness of the southwest side along the line of sight (see Appendix D). Keeping the same northeast profiles in the asymmetric model results in too strong intensities in the spectra because the molecular abundances are doubled. We aim to keep the continuity in the $n_{\rm H_{2}}$ and $T_{\rm kin}$ profiles across two sides, so we simply divided the abundances in the northeast sides by a common factor of $\sim$2 to make the modeled spectra fit with the observations. We note that our intention is to obtain a simultaneous fit that is generally consistent with the observational spectra since performing a complete search of the entire parameter space for the best possible fit is not achievable.

Figure 5 shows our best-fit profiles of $n_{\rm H_{2}}$, $T_{\rm kin}$, abundances of the above four cations, and the N₂H⁺ deuteration ratio, $X$(N₂D⁺)/$X$(N₂H⁺), along the main cut across the center of our asymmetric onion-like model, while Figure 4c shows the 3D visualization of the $n_{\rm H_{2}}$ distribution. We also fit our $n_{\rm H_{2}}$ profiles with the Plummer-like profile with $n_{\rm H_{2}}(r)=n_{0}/(1+(r/R_{0})^{\eta})$, where $n_{0}$, $R_{0}$, and $\eta$ are central density, flattening radius, and the power-law index, respectively. The fitted Plummer-like profiles are shown with green curves and the fitted parameters are annotated ($n_{0}=1.6\times 10^{5}$ cm^-3, $R_{0}=8.2$ kau and $\eta=2.23$ for the southwest profile, and $R_{0}=3.0$ kau and $\eta=2.30$ for the northeast profile) in the top row in Fig. 5. The above best-fit profiles are also numerically listed in Table 11. In addition, in Fig. 5, the $n_{\rm H_{2}}$, $T_{\rm kin}$, and $X$(N₂H⁺) profiles of the spherical model found by Tafalla et al. (2004) are also shown in purple. Since their model center was chosen at a different position with an offset of ($\Delta\mbox{RA}=-10$^′′, $\Delta\mbox{Dec}=-20$^′′) with respect to our model center, we compute these purple curves along our main cut (i.e., a secant line of their spherical model) using their parameters (see their Table 1 and 3). We will further discuss the discrepancy between our and their models in Sect. 5. Finally, the best-fit modeled spectra of the four cations along the main cut are shown in red in Fig. 3. In addition, Figure 9 shows our best-fit modeled N₂H⁺ (1–0) spectra overlaid on the data from Tafalla et al. (2004), while Figure 10 shows our best-fit modeled o-H₂D⁺ (1₁₀–1₁₁) spectra overlaid on our full o-H₂D⁺ data. Our reproduced spectra are in good agreement with the observations, suggesting that our model is globally valid for the L 1498 core.

We estimated the chemical timescale of each layer in the onion model of L 1498 with our derived abundance, $n_{\rm H_{2}}$ and $T_{\rm kin}$ profiles. Our method adopted a pseudo-time-dependent gas-phase deuterium chemical model (Pagani et al., 2009b) to simultaneously reproduce the deuteration ratio, $X$(N₂D⁺)/$X$(N₂H⁺), and the abundances of N₂H⁺, N₂D⁺, DCO⁺, and o-H₂D⁺ in each layer. This deuterium chemical model is specialized for starless cores in that each spin-state of H₂ and H${}_{3}^{+}$ isotopologues are included and the complete depletion in heavy species, except for CO and N₂, is assumed. In our model, we did not include the nitrogen chemistry. CO and N₂ would quasi-instantaneously reach chemical equilibrium with our observed species because our chemical network is very small. On one hand, CO and N₂ are the parent molecules reacting with the H${}_{3}^{+}$ isotopologues to form the N₂H⁺ and HCO⁺ isotopologues (i.e., the daughter cations); on the other hand, dissociative recombination of these daughter cations would convert them back to CO and N₂. The abundances of CO and N₂ are left as free parameters in each layer, so we could derive the current CO and N₂ abundance profiles. The other free parameters are the initial OPR of H₂ (OPR_intial(H₂)), the average grain radius ($a_{\rm gr}$), and the cosmic ray ionization rate ($\zeta$).

Figure 6 shows the best-fit chemical model solutions of $X$(N₂D⁺)/$X$(N₂H⁺), $X$(o-H₂D⁺), and $X$(DCO⁺) for each layer at the northeast and southwest sides. The outermost layers were not modeled simply because their large observational error bars cannot well-constrain the chemical solutions. For these chemical solutions, the initial OPR of H₂ is assumed to be its statistical ratio of 3. We adopted the $\zeta$ of $3\times 10^{-17}$ s^-1 throughout the chemical model, which is the best value found by Maret et al. (2013) for their chemical model to reproduce the C¹⁸O (1–0) and H¹³CO⁺ (1–0) integrated emission in L 1498. With the detection of coreshine, one would expect that larger dust grains ($a_{\rm gr}>0.1$ $\mu$m) appear in the core (Pagani et al., 2010; Steinacker et al., 2010; Lefèvre et al., 2014). Although Steinacker et al. (2015) cannot well constrain the maximum grain radius in L 1498 with the coreshine modeling on the IRAC1/2 images, they suggested that the maximum grain radius seems to be greater than $\sim$0.3 $\mu$m. On the other hand, Maret et al. (2013) found that in their chemical model, a power-law grain size distribution grown from the typical grain radius of 0.1 $\mu$m to a maximum grain radius of 0.15 $\mu$m is sufficient for reproducing the C¹⁸O and H¹³CO⁺ line emissions, suggesting that the grain growth does not yet significantly occur everywhere in the entire core. The apparent contradiction of the maximum grain radius reported by Maret et al. (2013) and that reported by Steinacker et al. (2015) arises because the low-$J$ transitions of C¹⁸O and H¹³CO⁺ are optically thick, and as such, they do not trace the core center but rather the low-density outskirts. Thus, these indicate that $a_{\rm gr}\gtrsim 0.3$ $\mu$m at the center of L 1498 while $a_{\rm gr}\approx 0.1$ $\mu$m at the outskirts. The bottom row in Fig. 7 shows the profiles of the grain radius in our chemical model. In order to make the chemical solutions match the observational abundances, we found that the grain radius should be at least 0.4 $\mu$m in the innermost layer, and 0.2 $\mu$m in the second layer. For the rest of the layers, we found matched chemical solutions with the typical grain radius of 0.1 $\mu$m. With the above conditions, we found that the observational abundances of all layers are matched with chemical solutions from 0.07 to 0.31 Ma on the northeast side, and from 0.08 to 0.31 Ma on the southwest side.

The top row in Fig. 7 shows the CO and N₂ abundance profiles from our chemical modeling. Their best-fit values are listed in Table 11. Error bars show the abundance ranges of the possible chemical solutions at each layer but not every combination of $X$(CO) and $X$(N₂) can make the chemical solutions fit the observational cation abundances. We assume a constant ¹²C¹⁶O/¹²C¹⁸O ratio of 560, which is comparable to the local ¹⁶O/¹⁸O ISM ratio of 557$\pm$30 (Wilson & Rood, 1994), to derive a C¹⁸O abundance profile. The C¹⁸O line emission from the central depleted region is too weak to be constrained by radiative transfer calculation. Our chemical model serves as an appropriate approach. At the outer layers, $X$(CO) is less constrained since deuteration is only an upper limit in the chemical models. By contrast, C¹⁸O is not heavily depleted and its $J$=2–1 line emission remains thinner and thus the radiative transfer modeling is applicable. We then used our C¹⁸O (2–1) observation to further determine the outer CO abundance (southwestern layer 4–8, and northeastern layer 2–4). With radiative transfer calculations, we found an outer $X$(C¹⁸O) of $1.9\times 10^{-8}$ can roughly make the C¹⁸O (2–1) spectrum models fit with observations (the last row in Fig. 3), validating our chemical modeling results. We then set the outer $X$(¹²CO) as $1.0\times 10^{-5}$.

We noted that the C¹⁸O (2–1) spectra around the pointing from ($-$10^′′, $-$10^′′) to ($-$40^′′, $-$40^′′) are somewhat lower in intensities compared to our modeled C¹⁸O spectra. It seems to be caused by a more considerable depletion along the line of sight toward these positions. It is possible that the C¹⁸O depletion center shifted outward from the core center that we determined by the emission peak position shared by 850 $\mu$m continuum, N₂H⁺ and o-H₂D⁺ data (also see the C¹⁸O $J$=1–0 and 2–1 integrated intensity maps shown in Fig. 3 from Tafalla et al., 2004). This would hint that the CO isotopologues have a more complex spatial cloud structure along the line of sight.

Similar to the procedure in Lin et al. (2020), we used the NIR/MIR images (Fig. 2) to derive the visual extinction map with the NICER method (Lombardi & Alves, 2001). We used the J, H, and $K_{s}$ images with the $R_{\rm V}$ = 3.1 dust models from Weingartner & Draine (2001) to make an envelope extinction map tracing the L 1498 envelope. On the other hand, we used the H, $K_{s}$, and IRAC2 images with the $R_{\rm V}$ = 5.5B dust models from Weingartner & Draine (2001) to make a core extinction map tracing the L 1498 core. However, the lack of sufficient stars with the color excess $E_{K_{s}-{\rm IRAC2}}$ measurements prevents us from deriving the extinction at the core center. Although many stars appear on the IRAC1/2 images, only a few stars are detected in the $K_{s}$ band in the core region, with which we found that the magnitude differences of $K_{s}$ and IRAC2 bands ($m_{K_{s}}-m_{\rm IRAC2}$) for these stars approach 2.5 mag toward the center. This magnitude difference is already larger than the difference between the completion magnitudes of $K_{s}$ and IRAC2 bands (20 and 18.5 mag, respectively), which explains the disappearance of the $K_{s}$-band counterparts. Therefore, we assume that the $K_{s}$-band magnitude of stars only detected in IRAC2 bands have $m_{K_{s}}\geq m_{\rm IRAC2}+2.5$ mag, and we generated the artificial $K_{s}$-band detection by

1. 1.

   selecting stars that are (a) detected in both IRAC1 and IRAC2 bands to ensure it is not due to contamination in a single band and (b) with the IRAC2 magnitude uncertainty less than 0.2 mag, and

2. 2.

   assigning their $K_{s}$-band magnitude as $m_{\rm IRAC2}+2.5$ mag, i.e., the brightest limit from the above assumption, resulting in a lower limit of $E_{K_{s}-{\rm IRAC2}}$ and $A_{\rm V}$.

After this process, the density of stars with color excess measurements across L 1498 allowed us to convolve the reddening data with a 50^′′-Gaussian beam to produce the core and envelope extinction maps. We combined the core map and envelope map by merging them at 5 mag and the final $A_{\rm V}$ map has been shown in Fig. 2g and the insert of Fig. 8. Stars only detected in the IRAC1/2 bands and with the artificial $K_{s}$-band detection are displayed as black dots on the insert of Fig. 8. Figure 8 also shows the radially averaged $A_{\rm V}$ profiles along the main cut within the 50^′′-wide strip. As a result, the $A_{\rm V}$ values in the central region should be interpreted as the lower limit and therefore our result suggests that $A_{\rm V}$ $\gtrsim 25$ mag at the 50^′′-beam.

With our non-LTE radiative transfer modeling on the N₂H⁺ emission line spectra, we find a central $n_{\rm H_{2}}$ of $1.6_{-0.3}^{+3.0}\times 10^{5}$ cm^-3 and $T_{\rm kin}$ of 7.5${}_{-0.5}^{+0.7}$ K, and a peak $N_{\rm H_{2}}$ of $4.5\times 10^{22}$ cm^-2 toward the L 1498 core. With the submm data, Tafalla et al. (2004) derived $N_{\rm H_{2}}=\mbox{(3--4)}\times 10^{22}$ cm^-2 with $T_{\rm d}=10$ K using the IRAM 30-m MAMBO 1.2 mm map. In contrast, our $N_{\rm H_{2}}$ peak value is slightly larger than that of Tafalla et al.. This discrepancy could be due to our different approach using the radiative transfer modeling of the line emission instead of the continuum emission. The crucial uncertainty in the non-LTE line modeling is the collisional coefficients, which is, however, much smaller than the uncertainty in the dust opacity at the submm wavelength. On the other hand, Tafalla et al. adopted a uniform $T_{\rm kin}$ of 10 K derived from NH₃ (1,1) and (2,2) lines as $T_{\rm d}$ by assuming the gas and dust are coupled. However, the lower critical density of these NH₃ lines ($\sim 2\times 10^{3}$ cm^-3; Pagani et al., 2007) compared to the N₂H⁺ $J$=1–0 line ($1.3\times 10^{5}$ cm^-3 at 10 K) would suggest that N₂H⁺ can trace the inner $T_{\rm kin}$ better than NH₃. Therefore, using 10 K may underestimate the density at the core center.

Another distinction is that our onion model is asymmetric, while Tafalla et al.’s model is spherically symmetric. The top row in Fig. 4 shows the comparison between our number density profiles and those of Tafalla et al.. We can see that our $n_{\rm H_{2}}$ profiles and theirs converge toward the outer southwestern region. While the spherical onion model adopted by Tafalla et al. provides a good approximation for deriving an averaged density profile of the elongated L 1498 core, Their model center ($\Delta\mbox{RA}=-10$^′′, $\Delta\mbox{Dec}=-20$^′′ with respect to ours) was based on the centroid of their 1.2 mm emission map rather than the emission peak. In contrast, our model center is situated at the JCMT SCUBA-2 850 $\mu$m continuum, N₂H⁺ (1–0), and o-H₂D⁺ (1₁₀–1₁₁) emission peaks, including their 1.2 mm continuum peak. In this case, our chosen model center would allow us to capture the peak $N_{\rm H_{2}}$ rather than average it out with the surrounding region.

Meanwhile, our NIR and MIR dust extinction measurements provide another constraint on the density in L 1498. Different from the dust emission measurement, the dust extinction is independent of $T_{\rm d}$, and the NIR/MIR dust opacity is better determined than the submm dust opacity (Pagani et al., 2015; Lefèvre et al., 2016). Figure 8 shows the comparison of the $A_{\rm V}$ profiles derived with the observational data and our onion model along the main cut. The orange/red squares represent the extinction at the 50^′′-beam of the whole L 1498 cloud (i.e., the core and envelope) derived from NIR/MIR images and the central $A_{\rm V}$ value is estimated as the lower limit ($A_{\rm V}$ $\gtrsim 25$ mag; see Sect. 4.4). The grey and blue step curves represent the $A_{\rm V}$ values contributed by the L 1498 core, which are converted from our asymmetric onion-like model by $N_{\rm H_{2}}/A_{\rm V}=9.4\times 10^{20}$ cm^-2 mag^-1 ($R_{\rm V}$ = 3.1; Bohlin et al., 1978), and the blue curve is convolved with the $A_{\rm V}$ beam size of 50^′′ for the comparison with the data. Our onion model represents the core region because it was constructed from the N₂H⁺ data, a molecule confined in the core region. We can see that the $A_{\rm V}$ profile of our onion model matches well with the data in the core region. The $A_{\rm V}$ profile of the data consists of the extinction from the L 1498 core and from an envelope with $A_{\rm V}$ $\approx$ 3 mag. Since our observations of N₂H⁺, N₂D⁺, DCO⁺, and o-H₂D⁺ are associated with the core as these molecules are mostly confined in the core region, we do not need to include the envelope to reproduce their observational spectra shown in Fig. 3. In Appendix C, we show that including an envelope with $A_{\rm V}$ = 3 mag in the radiative transfer model can still reproduce our C¹⁸O (2–1) spectra (the bottom row in Fig. 3), whereas it is necessary for reproducing the C¹⁸O (1–0) spectra obtained by Tafalla et al. (2004) with a lower critical density.

With the interferometer, the ALMA-ACA 1 mm continuum survey, FREJA (Tokuda et al., 2020), found no substructure at the 1 kau scale in the central region of L 1498, suggesting that the central density structure is very flat (characterized by a Plummer-like flattening diameter greater than 5 kau) and an upper limit on $n_{\rm H_{2}}$ of about $3\times 10^{5}$ cm^-3. This non-detection is in agreement with the $n_{\rm H_{2}}$ profile evaluated from our onion model with the central density of $1.6_{-0.3}^{+3.0}\times 10^{5}$ cm^-3 and the central flattened region size of $\sim$11 kau ($R_{0,\rm{SW}}+R_{0,\rm{NE}}$; see the Plummer-like parameters shown on the top row in Fig. 4). Therefore, our solution from the non-LTE N₂H⁺ radiative transfer modeling is compatible with the NIR/MIR dust extinction measurement as well as the 1mm continuum interferometry observation.

We find that N₂H⁺ is significantly depleted toward the L 1498 center with a depletion factor of 5.6${}_{-4.0}^{+2.1}$ with respect to the maximum $X$(N₂H⁺) at the fifth southwestern layer, and another depletion factor of 8.5${}_{-7.6}^{+0.5}$ with respect to the maximum at the second northeastern layer (see Fig. 5). Since the northeastern depletion feature is only resolved by two layers in the onion model, we take the depletion factor measured on the southwest side as the representative measurement in the following discussion (same for the other molecules). Our finding of the N₂H⁺ depletion is opposite to the N₂H⁺ enhancement by a factor of 3 reported by Tafalla et al. (2004) with their radiative transfer modeling. Our central N₂H⁺ abundance is $4.7\pm 1.7\times 10^{-11}$, while their value is $1.7\times 10^{-10}$. This discrepancy could be due to their different physical model and the approximation made in their radiative transfer calculations. Tafalla et al. (2004) derived a central $n_{\rm H_{2}}=9.4\times 10^{4}$ cm^-3, a factor of 1.7 lower than our $1.6_{-0.3}^{+3.0}\times 10^{5}$ cm^-3, and a constant $T_{\rm kin}$(NH₃) = 10 K, higher than our 7.5${}_{-0.5}^{+0.7}$ K (see their $n_{\rm H_{2}}$, $T_{\rm kin}$, and N₂H⁺ abundance profiles on Fig. 4). As previously mentioned, using the temperature derived from the NH₃ lines could be biased by the warmer outer layers and thus one could overestimate the central temperature. In addition, they omitted line overlap from their radiative transfer calculation in their customized version of the MC code, and the hyperfine-structure-resolved collisional coefficients were just not yet available. As a result, the peak intensities of their best-fit N₂H⁺ (1–0) spectra were too strong by $\sim$60% compared to data, but even with our current MC code, adopting their physical model would still result in $\sim$30% stronger intensities at the peaks. In terms of the line-of-sight-integrated measurements, our averaged $X$(N₂H⁺) and N₂H⁺ column density ($N_{\rm N_{2}H^{+}}$) toward the core center are $1.2\times 10^{-10}$ and $5.2\times 10^{12}$ cm^-2, respectively. Our $N_{\rm N_{2}H^{+}}$ value is of the same order of magnitude compared with $1.7\pm 0.7\times 10^{12}$ cm^-2 reported by Crapsi et al. (2005) despite their more simplistic LTE calculations.

Interestingly, both N₂H⁺ depletion and enhancement are seen from chemo-dynamical modeling. Aikawa et al. (2005) performed a self-consistent calculation with a quasistatically contracting Bonner-Ebert sphere, and found that N₂H⁺ is enhanced at a central $n_{\rm H_{2}}$ similar to L 1498 ($1.5\times 10^{5}$ cm^-3) and later becomes depleted at higher densities. Holdship et al. (2017) and Priestley et al. (2018) used an analytical approach for the Bonner-Ebert density evolution, but found that the N₂H⁺ depletion starts earlier in their calculation. Holdship et al. (2017) found an N₂H⁺ depletion factor of $\sim$50 when $n_{\rm H_{2}}$ reaches Tafalla et al. (2004)’s best-fit value for L 1498. Although this discrepancy might relate to their N chemistry details in the models, the depletion case would be preferred by our findings.

In starless cores, the depletion of the N₂H⁺ and HCO⁺ isotopologues is a result of the freeze-out of their parent molecules (N₂ and CO) in the core center. In L 1498, we find depletion factors of 5.6${}_{-4.0}^{+2.1}$ for N₂H⁺ and 17$\pm$11 for DCO⁺. However, in contrast, N₂D⁺ does not exhibit significant depletion toward the core center (see Fig. 5). Here, we do not interpret the inward decrease in $X$(N₂D⁺) on the northeastern side as significant depletion. This may be due to the limited coverage of our spectral observations in the northeastern area (see Fig. 3). The lack of significant depletion in N₂D⁺ is likely because the increase in deuteration of H${}_{3}^{+}$ outweighs the decrease in N₂ due to freeze-out. Regarding the deuteration of N₂H⁺, we have observed a $X$(N₂D⁺)/$X$(N₂H⁺) profile spanning from an upper limit of 0.07 in the outer region to a maximum of 0.27${}_{-0.15}^{+0.12}$ toward the center of L 1498. However, this maximum value is moderate compared to values obtained in the other two starless cores, L 183 (Pagani et al., 2007) and L 1512 (Lin et al., 2020). The authors found a maximum N₂H⁺ deuteration of 0.70$\pm$0.12 in L 183 and 0.34${}_{-0.15}^{+0.24}$ in L 1512. It is worth noting that the formation and destruction of N₂D⁺ and DCO⁺ follow the same chemical scheme in starless cores (Pagani et al., 2011). Given that N₂D⁺ and DCO⁺ are formed under the same conditions of H${}_{3}^{+}$ deuteration, the presence of DCO⁺ depletion suggests that the depletion of CO is greater than that of N₂.

Molecular nitrogen is the parent molecule reacting with the H${}_{3}^{+}$ isotopologues to form N₂H⁺ and N₂D⁺, while CO is the parent molecule also reacting with H${}_{3}^{+}$ isotopologues to form DCO⁺ (Pagani et al., 2011). In starless cores, N₂H⁺, N₂D⁺, DCO⁺ are formed via the above gas-phase paths. Therefore, as we presented in Sect. 4.3, $X$(N₂) and $X$(CO) are free parameters for fitting the observed abundances of their daughter cations. The top row in Fig. 7 shows the abundance profiles of N₂, ¹²CO, and C¹⁸O. Both ¹²CO and C¹⁸O have the same depletion factors of $\sim$20 compared with their maximum abundances of $1.0\times 10^{-5}$ for ¹²CO and $1.9\times 10^{-8}$ for C¹⁸O since we assumed a constant ratio of ¹²C¹⁶O/¹²C¹⁸O as 560 to derive the C¹⁸O profile from ¹²CO profile. By comparing with the standard ¹²CO abundance of $\mbox{(1--2)}\times 10^{-4}$ (Pineda et al., 2010), ¹²CO has a depletion factor of $\sim$200–400. For C¹⁸O, the depletion factor is $\sim$190 with respect to its standard abundance of $1.7\times 10^{-7}$ (Frerking et al., 1982). Conversely, the N₂ abundance does not decrease but even increases toward the core center (Fig. 7). This accounts for the discrepancy between the depletion of DCO⁺ and the enhancement of N₂D⁺ toward the core center.

In our chemical modeling, we find that the dust grain radius increases from 0.1 $\mu$m in the outer region to 0.4 $\mu$m in the core center. This grain growth slows down the depletion of CO and N₂ due to freeze-out onto dust grains at the core center by a factor of $\sim$3, as the depletion timescale is proportional to the dust grain radius (Maret et al., 2013). However, the freeze-out rates of CO and N₂ are suggested to be similar because of their identical mass and similar binding energies (Bisschop et al., 2006). Also, the cosmic ray ionization is too weak in the dense core to compensate for the depletion of CO and N₂. The discrepancy between the N₂ and CO abundance profiles may suggest that the N₂ production from the atomic N is still active such that the N₂ production can compensate for the depletion of N₂, similar to what Pagani et al. (2012) found in the core center of L 183. In contrast, the CO production has reached a steady-state in the diffuse cloud region (Oppenheimer & Dalgarno, 1975), making the freeze-out effect the dominant mechanism regulating the CO abundance in the dense core region.

With the envelope tracers of CO and HCO⁺, Maret et al. (2013) performed pseudo-time-dependent chemical modeling to study their depletion in L 1498. They constrained their model by fitting the integrated intensity profiles and the central spectra of C¹⁸O (1–0) and H¹³CO⁺ (1–0). To reproduce their observations, they found a maximum grain radius of 0.15 $\mu$m in the core center. Our maximum grain radius of 0.4 $\mu$m is greater than their value, likely because the four molecular cations (N₂H⁺, N₂D⁺, o-H₂D⁺, and DCO⁺) we used are not heavily depleted and their emissions are optically thinner than the above envelope tracers. Therefore, our chemical modeling is more sensitive to the grain radius in the central region. Despite this difference, when they included such grain growth in their chemical modeling, they derived a central CO abundance of $\sim 6\times 10^{-7}$, which is consistent with our value of 5.0${}_{-1.0}^{+26.6}\times 10^{-7}$. In addition, their chemical network includes nitrogen chemistry. They found that a certain fraction of the nitrogen should be locked in NH₃ ices to reduce the gas-phase abundances of N₂, and atomic N. As a result, they obtained a N₂H⁺ abundance of $4\times 10^{-10}$, which is closer to $1.7\times 10^{-10}$ derived by Tafalla et al. (2004) (also see Sect. 5.2 for the comparison of our N₂H⁺ abundance profile with Tafalla et al.). However, our central $X$(N₂H⁺) is even lower at $4.7\pm 1.7\times 10^{-11}$, and our central $X$(N₂) at 4.0${}_{-2.0}^{+2.3}\times 10^{-6}$ is larger than Maret et al. (2013)’s $2\times 10^{-7}$. Therefore, more detailed modeling is necessary to understand the nitrogen chemistry in L 1498.

Figure 6 shows that between 0.07 to 0.31 Ma, the chemical solution of each layer appears at the observed abundances within the uncertainties. One important characteristic of deuterium chemistry in starless cores is that the deuteration fractionation can serve as an effective chemical clock. This is because it primarily depends on time and is influenced by a few initial parameters, including the initial ortho-to-para ratio of H₂ (OPR_intial(H₂)), the cosmic ray ionization rate ($\zeta$), and the average grain radius ($a_{\rm gr}$). It is less affected by temperature and dynamical parameters such as the turbulent Mach number and the mass-to-magnetic flux ratio (Pagani et al., 2013; Kong et al., 2016; Körtgen et al., 2017). Therefore, given reasonable values/ranges of OPR_intial(H₂), $\zeta$, and $a_{\rm gr}$, the timescale of deuteration fractionation (i.e., the chemical timescale) provides a means to measure the age of a contracting core.

With pseudo-time-dependent chemical analysis, we have determined the chemical timescales of the deuteration in each layer in Sect. 4.3. Since the physical structure is kept constant at the present state in the pseudo-time-dependent models during the evolution of the chemical composition, our chemical timescales do not include the physical evolution of the core. Therefore, the determined chemical timescales at the fixed densities serve as a lower limit to the the core age. Particularly, we take the chemical timescale of a less evolved outer layer, where $X$(N₂D⁺)/$X$(N₂H⁺) ratio and $X$(o-H₂D⁺) are measured, to represent the lower limit of the core age. It follows that the lower limit on the core age of L 1498 is $\sim$0.16 Ma determined by the fourth southwestern layer, which is the last layer with the N₂D⁺/N₂H⁺ ratio measurement, and this lower limit is conservatively chosen at the beginning of the time range of this layer. We would like to note that we use the $X$(N₂D⁺)/$X$(N₂H⁺) ratio and $X$(o-H₂D⁺) to determine the chemical timescale of a layer, while we use $X$(DCO⁺) to derive the CO abundance. Since the fourth southwestern layer has no o-H₂D⁺ detection but only an upper limit on its abundance, the interpretation for this time limit is the time duration for reaching $X$(N₂D⁺)/$X$(N₂H⁺) under the assumption of $X$(o-H₂D⁺) at the maximum of its possible range, that is the fastest time of that layer.

Although deuteration fractionation depends on OPR_intial(H₂), $\zeta$, and $a_{\rm gr}$, we assumed OPR_intial(H₂) to be the statistical ratio of 3. Additionally, we adopted $\zeta=3\times 10^{-17}$ s^-1 from Maret et al. (2013) and the typical $a_{\rm gr}$ of 0.1$\mu$m for the the fourth southwestern layer (see Sect. 4.3). Therefore, L 1498 is presumably older than 0.16 Ma. If we assume another OPR_intial(H₂) of 0.5 as measured in the diffuse medium (Crabtree et al., 2011), we find another, but similar, age lower limit of 0.13 Ma. In the fast collapse model, the contraction timescale of the core can be characterized by the free-fall time, $t_{\rm ff}$.

where $G$ is the gravitational constant and $\rho$ represents the mass density, which can be expressed as $\rho=\mu_{\rm H_{2}}m_{\rm H}n_{\rm H_{2}}$, where $m_{\rm H}$ is the mass of atomic hydrogen and $\mu_{H_{2}}=2.8$ denotes the mean molecular weight per H₂. Considering that dense starless cores are typically contracted from lower density gas at $n_{\rm H_{2}}=10^{4}$ cm^-3, $t_{\rm ff}=0.31$ Ma, which is compatible with the lower limit of the core age of L 1498 at 0.16 Ma. On the other hand, the core contraction timescale in the slow collapse model can be associated with the ambipolar diffusion timescale in the presence of the magnetic field support. The ambipolar diffusion timescale can be longer than a free-fall time by an order of magnitude (Shu et al., 1987; Tassis & Mouschovias, 2004; Mouschovias et al., 2006; McKee & Ostriker, 2007). Although a lower limit of the core age at 0.16 Ma does not conclusively rule out the possibility of a much longer actual core age, such a long contraction timescale does not seem to be favored by our findings. Particularly, significant inward motions have been detected toward the L 1498 envelope via CS lines (Lee et al., 2001; Lee & Myers, 2011) and HCN lines (Magalhães et al., 2018). Therefore, our result suggests that the contraction of L 1498 follows the fast collapse model, indicating that L 1498 likely formed rapidly.

Maret et al. (2013) and Jiménez-Serra et al. (2021) also derived the chemical timescales of L 1498 but with different molecules. Maret et al. (2013) derived a chemical timescale of L 1498 with a range of 0.3–0.5 Ma based on the CO freeze-out modeling with the C¹⁸O (1–0) line, while Jiménez-Serra et al. (2021) derived a chemical timescale as 0.09 Ma from their COMs observation toward the methanol emission peak at the outer layers of L 1498. Because both studies employed the pseudo-time-dependent chemical analysis, their chemical timescales are also lower limits on the age of L 1498. It is noteworthy that although the timescale derived by Maret et al. is longer than our deuteration timescale, using their value as the lower limit of the core age would still favor the fast collapse model. In fact, $X$(CO) at the core center derived from our deuteration chemical modeling aligns with their derived CO value (see Sect 5.3). However, assessing CO depletion toward highly embedded and more evolved cores can be challenging due to the contamination along the sightline and the severe CO depletion. In such cases, the deuteration fraction is easier to measure and can provide valuable insights into the chemical evolutionary history of the core contraction. Therefore, it is demonstrated that different chemical clocks provide a comprehensive point of view on the core ages. On the other hand, in the chemical modeling of COMs, the surrounding environment could be more important than chemical evolutionary effects (Lattanzi et al., 2020). Then COMs may not act as a chemical clock but as an indicator of the environmental conditions.

In this work, we discuss the low-mass core formation in the context of contrasting fast and slow collapse models of relatively isolated entities in the parent clouds. However, we note that the interactions between starless cores and their surrounding clumps/filaments may lead to more dynamical scenarios, such as fragmentation into smaller entities and mass accumulation/accretion from parent filamentary structures, suggesting a hierarchical nature (e.g., Vázquez-Semadeni et al., 2019). This complexity potentially aligns with scenarios favoring fast collapse. Recently, Bovino et al. (2021) performed a 3D turbulent MHD simulation coupled with deuterium chemistry to study a collapsing filament, fragmenting into dense cores. They proposed the o-H₂D⁺ to p-D₂H⁺ ratio as an alternative chemical clock to the N₂D⁺/N₂H⁺ ratio for tracing OPR(H₂). Their APEX o-H₂D⁺ (1₁₀–1₁₁) and p-D₂H⁺ (1₁₀–1₀₁) observations of six starless cores in the Ophiuchus region demonstrate that fast collapse models align well with the observed o-H₂D⁺/p-D₂H⁺ ratios with core ages at $\sim$0.05–0.2 Ma, indicating that these Ophiuchus cores can form more rapidly.