{'fraction_non_alphanumeric': 0.04301530963606505, 'fraction_numerical': 0.04089632886581554, 'mean_word_length': 3.951215247420878, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Physical and chemical modeling of the starless core L 1512 ?

Context. The deuterium fractionation in starless cores gives us a clue to estimate their lifetime scales, thus allowing us to distinguish between dynamical theories of core formation. Cores also seem to be subject to a differential N2 and CO depletion, which was not expected from the models. Aims. We aim to create a survey of ten cores to estimate their lifetime scales and depletion profiles in detail. After describing L 183, located in Serpens, we present the second cloud of the series, L 1512, from the star-forming region Auriga. Methods. To constrain the lifetime scale, we performed chemical modeling of the deuteration profiles across L 1512 based on dust extinction measurements from near-infrared observations and nonlocal thermal equilibrium radiative transfer with multiple line observations of N2H, N2D, DCO, C18O, and 13CO, plus H2D (110–111). Results. We find a peak density of 1.1×105 cm−3 and a central temperature of 7.5±1 K, which are higher and lower, respectively, compared with previous dust emission studies. The depletion factors of N2H and N2D are 27+17 −13 and 4 +2 −1 in L 1512, which are intermediate between the two other more advanced and denser starless core cases, L 183 and L 1544. These factors also indicate a similar freeze-out of N2 in L 1512, compared to the two others despite a peak density one to two orders of magnitude lower. Retrieving CO and N2 abundance profiles with the chemical model, we find that CO has a depletion factor of ∼430–870 and the N2 profile is similar to that of CO unlike that toward L 183. Therefore, L 1512 has probably been living long enough so that N2 chemistry has reached steady state. Conclusions. N2H modeling is necessary to assess the precise physical conditions in the center of cold starless cores, rather than dust emission. L 1512 is presumably older than 1.4 Myr. Therefore, the dominating core formation mechanism should be ambipolar diffusion for this source.


Introduction
Starless cores are the sites of future star and planet formation. Although we know gravitational collapse plays a main role during star formation, the details of the collapse process, particularly the timescale, are not yet well understood. There are two major dynamical theories with different dominating mechanisms. One theory is the ambipolar diffusion-controlled scenario, advocating that molecular clouds evolve quasi-statically to form starless cores after 1-10 Myr (Shu et al. 1987;Tassis & Mouschovias 2004;Mouschovias et al. 2006;McKee & Ostriker 2007). The other theory is the supersonic turbulenceinduced scenario, in which cores can form faster through gravoturbulent fragmentation with a free-fall timescale less than 1 Myr (Mac Low et al. 1998;Padoan & Nordlund 1999;Hartmann et al. 2001;Federrath & Klessen 2012). Measuring the age of the Based in part on observations obtained with WIRCam, a joint project of CFHT, Taiwan, Korea, Canada, France, at the Canada-France-Hawaii Telescope (CFHT), which is operated by the National Research Council (NRC) of Canada, the Institute National des Sciences de l'Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. The observations at the Canada-France-Hawaii Telescope were performed with care and respect from the summit of Maunakea which is a significant cultural and historic site. cores would therefore allow us to select the correct scenario but the task is difficult. Time-dependent chemical analysis has long been thought to be a solution to this problem but unknown initial conditions too often hamper a precise determination of the elapsed time from the molecular cloud assembly to the formation of the core (see, e.g., Nilsson et al. 2000 for a discussion on the sensitivity of a chemical model to initial conditions). To circumvent the initial condition problem, Pagani et al. (2009bPagani et al. ( , 2013 proposed an approach based on the N 2 D + /N 2 H + deuteration profile across the starless core. These works showed that the profile itself is mostly dependent on the time it took to form the starless core, showing a clear distinction between fast collapsing cores with steep profiles and slow collapsing cores with much higher deuteration level and much flatter profiles. This is possible thanks to the slow decrease of the ortho/para-H 2 (thought to be originally produced on grains in the 3:1 ratio), which is comparable to the lifetime of starless cores and because a high abundance of ortho-H 2 prevents the deuteration to proceed (Pagani et al. 1992a(Pagani et al. , 2009b(Pagani et al. , 2011. Varying the initial conditions did not change the clear-cut separation between the two types of profiles. L 183 was shown to clearly belong to the fast collapse timescale (in a free-fall time; Pagani et al. 2013). More recently, other dynamical studies also based on deuteration (but Article number, page 1 of 22 arXiv:2002.01346v2 [astro-ph.GA] 18 Feb 2020 not on N 2 H + deuteration radial profile evolution) have been published (Kong et al. 2015;Lackington et al. 2016;Körtgen et al. 2017;Körtgen et al. 2018;Giannetti et al. 2019); but these works have contradictory conclusions concerning the typical lifetime of cores.
The depletion effect on heavy species (e.g., CO, CS, SO) makes these species ineffective tracers of n H 2 and T kin in starless cores. Undepleted high-density tracers with low excitation energies are needed to reveal the inner structure of starless cores. NH 3 and N 2 H + line emissions are commonly found to be spatially anti-correlated with CO line emission among low-and high-mass starless cores, which suggests that light nitrogenbearing species are better tracers of the inner structure of starless cores (e.g., Tafalla et al. 2002;Fontani et al. 2006;Pagani et al. 2007). The N 2 H + (J=1-0) rotational line provides a better constraint on both n H 2 and T kin in dense cores than the NH 3 (J, K)=(1,1) inversion line, since the former has a higher critical density (5 × 10 5 cm −3 ) compared with the low critical density (2 × 10 3 cm −3 ) of the latter (Pagani et al. 2007). However, when density is higher than ∼10 6 cm −3 , N 2 H + and its isotopologs start to decrease in abundance, which suggests that their parent molecule, N 2 , has frozen out onto dust (e.g., Bergin et al. 2002;Pagani et al. 2005Pagani et al. , 2007Pagani et al. , 2012. In this case, since the CO-(and N 2 -) depleted environment enhances the deuterium fractionation, H + , H + 3 , and their isotopologs would be the most abundant ions. Then H 2 D + becomes the best tracer for the innermost region in starless cores owing to the low excitation energy (∆E/k = 18 K) of the ortho-H 2 D + ground transition, J K a K c =1 10 -1 11 (Khersonkii et al. 1987;Pagani et al. 1992b;Bergin et al. 2002;van der Tak et al. 2005;Caselli et al. 2008;Pagani et al. 2009a).
Many studies have estimated the depletion level in starless cores from column density maps, which imply a uniform density along the line of sight (e.g., Caselli et al. 1999;Bacmann et al. 2002;Bergin et al. 2002;Kong et al. 2016). However, the depletion level at the center of cores could be underestimated because of line-of-sight integration of their undepleted outer regions. Pagani et al. (2007) performed spherically symmetric "onion shell" modeling of N 2 H + data toward L 183, and showed a volume depletion of N 2 H + by a factor of 6 +13 −3 , which is larger than the depletion factor of <2 toward the same source estimated from the column density in their earlier work . Similarly, Pagani et al. (2012) show that CO depletion reaches a factor of 2000 instead of the lower limit of 43 they found in an earlier work .
With the tools we developed (Pagani et al. 2007(Pagani et al. , 2009b(Pagani et al. , 2012(Pagani et al. , 2013 we can now expand our work, solely based on L 183 until now, to other clouds. We have selected a list of ten clouds spread among several star-forming regions (Auriga, Taurus, Serpens, and Rho Ophiuchus) for which we would like to study the dust content and dust growth (via the coreshine effect, Pagani et al. 2010b;Steinacker et al. 2010), the CO and N 2 depletion profiles and the N 2 D + /N 2 H + deuteration profile. Dust growth, depletion profile, and deuteration profile are three relatively independent ways to measure the age of starless cores and therefore reinforce our confidence in the age determination of each of these cores. This allows us to address the discrepancies revealed by the various works on dynamical deuteration mentioned above. Each core is studied in detail. From the comparison of all the cores, we gain a better understanding of starless core formation and beyond, of low-mass star formation. This first paper focuses on the gas phase of L 1512.
L 1512 is a relatively isolated nearby starless core that is located at a distance of 140 pc near the edge of the Taurus (Myers et al. 1983;Launhardt et al. 2013). At 1 pc scale, L 1512 is surrounded by CO gas elongated in the north-south direction with a sharp edge in the southern region (Falgarone et al. 1998(Falgarone et al. , 2001. The L 1512 envelope is found to have no significant infall but probably oscillatory (or expansion) motion (Lee et al. 2001;Sohn et al. 2007;Lee & Myers 2011). Kim et al. (2016) further propose that the southern sharp edge is caused by external radiation from an A-type main-sequence star located at 97 pc away from us; this means a distance of ∼40 pc between the two if L 1512 is at the Taurus-Auriga distance, which is a separation too large for a significant impact. The L 1512 core has a globular shape with a size of ∼0.08 pc, and is quiescent with line widths of ∼0.2 km s −1 (Cox et al. 1989;Fuller & Myers 1993;Caselli et al. 1995Caselli et al. , 2002. The comparison between the HC 3 N and NH 3 line widths shows a small nonthermal velocity dispersion of 0.038 km s −1 and an averaged gas kinetic temperature of 11.6 K (Fuller & Myers 1993). N 2 H + observations reveal a nondepleted morphology at 0.1 pc scale, but is found to be depleted in the innermost region of L 1512 with an offset between the N 2 H + and dust emission peaks (Lee et al. 2003).
In this paper, we describe our observations in Sect. 2 and present them in Sect. 3. In Sect. 4, we perform the 1D spherically symmetric nonlocal thermal equilibrium (non-LTE) radiative transfer model and time-dependent chemical model. In Sect. 5, we discuss the lifetime and the possible formation mechanism of L 1512. We summarize our results in Sect. 6.

Observations and data reduction
Observational parameters are summarized in Table 1

IRAM 30 m observations
We observed L 1512 using the IRAM 30 m telescope in December 2013, May and October 2014, and September 2017. Five molecular lines were observed in frequency-switching mode using the dual polarization Eight MIxer Receiver (EMIR), the Versatile Spectral Autocorrelator (VESPA), and the Fourier Transform Spectrometer (FTS); these lines are N 2 H + (1-0), N 2 H + (3-2), N 2 D + (2-1), DCO + (2-1), and DCO + (3-2). The frequency resolution was adapted to keep the velocity sampling close to 50 m s −1 . Both polarizations are systematically averaged. The spectra were observed on a 12 grid in an east-west strip 12 south of the core center (RA, Dec) J2000 = (5 h 04 m 07 s .5, +32 • 43 25 . 0) and along a north-south strip across the core center. Pointing was regularly checked every 1.5 hours and pointing error was monitored to be less than 3 . Data were subsequently folded and baseline subtracted with CLASS 1 . The original cut passes 12 south of the true center because when the project started the CFHT and SCUBA-II maps were not yet available and the core center position was not accurately known. Complementary observations in September 2017 (Director Discretionary Time) were performed. Two short east-west strips across the core center and 12 above the center, and one short strip along the north-south strip (mentioned above) for N 2 H + (1-0) and DCO + (2-1) allowed us to check that the modeling of the core is still consistent with the included offset. Figure 1 shows the pointings of spectral line observations, and the observational parameters are summarized in Table 1. We complemented our observations with on-the-fly (OTF) maps of N 2 H + (1-0) and C 18 O (2-1) obtained by Lippok et al. (2013), and multipointing mosaic maps of C 18 O (1-0) and 13 CO (1-0) obtained by Falgarone et al. (1998). Lippok et al. (2013) and Falgarone et al. (1998) give observational details. The above data are regridded to supplement two orthogonal strips (Figs. 3 and 4) and measure the spectral center velocity (Fig. 5d) using a Gaussian fit routine in CLASS. We also produce the moment maps of N 2 H + (1-0) and C 18 O (1-0) (Figs. 2h,i and Figs. 5a,b).
1 http://www.iram.fr/IRAMFR/GILDAS Among these, the N 2 H + moment maps include our observations and the OTF data.

JCMT observations
The N 2 H + (4-3) and H 2 D + (1 10 -1 11 ) observations were carried out in December 2015 using the JCMT 15 m telescope and the 16 pixels HARP receiver (two pixels of which are non-functioning) in frequency-switching mode. Thanks to the proximity of the two lines (250 MHz), they are routinely observed together with two sub-bands from the Auto Correlation Spectral Imaging System (ACSIS) spectrometer. Data are converted to the CLASS format to be reduced with CLASS (folding and baseline subtraction). Observational parameters are summarized in Table 1, and the pointings of spectral line observations are shown in Fig. 1. SCUBA-II observations were retrieved from the JCMT archive 2 . They are part of the M13BC01 project (PI, James di Francesco). As far as we know, these data have not been published previously.

GBT observations
We observed L 1512 using the GBT 100 m off-axis telescope in November 2014. We used the MM1 (67-74 GHz) and MM2 (73-80 GHz) W-band dual polarization sub-band receivers in in-band frequency-switching mode with the Versatile GBT Astronomical Spectrometer (VEGAS) backend. Pointing was checked every hour and each session started with a surface shape verification by using Out Of Focus Holography (OOF) on a strong source. Data were preprocessed (calibrated, folded, and baseline subtracted) in the GBTIDL data reduction program and converted to CLASS format for subsequent reduction. Thanks to the very smooth baseline of the receiver (due to the use of large band amplifiers instead of mixers), the preprocessing of the data allows us to treat separately the ON and OFF observations (i.e., total power mode) and gain a factor of √ 2 on the noise (the method will be presented in a subsequent paper).

CFHT observations
The CFHT Wide InfraRed CAMera (WIRCAM) was used with the wide filters J, H, and K s to observe the source on the night of 26 December 2013. To preserve the extended emission (scattered light known as cloudshine ;Foster & Goodman 2006), we used a SKY-TARGET-SKY nodding mode to subtract the atmospheric infrared emission. Seeing conditions were typically 0 . 8. Data reduction was performed at the (now closed) TERAPIX center using a specific reduction method to recover the extended emission. The integration times were typically of 0.5 to 1 hour on-source per filter, to reach a completude magnitude of 21.5 (J band) to 20 (K s band). The extended emission will be analyzed in a subsequent paper.

Spitzer observations
Spitzer observations are collected from the Spitzer Heritage Archive (SHA) 3 . L 1512 was observed in two programs with Spitzer InfraRed Array Camera (IRAC), Program id. 94 (PI, Charles Lawrence) and Program id. 90109 (PI, Roberta Paladini); the second program occurred during the warm period of the mission (only 3.6 and 4.5 µm channels still working). The P94 data were already discussed in Stutz et al. (2009), while the warm mission data are analyzed in Steinacker et al. (2015). These two papers provide the presentation of the observations. During the warm mission, deep observations of the cloud were performed and a completed magnitude of 18 in band I2 (4.5 µm) was reached. Figure 2 shows the continuum maps of L1512 at near-infrared (NIR), mid-infrared (MIR), and submillimeter (submm) wavelengths. The NIR maps of L1512 at J, H, and K s bands (Figs. 2ac) reveal extended emission that is distributed from an annular shape at J band and progressively merges to a concentrated shape at K s band. These structures are caused by the scattered light of ambient interstellar radiation field by dust grains (a gr ∼ 0.1 µm) at the periphery of L1512. Similar structures in dark clouds in the Perseus molecular complex have been found and studied by Foster & Goodman (2006) and Juvela et al. (2006). Foster & Goodman (2006) introduced this scattering effect as cloudshine, and proposed to use it as a general tracer of the column density of dense clouds. Since the wavelength of the maximum scattering cross section increases with the size of dust grains, the different spatial spread of the cloudshine of L1512 at the three bands also indicates that the dust grains have the smallest size at the outskirts, and become larger toward the inner region as discussed in Steinacker et al. (2010) and Lefèvre et al. (2014Lefèvre et al. ( , 2016.

Continuum maps
On the other hand, the MIR maps of L1512 from IRAC1 and IRAC2 bands (Figs. 2d,e) show a more compact emission feature, while the IRAC4 map (Fig. 2f) shows absorption instead. The absorption feature is mainly caused by the bright PAH emission from the diffuse medium, the extinction of which becomes superior to the inner scattered light contribution (Lefèvre et al. 2016). The similar spatial distributions between the compact 3 https://sha.ipac.caltech.edu/applications/Spitzer/SHA/ emission and absorption indicate that the emission also comes from the central core. Pagani et al. (2010b) and Steinacker et al. (2010) firstly recognized this phenomenon in their sample of nearby starless and Class 0/I protostellar cores. By analogy to cloudshine, they named this effect coreshine. With the benefit of the smaller beam sizes of IRAC maps compared with those of JCMT SCUBA maps, we determined the center of L1512 core at (RA, Dec) J2000 = (5 h 04 m 07 s .5, +32 • 43 25 . 0) using IRAC4 band (8 µm) with a beam size of 2 , and denote it as a cross shown in each panel of Fig. 2. Since L1512 is a quiescent starless core, each hyperfine component of N 2 H + (1-0) and N 2 D + (1-0) is well separated with line widths of ∼0.2 km s −1 . Along the main horizontal cuts, the coverage of the N 2 H + (1-0) spectra includes the whole starless core with the intensity peak toward the central region. Likewise, the spectra of N 2 H + (3-2), N 2 D + , and DCO + also peak toward the central region, which suggests that these four cations trace the dense core region. In addition, our highest observed transition of N 2 H + , the J=4-3 line, shows no signal according to the 3σ criterion with a rms noise level of ∼30 mK because of the low kinetic temperature in the core region and low density compared to the requested critical density for this transition (two orders of magnitude lower). Among the four cations, the ortho-H 2 D + (1 10 -1 11 ) line is detected in a smaller spatial region (∆RA = −17.2 to 42.8 ), which indicates that ortho-H 2 D + traces the innermost region of the L 1512 core. On the other hand, the intensity of the C 18 O (1-0) and (2-1) spectra are relatively constant across the inner core, which is a clear indication that C 18 O is depleted across the whole core.

Molecular emission lines
Along the main vertical cuts, Fig. 4 shows that N 2 H + , N 2 D + , and DCO + are also depleted in the innermost region. Their spectral peaks drop between the region from ∆Dec = −24 to +12 for which the intensity decrements range from 0.04 to 0.15 K with ∼2-5σ significant levels, while they should increase to follow the increase of column density toward the center if there were no depletion. Similarly, the ortho-H 2 D + spectra show a roughly constant intensity across the peak, which is a sure sign of depletion since the line is optically thin and should be roughly proportional to the H 2 column density. The depletion effects are also found in the integrated intensity maps of N 2 H + (1-0) and C 18 O (1-0) shown in Figs. 2h and i. The C 18 O emission clearly shows depletion toward the center of the cloud because the peak emission is situated on the outskirts, while an offset between N 2 H + emission peak and dust peak suggests that N 2 H + also depletes in the innermost region.

Velocity structure
Figures 5a and b show the mean velocity maps of L 1512, which reveals a velocity gradient in the north-south direction. The constant gradient inside the core itself is compatible with a rigid-like rotation feature or could be due to an oscillation of the whole filament. Figure 5d shows the central velocity of each spectral line along the main vertical cuts fitted with the CLASS Gaussian fitting task. The hyperfine structures of the N 2 H + , N 2 D + , and DCO + rotational lines are taken into account for the fitting. We see that the N 2 H + (1-0) emission appears to have a rather uniform velocity gradient of 2.26 ± 0.04 km s −1 pc −1 between ∆Dec = −48 and +36 , which is compatible with a rigid body rotation. The N 2 H + (3-2), N 2 D + (2-1), DCO + (2-1), (and possibly ortho-H 2 D + (1 10 -1 11 ) although its velocity measurement is hampered by a low signal-to-noise ratio) emissions also have similar dynamical behaviors. The C 18 O (1-0) and 13 CO (1-0) data from Falgarone et al. (1998) partially follow the ionic species velocity drift but with an offset of ∼50 m s −1 toward the center. Because the 13 CO emission traces a more extended region than C 18 O, the similarity between the C 18 O and 13 CO emissions might be due to C 18 O depletion in the core region and its emission being dominated by the L 1512 envelope. The CO isotopolog redshift by ∼50 m s −1 toward the center could be explained by a contraction of the envelope, if the line is optically thick enough to mask the corresponding blueshifted component that we would expect on the back side of the cloud; however, this seems difficult to justify for the C 18 O (1-0) transition. The other possible explanation is that the envelope is experiencing an oscillation in which the core is coupled to it with a delay and is presently toppling over. This kind of oscillation has been proposed to explain the Orion and the G350.54+0.69 filaments (Stutz et al. 2018;Liu et al. 2019) and has also been invoked for L 1512 itself (Lee et al. 2001;Sohn et al. 2007;Lee & Myers 2011).

Analysis
Our first step is to analyze the dust absorption maps from 1.2 µm to 4.5 µm to estimate the total column density of gas+dust all over the cloud. We converted this map into a sphere on firstorder approximation to analyze the line emissions. We adopted a correspondingly 1D spherically symmetric non-LTE radiative transfer code (Bernes 1979;Pagani et al. 2007) to reproduce our observed emission line spectra. Since the shape of L 1512 is globular, we can assume that its distributions of volumetric number density (n H 2 ), gas kinetic temperature (T kin ), volumetric relative abundances with respect to H 2 of the observed species (X species = n species /n H 2 ) and turbulent velocity (V turb ) can be approximated by an onion-shell structure composed of multiple concentric homogeneous layers. We note however that the cloud is slightly dissymmetric and we used two slightly different models to describe the east and west sides of the cloud. In addition, rotational velocity field (V rot ) and radial velocity field (V rad ) can also be applied to each layer. We chose the widths of the layers to be 12 (=1,680 AU), which is the smaller spacing of the two pointing-grids (Fig. 1). The advantage is that we can determine the physical parameters and chemical abundances at each layer sequentially from the outermost to the innermost layer by sampling sightlines of progressively decreasing radius along a cut across L 1512.
Our procedure to determine n H 2 , T kin , X species , and V turb of each layer is as follows. First, the n H 2 profile was determined independently from an extinction map. Second, we used N 2 H + data to determine their T kin , X(N 2 H + ), and V turb profiles. Third, we assumed that all the other species have the same T kin and V turb profiles as those of N 2 H + to determine the abundance profiles of N 2 D + , DCO + , and ortho-H 2 D + . Given the abundance profiles of the various ions, we adopted a pseudo time-dependent chemi-Article number, page 7 of 22  cal model (Pagani et al. 2009b) to estimate the lifetime scale of each layer with our derived n H 2 and T kin , and we also derived the abundance profiles of their parent species, CO and N 2 .

Visual extinction
We adopted the Near-Infrared Color Excess Revisited technique (NICER; Lombardi & Alves 2001) to derive the visual extinction, A V , with color excesses of stars from NIR and MIR bands.
Since the NIR data preferentially trace the diffuse region, we used the R V = 3.1 model from Weingartner & Draine (2001) to derive the A V map of the envelope (Fig. 6a) using J, H, and K s bands. On the other hand, the MIR data preferentially trace the dense region. We alternatively used the R V = 5.5 case B model from Weingartner & Draine (2001) to derive A V map of the core region (Fig. 6b)  gion better (by including spherical grains up to 10 µm in size), although not perfectly (Cambrésy et al. 2011;Ascenso et al. 2013). Thanks to the deep sensitivity of the CFHT and Spitzer observations, the density of stars across the cloud core is high enough to convolve the reddening data with a 2D Gaussian profile of half-power beam width (HPBW) as small as 30 to construct the maps. The dense core map traces the central region better than the diffuse envelope map but could not properly trace the outer low-A V region, and vice versa. In order to combine these two maps, we calculated the azimuthally-averaged A V profiles of the envelope and core maps shown in Fig. 6d from the 1 horizontal strip across the center (Figs. 6a,b) to avoid the cometary tail. We could see that the two profiles merge at ∼2.5 mag and become similar in the outer region. Therefore, an A V threshold of 2.5 mag is used to combine the central region from the core map and the outer region from the envelope map, and the result is shown in Fig. 6c.

Density profile
The volume density (n H 2 ) profile of L 1512 is modeled with the spherically symmetric Plummer-like profile, where n 0 is the central density, R 0 is the characteristic radius, η is the power-law index of n H 2 profile as r R 0 , and R edge is set as 320 (=44,800 AU) to cover the whole core. Since the contribution from the diffuse envelope is small compared with the dense core, we simply adopted the R V = 5.5B model for the conversion of the extinction, A V , into column density, N H 2 . Therefore we obtain a relation between A V and N H 2 (Bohlin et al. 1978), In order to derive n 0 , R 0 and η in the volume density model of L 1512, we produced an A V model map by (1) integrating a Plummer-like n H 2 model along the line of sight and convolving it with the same beam size as the combined A V map (Fig. 6c) to obtain an N H 2 model map, and then (2) derive an A V model map using Eq. 2. Therefore, we could perform a fitting on the azimuthally-averaged profiles of the A V model map and the combined A V map with the Levenberg-Marquardt algorithm included in the Python scipy.optimize.curve_fit function. Figure 7 shows the fitting results for the east and west sides of L 1512 using the A V computed in the strip shown in Fig. 6c, which spatially covers the main horizontal cuts. In order to ensure that the volume density models of both sides agree with each other at the center, n 0 was previously fitted by an averaged profile from both sides. Then n 0 was fixed in the individual fitting for each side. These determined parameters and volume density profiles are shown in the first row in Fig. 8. We see that the horizontal n H 2 profile is skewed in such a way that it is steeper on the west side than the other side. The volume densities drop to ∼10 4 cm −3 at both edges of our spectral observation coverage.

Radiative transfer applied to the onion-shell model
We assumed an onion-shell model comprised of nine and six 12 layers in thickness on the east and west sides, respectively.
Then we reproduced the spectra of N 2 H + (1-0, 3-2, and 4-3) along the main horizontal cuts with the 1D spherically symmetric non-LTE radiative transfer code (Pagani et al. 2007) by varying T kin , X(N 2 H + ), V turb , V rot , and V rad at each layer. Here we used the para-H 2 -N 2 H + collisional coefficients instead of the original He-N 2 H + hyperfine collisional coefficients . We also used specific para-H 2 -N 2 D + collisional rate coefficients for the deuterated isotopolog (see Appendix B). One of the important features of this code is that the line overlap between hyperfine transitions at close frequencies is taken into account to calculate the correct spectral line shape and excitation. Among the input parameters, V turb , V rot , and V rad are determined from observations. We set both V rot and V rad to zero, since the L 1512 core has no obvious infall/expansion motion and only a vertical rotational component, while our model follows the RA (horizontal) structure. For V turb , we found that a uniform turbulent velocity of 80 m s −1 could reproduce the spectral line widths along the cuts. Afterward, we performed an iterative spectral fitting procedure to obtain T kin and X(N 2 H + ) profiles, taking care to eliminate nonphysical solutions such as parameter oscillations between layers, even if they provide better χ 2 .
An important feature of our approach is that instead of converting the data into the T MB scale for direct comparison with the non-LTE model, which is absolutely not valid for extended regions such as starless clouds, we reproduce the observations by computing the telescope coupling to the sky, based on its main beam, error beams, and forward scattering and spillover efficiency (η fss ) 4 . We then obtained our model predictions directly calibrated in the T * A scale, whereas the T MB scale is overcorrected for low η MB values such as those found for the IRAM 30 m in the 1.3-0.8 mm range since the error beams pick up a substantial fraction of the signal when the source is extended. This leads to wrong line ratios and overestimates of the temperature and/or density of the cloud and also has an impact on the species abundance.
Once the density, temperature, and kinematic properties of the onion-shell model were constrained by the N 2 H + modeling, we reproduced the spectra of N 2 D + (1-0, and 2-1), DCO + (1-0, 2-1, and 3-2), and ortho-H 2 D + (1 10 -1 11 ) by varying only their abundances and assuming that they share the same physical parameters as N 2 H + . Similarly, X(N 2 D + ), X(DCO + ), and X(ortho-H 2 D + ) profiles were obtained using the above iteratively spectral fitting procedure. Since we lacked N 2 D + observations at large radii, we set upper limits of X(N 2 D + ) in the outer layers from X(N 2 H + ) multiplied by the minimum observed X(N 2 D + )/X(N 2 H + ) = 0.05. For the non-detection observations of ortho-H 2 D + , we also estimated their abundance upper limits. Figure 8 shows the profiles of T kin , V turb , abundances of the above four cations, and the N 2 H + fractionation, X(N 2 D + )/X(N 2 H + ). The above best-fit profiles are also numerically listed in Table C.1. The error bar in each layer is determined by the range of quantity, where ∆χ 2 = χ 2 (quantity) − χ 2 (the best-fit quantity) ≤ 1. ( We only determined the upper limit of X(ortho-H 2 D + ) toward the innermost layer. A higher signal-to-noise ratio of ortho-H 2 D + would be needed to derive its lower limit because the volume of the innermost layer is the smallest and its contribution to the emergent intensity is also small. Finally, the best-fit modeled spectra of the four cations along the main horizontal cuts are shown in red in Fig. 3. The fit with the observed spectra is remarkably good. To see if the observed    (Fig. 6c), where the best-fit parameters are annotated. The profiles of kinetic temperature (T kin ), turbulent velocity (V turb ), abundances (X species ), and the abundance ratio of N 2 D + /N 2 H + are the best-fit results from the non-LTE radiative transfer calculations. spectra along the main vertical cut could also be reproduced by our onion-shell model, we simply assumed that the physical parameters and abundance profiles were the same as those of the east side except for V rot profile that is now not zero. Since the L 1512 core has a vertical motion compatible with a solid-body rotation, we adopted the V rot profile of N 2 H + , N 2 D + , and DCO + Article number, page 11 of 22 A&A proofs: manuscript no. main shown in Fig. 5c, which averages the central velocities of the vertical spectra (Fig. 5d) between the north and south sides. For ortho-H 2 D + , we simply adopted the V rot profile of N 2 H + because of its weak detection. The above V rot profiles are given at layer interfaces and the rotational velocity inside a layer is linearly interpolated between the two given velocities at the inner and outer sides of that layer. The modeled spectra of the four cations are shown in red in Fig. 4, and are again in good agreement with the observations. Our results suggest that the 1D spherically symmetric assumption is globally valid for L 1512.

Time-dependent chemical model
We adopted a pseudo time-dependent gas-phase chemical model (Pagani et al. 2009b) to reproduce the deuteration ratio, X(N 2 D + )/X(N 2 H + ), and abundances of N 2 H + , N 2 D + , DCO + , and ortho-H 2 D + in each layer. In the pseudo time-dependent approach, the physical properties (n H 2 , T kin ) are kept constant while the chemical abundances (X species ) evolve. We considered a chemical network, where (1) different spin states (ortho, para, and meta) of H 2 and H + 3 , and their isotopologs are all included, and (2) heavy molecules are totally depleted except for CO and N 2 which are partially depleted.
The initial ortho-H 2 /para-H 2 ratio is assumed to be their statistical weights of 3:1 and H + 3 is formed later via the cosmic ray ionization of H 2 . After the ortho/para ratio (OPR) of H 2 has dropped to values low enough to prevent destruction of H 2 D + by ortho-H 2 (OPR 1% or less; Pagani et al. 1992aPagani et al. , 2009b, the deuterium fractionation of H + 3 enhances in the low T kin and highly depleted environment, forming H 2 D + , and subsequently D 2 H + and D + 3 . Then H + 3 and its isotopologs transfer protons/deuterons to CO and N 2 and produce HCO + , DCO + , N 2 H + , and N 2 D + . HCO + is not a suitable tracer for starless cores because its low rotational transition emission is optically thick and this ion is present in both the envelope and the core (which is also true for H 13 CO + and HC 18 O + ), while ortho-H 2 D + , DCO + , N 2 H + , and N 2 D + emissions are optically thinner and most importantly are confined to the core (Pagani et al. 2012). Since we evaluated the quasi-instantaneous conversion of CO and N 2 into the aforementioned ions, we can assume that their abundances are in equilibrium with these ions locally. Therefore, the gasphase abundances of these two neutral molecules are free parameters in our chemical model. The other free parameters are the average grain size (a gr ) and the cosmic ray ionization rate (ζ). Figure 9 shows the chemical model solutions of X(ortho-H 2 D + ), X(N 2 D + )/X(N 2 H + ), and X(DCO + ) for each layer toward the east and west sides. Typical ζ of 10 −17 s −1 was adopted throughout the chemical model. For a gr , although Steinacker et al. (2015) find that the maximum spherical grain size can reach ∼0.5 µm in this source, the majority of the grains are still small and the larger grains have fluffy surfaces (e.g., Tazaki et al. 2016), therefore limiting the decrease of the total grain surface cross section. The impact on the chemistry and the freeze-out time should be negligible compared to the standard 0.1 µm grain size usually adopted in chemical models. In each layer, X(CO) and X(N 2 ) were adjusted to make the chemical solutions of the four cations fit the observed abundances (or remain below the observed upper limits) simultaneously. We found that the time ranges for all layers meeting the observations are 0.2-2.6 Myr for the east side and 0.2-1.9 Myr for the west side. Figure 10 shows the derived CO and N 2 abundance profiles, while their values are listed in Table C.1. Although the isotopologs of CO are not considered in our chemical network, we can obtain the X(C 18 O) profile through a 12 C 16 O/ 12 C 18 O abundance ratio in the layers in which C 18 O cannot be measured directly. In the outer region, where radiative transfer modeling is applicable, we find that the C 18 O abundance is 1 × 10 −7 up to a radius of 3 (=25,200 AU). Since X(C 18 O) is limited to 1 × 10 −7 toward the sixth eastern layer, where the 12 CO abundance reaches its maximum of 5 × 10 −5 , we find an isotopolog ratio of 12 C 16 O/ 12 C 18 O = 500. This is compatible with the local 16 O/ 18 O ISM ratio of 557±30 (Wilson & Rood 1994). As there is no reason to vary this ratio across the cloud, we obtained the inner C 18 O abundance profile between the fifth eastern and the third western layers by applying the same ratio of 500 to X( 12 CO) profile. On the other hand, the outer 12 CO abundance profile is kept as high as 5 × 10 −5 because X( 12 CO) is poorly determined from the chemical model when both X(N 2 D + )/X(N 2 H + ) and X(ortho-H 2 D + ) are only upper limits. As a result, we were able to reproduce the C 18 O spectra via the non-LTE radiative transfer calculation. Owing to the nonsymmetrical velocity shift of C 18 O, we adjusted the velocity offset of the modeled spectra individually on each spectrum based on the observations. The impact on the line shape is too small to be detectable. The northern C 18 O abundance profile (positive offsets) is different from the east-west profile and had to be decreased by a factor of 2-16 to fit the observations. This suggests a larger CO depletion toward the northern tail. It turns out that the observations can be fitted by our model results of C 18 O (1-0) and (2-1) along the main horizontal cut (Fig. 3) and main vertical cut (southward only, Fig. 4) computed with the abundances derived from the chemical model, thereby validating our chemical model results.

Discussion
As shown in Sect. 4, we built an onion-shell model of L 1512 to represent its physical structure (n H 2 , T kin , and V turb are shown in Fig. 8; V rot in Fig. 5c) and chemical abundances (N 2 H + , N 2 D + , DCO + , and ortho-H 2 D + are shown in Fig. 8; CO, C 18 O, and N 2 in Fig. 10). Our model reproduces all the observed spectra along the main horizontal and vertical cuts (Figs. 3 and 4) via the non-LTE radiative transfer calculations. In this section, we first compare our findings with other studies and then address the lifetime of L 1512.

Density and kinetic temperature
By integrating the Plummer profile to the edge of the cloud where it merges in extinction with the background (∼320 from center), we find a peak N H 2 of 2.3 × 10 22 cm −2 toward the center of L 1512, which is higher than previously reported (0.8 − 1.3 × 10 22 cm −2 ; Lee et al. 2003;Lippok et al. 2013Lippok et al. , 2016. From the non-LTE modeling of both C 18 O and N 2 H + we also find 2.1 ×10 22 cm −2 in a 180 radius, which is a lower limit since the cloud is truncated to the C 18 O emission extension and is consistent with the Plummer profile. This discrepancy is due to our different approach to derive the total column density. We use dust extinction, which is independent of the dust temperature (T dust ), and N 2 H + modeling for which the uncertainty on the collisional coefficients is extremely small compared to the uncertainty on dust properties in emission at long wavelengths. On top of dust properties uncertainties, it is well known that the problem of quantifying dust in emission is degenerate. This is because of the unconstrained temperature variation along the line of sight that can hardly be disentangled from dust property variations and the nonlinearity of the blackbody intensity with temperature for   A&A proofs: manuscript no. main the combination of wavelengths and temperatures of concern in starless cores . Indeed, Lee et al. (2003) and Lippok et al. (2013Lippok et al. ( , 2016 performed their estimations with dust thermal emission by assuming a constant T dust of 10 K or using a decreasing T dust toward the core center. However, their N H 2 are less than ∼50% of our result. Pagani et al. (2015) suggest that if only dust continuum data are used, T dust can be overestimated and consequently the column density of starless cores can be underestimated. This is because warmer dust emission can easily dominate over cold dust emission in the spectral energy distribution fitting.
For the kinetic temperature, we find a central T kin of 7.5±1 K on the west side, which is better constrained that the central temperature for the east side (7.5 +4 −2 K) . We also note that the second central layer on the east side is better constrained too (8.2 +1 −2 K). Thus, we conclude that T kin is down to 7.5±1 K in the center and rises to ∼10-17 K in the outer region. Because dust-gas coupling is efficient enough for n H 2 > 10 4.5 cm −3 (Goldsmith 2001), we can assume that T dust = T kin in the core center based on our modeling. The central T dust is therefore lower than Lippok et al. (2013Lippok et al. ( , 2016's 9.8-11.3 K, which explains our peak column density discrepancy as already noticed in other sources by Pagani et al. (2015).

Cation abundance profiles
We find that N 2 H + is significantly depleted at the center of L 1512 (n H 2 = 1.1 × 10 5 cm −3 ) and its abundance starts to decrease at n H 2 of 3 (east) to 8 (west) ×10 4 cm −3 . This is one order of magnitude lower than toward L 183. By averaging the depletion factors calculated for the east and west sides, we obtain an averaged N 2 H + depletion factor of 27 +17 −13 , which is larger than the factor of 6 +13 −3 in L 183 previously found by Pagani et al. (2007) while L 183 even has a denser central n H 2 ≥ 2×10 6 cm −3 . Strangely, Lee et al. (2003) suggests that N 2 H + may be depleted in L 1512 but not in another advanced and denser starless core, L 1544 (n H 2 ≥ 10 6 cm −3 ) whereas Redaelli et al. (2019) present new observations and analysis of L 1544 for which they report a large volumetric depletion factor (> 100). For N 2 D + , we find a depletion factor of 4 +2 −1 , which is also higher than the factor of 2-2.5 found in L 183 (Pagani et al. 2007) but lower than the one reported for L 1544 by Redaelli et al. (2019) 5 , ∼15. It is surprising that the depletion factors do not correlate with the peak density among these three sources. One of the possibilities is that the dust grain size in L 1512 might be smaller than in the other cores keeping the depletion timescale shorter (Lee et al. 2003). Another possibility is the time it took to form the cores. L 183 seems to have contracted more rapidly than the two other cores (≤0.7 Myr for L 183, Pagani et al. 2013;∼1 Myr for L 1544, Redaelli et al. 2019≥1.4 Myr for L1512, see below Sect. 5.4), though this would have given more time for grain growth in L 1512. Many coupled dynamical and chemical models also cannot simulate such high depletion factor at n H 2 10 5 cm −3 (e.g., Aikawa et al. 2001;Li et al. 2002;Pagani et al. 2013), which might be caused by their too simplistic approach. Consequently, more detailed modeling would be necessary to understand the physical and chemical evolution of L 1512.
From N 2 H + and N 2 D + , we can obtain the deuteration ratio, X(N 2 D + )/X(N 2 H + ), which is increasing from an upper limit of ≤0.05 at large radii to the maximum of 0.34 +0.24 −0.15 at the center of L 1512. The higher deuteration ratio at the center seems contradictory with the decreasing abundance of ortho-H 2 D + toward the center (by a factor of 10). It can be explained if the ortho-H 2 D + depletion is due to the further deuteration of H + 3 . Namely, H 2 D + is converted to D 2 H + and D + 3 , which is supported by our chemical modeling best-fit result, X(D + 3 )/X(ortho-H 2 D + ) ≈ 64 and X(D 2 H + )/X(ortho-H 2 D + ) ≈ 15 at the innermost layer (Fig. 11). The contrast is higher than toward L 183 (Pagani et al. 2009b(Pagani et al. , 2013 where D 2 H + is comparable to H 2 D + in abundance and D + 3 less than a factor of 10 higher than H 2 D + . Moreover, in L 183, the ortho-H 2 D + does not show depletion toward the center. This can be linked to the longer timescale of evolution of L 1512 as discussed in Pagani et al. (2013) appendix B, where it is shown that the para-D 2 H + /ortho-H 2 D + is much higher in the slow collapsing case than in the fast collapsing case (the reader should be aware that para-D 2 H + is only a small fraction of total D 2 H + when comparing the two papers).
Remaining gas-phase CO can react with deuterated H + 3 and form DCO + (Pagani et al. 2011), we find that DCO + happens to deplete at n H 2 of 3 (east) to 7 (west) ×10 4 cm −3 , and reaches a depletion factor of 9 +21 −3 at the center, which is comparable to the factor of ≥ 17 found in L 183 (Pagani et al. 2012). Like in L 183, DCO + depletes because its parent molecule, CO, freezes out onto dust grains much more than what the growth of H + 3 deuteration can compensate for (Pagani et al. 2012). Likewise, N 2 H + and N 2 D + deplete because of the freeze-out of N 2 . CO and N 2 depletion are discussed hereunder.

CO and N 2 depletion
It is impossible to measure X(CO) and X(N 2 ) directly in a dense core. This is because the emission of the low-J transitions of CO (including its isotopologs) usually becomes optically thick even for rarer isotopologs such as C 18 O and C 17 O, for which the emission is dominated by the envelope and N 2 has no transition detectable 6 in starless cores (Pagani et al. 2012).
In Sect. 4.4, we presented a chemical model to relate the four observed cations (N 2 H + , N 2 D + , DCO + , and ortho-H 2 D + ) to CO and N 2 in the depleted region to retrieve their abundance profiles. Figure 10 shows the abundance profiles of 12 CO, C 18 O, and N 2 . We find that N 2 has a depletion factor of ∼200 compared with its maximum values of 3 × 10 −5 . For CO isotopologs, 12 CO has a depletion factor of ∼220 compared with its maximum values of 5 × 10 −5 . Compared with the standard 12 CO abundance of 1-2 × 10 −4 (Pineda et al. 2010), 12 CO has a depletion factor of ∼430-870. C 18 O would also have a depletion factor of ∼220 compared with its maximum value of 1 × 10 −7 in the outer region since we fix the 12 C 16 O/ 12 C 18 O ratio of 500, and a factor 7 of ∼370 compared to the standard C 18 O abundance of 1.7×10 −7 (Frerking et al. 1982), which are much larger than the factors of 30 and 25 found by Lippok et al. (2013) 8 and Lee et al. (2003), respectively. Although both studies try to recover volumetric abundances of C 18 O and not simply column-density averages, the abundance 6 The lowest N 2 rotational transition J=2-0 at 357.808 GHz is a very weak electrical quadrupolar transition, and in any case would be blocked by the telluric opacity if it were detectable. 7 The discrepancy between the two depletion factors, 430-870 of 12 CO and 370 of C 18 O, is because the standard 12 CO abundance of 1-2 × 10 −4 and C 18 O abundance of 1.7 × 10 −7 are not compatible with a standard 12 CO/C 18 O ratio of 500. 8 However, we could not reproduce the Lippok et al. (2013) model results (the emergent intensities being typically a factor of 2 too low based on their values) and the discrepancy has not been identified, but could be linked to a possible confusion between H and H 2 densities in the Lippok et al. (2013)  isotopolog abundances derived from the chemical model constrained by the observations (in particular of ortho-H 2 D + , which explains why we separately show the two spin states for this species for the sake of clarity).
drop they derive can only be a lower limit because the line from the center is too weak to constrain the models. A factor of 30 already indicates that the contribution of the central part to the total intensity is less than ∼10%, which would require a very high signal-to-noise ratio to evaluate and still does not guarantee that the measured variation is not simply due to irregularities in the cloud shape rather than abundance variations (see Pagani et al. 2010a for a similar case).
The abundance of CO and N 2 in starless cores is a matter of long debate (see, e.g., Pagani et al. 2012;Nguyen et al. 2018) since the daughter species of N 2 , N 2 H + , is present at much higher densities than CO itself in starless cores, while CO and N 2 have identical masses and similar sticking coefficients (Bisschop et al. 2006) and their freeze-out rates are thus comparable. For L 1512, 12 CO and N 2 profiles are similar to each other within their error bars. This is different from the results reported by Pagani et al. (2012) who found a different N 2 profile in L 183, where X(N 2 ) is less than X( 12 CO) by about two orders of magnitude at the edge of the starless core. These authors suggest that N 2 might still be forming from atomic N in L 183. In that case, the depletion and production rates of N 2 would partly compensate each other. In contrast with L 183, we find that L 1512 has chemically evolved on a much longer timescale so that N 2 chemical production from atomic N has reached steady state. This could explain the similar depletion profiles. The fact that the CO depletion in L 183 is ∼3 times deeper than in L 1512 is mostly due to its peak density being 20 times higher, making collisions with dust grains much more frequent, but the difference is partly compensated by a much longer evolution time in L 1512. The L1512 case is interesting because, by showing that N 2 and CO profiles are identical while N 2 H + is tracing the core and not CO, it demonstrates that the N 2 /CO problem is not correctly formulated. N 2 H + cannot be directly compared to CO because it is a daughter species of N 2 and its abundance therefore also depends on the abundance of H + 3 . We would need to compare HCO + to N 2 H + , which is not possible in a simple way. Also, N 2 H + is a trace species, typically at the 10 −10 abundance level, at least two orders of magnitude less abundant than CO even depleted by a factor 2000 as in L 183. A correct approach is therefore to compare DCO + and N 2 D + ; in most clouds, it is clearly visible that DCO + is more abundant than N 2 D + because its emission is usually as extended as that of N 2 H + itself. There is therefore no N 2 /CO crisis. However, tracing the abundance of both species in starless cores remains important since it helps the understanding of the chemical evolution of the core. Figure 9 shows that the chemical solution of each layer meets the observed abundances within ∼0.2-2.6 Myr. Since the physical structure is kept constant in our pseudo time-dependent chemical modeling, we see that the chemical process is relatively accelerated in the inner layers compared with the outer layers. Therefore, we underestimate the timescale for the inner layers, which were not very dense at the beginning of contraction. In addition, only the deuteration ratio of N 2 D + /N 2 H + and the abundance of ortho-H 2 D + are used to constrain the chemical age, while X(DCO + ) is used to determine the CO abundance. It follows that our result gives a lower limit on the lifetime scale of L 1512 (∼1.43 Myr) that is set by the third western layer, which is the last one in our chemical model to reach the observed N 2 D + /N 2 H + ratio. As a lower limit is required, this value is set at the beginning of the time range of the third western layer as shown on the top right panel (the chemical solutions for the western deuteration ratios) in Fig. 9. We also note that since the third western layer has no ortho-H 2 D + detection but only an upper limit on its abundance, the time limit for this layer is the time it takes to reach X(N 2 D + )/X(N 2 H + ) under the assumption of X(ortho-H 2 D + ) at the maximum of its possible range, that is the fastest time of that layer. Therefore, L 1512 is presumably older than 1.43 Myr. We might ask whether the initial OPR of 3 is warranted. If we set an initial OPR as low as 0.01, we cannot find a good fit for both sides. On the other hand, if we use an OPR of 0.5 as the initial value, as measured in the diffuse medium (Crabtree et al. 2011) toward a few lines of sight, we can only find very marginal fits for both sides and the time remains somewhat above 0.5 Myr. This suggests that the OPR cannot be lower than ∼1. Pagani et al. (2009b) applied the same pseudo timedependent chemical modeling on L 183 and found the lower limit on L 183 lifetime scale to be 0.15-0.2 Myr. Later, Pagani et al. (2013) conclude that the contraction of L 183 follows the free-fall timescale (∼0.5-0.7 Myr). Namely, L 1512 is older than the L 183 starless core at least by a factor of ∼2-3. This is in agreement with our conclusion from Sect. 5.3 that L 1512 may be chemically more evolved than L 183, while physically L 183 has reached a higher density in a shorter time. Consequently, ambipolar diffusion may have slowed down the contraction of L 1512 or even halted it to the present state, while it has had no impact on L 183. This would imply that the magnetic field is stronger in L 1512 than in L 183. This will be the subject of a follow-up study of L 1512.

Lifetime scale
Article number, page 15 of 22 A&A proofs: manuscript no. main

Conclusions
We established the extinction map of L 1512 from deep JHK s and Spitzer/IRAC data. We performed a non-LTE radiative transfer calculation with an onion-shell model to reproduce the observed spectra of L1512. We obtained separately, for the east and west sides, 1D spherically symmetric profiles of the physical parameters and chemical abundances of N 2 H + , N 2 D + , DCO + , and ortho-H 2 D + . Afterward, we used a time-dependent chemical model to estimate the lower limit on the lifetime scale of L 1512 by fitting the chemical solutions to the observed abundances. Our chemical model captures the main reactions in the H + 3 deuteration fractionation and assumes that heavy species are totally depleted except for CO and N 2 , which are partially depleted. Thus, by relating the observed abundances to CO and N 2 in the depleted region, we could also obtain the CO and N 2 profiles and their depletion factors. We summarize our results as follows: 1. We derived the central molecular hydrogen density of 1.1 × 10 5 cm −3 from the dust extinction measured from NIR and MIR maps and the central kinetic temperature of 7.5±1 K from N 2 H + observations. At such density, gas is thermalized with dust. 2. The depletion factors of N 2 H + , N 2 D + , DCO + , and ortho-H 2 D + are 27 +17 −13 , 4 +2 −1 , 9 +21 −3 , and 10. Compared with N 2 H + , the smaller N 2 D + depletion factor is due to the deuterium fractionation enhanced from a upper limit of ≤ 0.05 at large radii to 0.34 +0.24 −0.15 in the center. The depletion of ortho-H 2 D + suggests that ortho-H 2 D + is further deuterated to D 2 H + and D + 3 . 3. The depletion of N 2 H + and N 2 D + traces the freezing out of N 2 toward the center of L 1512. We find that N 2 has a depletion factor of ∼200 compared to its maximum abundance of 4.7 × 10 −5 . Likewise, the depletion of DCO + indicates the freeze-out of CO. We find that 12 CO and C 18 O have a depletion factor of ∼220 between internal and external layers. If we compare their minimum abundance to their standard (literature) abundance, we find the depletion factor is about 400 for both isotopologs. 4. The similarity of CO and N 2 profiles suggests that L 1512 may have chemically evolved long enough and in particular that N 2 chemistry has reached its steady state. This could explain that L 1512 has a higher depletion factor of N 2 compared to the denser starless core, L 183. 5. L 1512 is presumably older than 1.4 Myr. We conclude that ambipolar diffusion is the dominant core formation mechanism for this source.