Gas phase Elemental abundances in Molecular cloudS (GEMS) V. Methanol in Taurus

Methanol, one of the simplest complex organic molecules in the Interstellar Medium (ISM), has been shown to be present and extended in cold environments such as starless cores. We aim at studying methanol emission across several starless cores and investigate the physical conditions at which methanol starts to be efficiently formed, as well as how the physical structure of the cores and their surrounding environment affect its distribution. Methanol and C$^{18}$O emission lines at 3 mm have been observed with the IRAM 30m telescope within the large program"Gas phase Elemental abundances in Molecular CloudS"(GEMS) towards 66 positions across 12 starless cores in the Taurus Molecular Cloud. A non-LTE radiative transfer code was used to compute the column densities in all positions. We then used state-of-the-art chemical models to reproduce our observations. We have computed N(CH$_3$OH)/N(C$^{18}$O) column density ratios for all the observed offsets, and two different behaviours can be recognised: the cores where the ratio peaks at the dust peak, and the cores where the ratio peaks with a slight offset with respect to the dust peak ($\sim $10000 AU). We suggest that the cause of this behaviour is the irradiation on the cores due to protostars nearby which accelerate energetic particles along their outflows. The chemical models, which do not take into account irradiation variations, can reproduce fairly well the overall observed column density of methanol, but cannot reproduce the two different radial profiles observed. We confirm the substantial effect of the environment onto the distribution of methanol in starless cores. We suggest that the clumpy medium generated by protostellar outflows might cause a more efficient penetration of the interstellar radiation field in the molecular cloud and have an impact on the distribution of methanol in starless cores.


Introduction
Methanol, CH 3 OH, and even more complex organic molecules (COMs) have been widely observed in star-forming regions. COMs are defined as organic molecules with ≥ 6 atoms (van Dishoeck 2009), and their formation is considered to be the first step in the chemical complexity that will eventually be inherited by forming planets in the process of star and planetary system formation. In the past decade several papers have reported on the detection of methanol and other COMs towards starless cores showing that COMs are present and their emission is extended in cold and shielded environments (Bacmann et al. 2012;Bizzocchi et al. 2014;Vastel et al. 2014;Jiménez-Serra et al. 2016;Scibelli & Shirley 2020;Scibelli et al. 2021). Starless cores provide the ingredients that will eventually build-up a planet, both the organic as well as the refractory material. Therefore, it is valuable to understand the chemistry that drives the build-up of chemical complexity in stellar nurseries such as starless cores. COMs have been widely observed towards hot cores and hot corinos, both characterised by an active chemistry driven by ice sublimation at T≥100 K (e.g. Caselli et al. (1993); Cazaux et al. (2003); Caux et al. (2011); Jørgensen et al. (2012); Belloche et al. (2014)). COMs are mostly formed on the surface of dust grain via: i) hydrogenation of carbon monoxide (producing small COMs like formaldehyde and methanol, Charnley et al. 1995;Watanabe & Kouchi 2002); and ii) radical-radical (associative) reactions (producing larger COMs like methyl formate and dimethyl ether, Garrod & Herbst 2006). In order for radical-radical COM formation to happen, it is crucial to have a central heating source like a young stellar object (YSO) because the mobility of radicals on the surface appears at T≥30 K. Because of the lack of a central heating source, the detection of COMs toward starless cores challenged our understanding of COMs in the interstellar medium. Large COMs like acetaldehyde, dimethyl ether and methyl formate have been observed towards dark cloud cores and starless cores (Bacmann et al. 2012;Vastel et al. 2014;Jiménez-Serra et al. 2016;Nagy et al. 2019;Scibelli & Shirley 2020;Scibelli et al. 2021;Jimenez-Serra et al. 2021).
The high-level of molecular complexity revealed by these studies towards cold cores has been a challenge for chemical models (Öberg et al. 2010). Recent laboratory work has shown that COMs like acetaldehyde, vinyl alcohol, and also glyicine is possible in interstellar ice analogs with non-energetic processes (Chuang et al. 2020;Ioppolo et al. 2021). The formation of COMs on the grains depends on the efficiency for the large radicals to diffuse on the surfaces. Although it is an uncertain parameter, a larger diffusion efficiency will produce more COMs on the surfaces (Ruaud et al. 2015;Walsh et al. 2016). The formation in the ices also depends on the chemical network itself, which is probably not complete for the larger molecules. Ruaud et al. (2015) for instance have proposed new routes to form COMs at low temperature (and low diffusion efficiency) via reactions induced by the formation of complexes with the main ice components. New gas-phase pathways have also been proposed but still require the efficient production and desorption of methanol (Vasyunin & Herbst 2013;Balucani et al. 2015). If the efficiency on the surface is sufficient, remains the problem of desorbing these molecules strongly bound to the surfaces. The most promising non-thermal desorption mechanisms seem to be chemical reactive desorption from CO-rich surfaces (Minissale et al. 2016;Vasyunin et al. 2017), chemical explosions (Rawlings et al. 2013;Ivlev et al. 2015), and grain sputtering induced by cosmic-rays (Wakelam et al. 2021). Methanol is one of the simplest COMs and the precursor of many of the O-bearing COMs, hence understanding the processes responsible for its release in the gas phase is of pivotal importance to constrain chemical models and reproduce the chemical complexity that we observe in starless cores.
In this paper we present our results on the emission of methanol towards 12 starless cores as well as towards the clouds where they are embedded. We study the dependency, or the lack thereof, between the emission of methanol and physical parameters such as visual extinction, volume density, as well as environment. The paper is structured as follows: Section 2 describes the observations, Section 3 presents the dataset as well as the TMC-1 molecular cloud and the B213 filament, where the 12 observed starless cores are located. The analysis of the observed data and its results are in Sections 4 and 5. The results of the chemical modelling are presented in Section 6.2. Our conclusions are discussed in Section 7. More detailed information on each source is presented in the Appendice.

Observations
The molecular transitions used in this study are reported in Table 1. The 3 mm methanol lines were observed within the "Gas phase Elemental abundances in Molecular CloudS" (GEMS) (PI: Asunción Fuente) IRAM 30m Large Program and the observing procedure and receiver setups were described by Fuente et al. (2019). The HPBW is varying with the frequency as HPBW(")=2460 /ν where ν is in GHz. The observing mode was frequency switching with a frequency throw of 6 MHz well adapted to remove standing waves between the secondary mirror and the receivers. The Eight MIxer Receivers (EMIR) and the Fast Fourier Transform Spectrometers (FTS) with a spectral resolution of 49 kHz were used for these observations. The intensity scale is T M B , which is related with T Table B.1 in Fuente et al. 2019). GEMS observations in TMC 1 were completed with C 18 O 2→1 observations carried out with the IRAM 30m telescope in a previous project.
For TMC 1, we use observations of the CH 3 OH J Ka,Kc = 1 0,1 →0 0,0 line carried out with the Yebes 40m radiotelescope in Q-Band (Tercero et al. 2020) during March-April 2018. The 40m telescope is equipped with HEMT receivers for the 2.2-50 GHz range, and a SIS receiver for the 85-116 GHz range. Single-dish observations in K-band (21-25 GHz) and Q-band (41-50 GHz) were performed simultaneously. The backends consisted of FFTS covering a bandwidth of ∼2 GHz in band K and ∼9 GHz in band Q, with a spectral resolution of ∼38 kHz. Central frequencies were 23000 MHz and 44750 MHz for the K and Q band receivers, respectively. The observing procedure was position-switching, and the OFF-positions are RA(J2000) = 04 h 42 m 24 s .24 Dec (2000):25 • 41 27 .6 for TMC 1-CP, RA(J2000) = 04 h 42 m 29 s .52 Dec (2000):25 • 48 07 .2 for TMC 1-NH 3 , RA(J2000) = 04 h 42 m 32 s .16 Dec(J2000):25 • 59 42 .0 for TMC 1-C. These positions were checked to be empty of emission before the observations. The intensity scale is T M B with conversion factors of 4.1 Jy/K in band K (T M B /T * A =1.3) and in 5.7 Jy/K in band Q (T M B /T * A =2.1). The HPBW of the telescope is 42" at 7 mm and 84" at 1.3 cm. As an example of our dataset, Figure 1 shows to the left the H 2 column density map of the TMC-1 molecular cloud (upper panel) and the B213 filament (lower panel) computed from Herschel and P lanck data (Palmeirim et al. 2013;Rodríguez-Baras et al. 2021), and on the right a zoom-in into the H 2 column density map of B213-C1 where the observed offsets are marked (upper panel) and the spectra of the 2 1,2 -1 1,1 (E 2 ) transition of methanol are overlaid with the 1-0 transition of C 18 O (lower panel). The spectra observed towards the other sources in our sample are shown in Figures A.1 Notes. (a) Frequencies and uncertainties from the CDMS (Müller et al. 2005) (b) n * is the critical density, calculated at 10 K.
(c) Energy relative to the ground 00,0, A rotational state. (d) Calculated with the rate coefficient reported in Bizzocchi et al. (2014)

Dataset
We have observed several lines of CH 3 OH and two lines of C 18 O towards the GEMS sources in Taurus with the aim of understanding how these two chemically related molecules behave towards several starless cores within the same environment. The GEMS dataset allows us to study the abundance of both molecules, as well as their ratios, towards different offsets within and outside the cores, and hence tracing gas with different physical properties such as volume density and temperature. Our source sample is listed in Table 2. All transitions of methanol studied in this paper, except the J Ka,Kc = 1 0,1 -0 0,0 (A + ) transition at 48 GHz, have been observed towards every offset in our source sample (although they were not always detected). Both C 18 O transitions have been observed towards every offsets in the starless cores in TMC-1, while only the J = 1-0 transition has been observed towards the starless cores in the B213 filament. Figures A.1-A.11 show the spectra of the 2 1,2 -1 1,1 (E 2 ) transition of methanol overlaid with the 1-0 transition of C 18 O towards all offsets in our dataset, with the exception of B213-C1, shown in Figure 1.

TMC-1 cloud
Taurus Molecular Cloud-1 (TMC-1) is a cold and dense cloud with a filamentary structure of about 5 ×40 in the centre of the Taurus molecular cloud at a distance of 140.2 +1.3 −1.3 pc (Galli et al. 2019). The chemistry of TMC-1 has been the object of many studies in the literature (e.g. Guelin et al. 1982;Pratap et al. 1997;Saito et al. 2002;Schnee et al. 2007;Suutarinen et al. 2011). The most peculiar feature of TMC-1 is the chemical differentiation between the southern part, the cyanopolyyne peak (TMC-1 CP), where carbon-chain molecules are more abundant, and the northern part, the ammonia peak (TMC-1 NH 3 ), where NH 3 and N 2 H + are more abundant. The origin of this differentiation has been extensively discussed, and it can be explained by either a small variation in the chemical timescale (∼ 10 5 yr), and/or by the different density (with the ammonia peak being denser) (Hirahara et al. 1992;Suzuki et al. 1992). Also a variation in the C/O ratio can reproduce an abundance gradient similar to that observed (Pratap et al. 1997). We present here the observations of methanol and C 18 O towards the cyanopolyyne and ammonia peaks (TMC-1 CP and NH 3 ). We also present observations towards the pre-stellar core TMC1-C (Schnee et al. 2007(Schnee et al. , 2010. Several observational evidences suggest that this core is a later evolutionary stage than TMC1 (CP), presenting higher deuteration fractions, more similar to those in L1544 (Crapsi et al. 2005;Navarro-Almaida et al. 2021). The spectra of the 2 1,2 -1 1,1 (E 2 ) transition of methanol overlaid with the 1-0 transition of C 18 O towards TMC-1C, TMC-1 CP and TMC-1 NH 3 are shown in Figures A.1-A.3.

B213 filament
B213 (L1495) is one of the longest and most prominent filaments in the Taurus molecular cloud, as it extends for over 80 , and it is located towards the north-west of the cloud at a distance of 129.9 +0.4 −0.3 pc (Galli et al. 2019). The gas towards B213 presents a very complex velocity structure that has been interpreted as the result of colliding filaments (Duvert et al. 1986). B213 has been observed in CO isotopologues with large scale maps (Onishi et al. 1996;Goldsmith et al. 2008). The denser gas has been studied using ammonia, H 13 CO + , N 2 H + and SO (Benson & Myers 1989;Hacar et al. 2013). Furthermore, some dense cores in B213 contain young stellar objects (Rebull et al. 2010). We present here the observations of methanol and C 18 O towards nine starless cores in B213. The spectra of the 2 1,2 -1 1,1 (E 2 ) transition of methanol overlaid with the 1-0 transition of C 18 O towards the starless cores we observed in B213 are shown in Figure 1 and Figures

Analysis
We have observed the three rotational transitions belonging to the 2-1 transitions of methanol at 96 GHz (symmetry A, E1 and E2) as well as the 1-0 transition at 108 GHz (symmetry E2-E1), in all cores of our sample. In addition, we have observed the J Ka,Kc = 1 0,1 -0 0,0 (A + ) transition of methanol at 48 GHz towards the three cores in TMC-1. The intensity ratios among the four lines at 3 mm in our dataset cannot be reproduced assuming LTE, implying that we have to use a non-LTE code if we want to derive a reliable column density of methanol in all positions within our dataset. Since methanol lines have a relatively high critical density (between 0.7 and 3 ×10 4 cm −3 , see Table 1), the correspondent column density is very sensitive to variations in the H 2 volume density. The H 2 volume density towards all observed positions within the GEMS catalogue has been calculated by modelling the emission of CS, C 34 S, and 13 CS with the radiative transfer code RADEX (van der Tak et al. 2007) in Rodríguez-Baras et al. (2021). However, the lines of methanol and CS could come from different velocity components and/or different layers within the same cloud. The line profiles of C 18 O and CH 3 OH are different in many positions. Therefore, in order to be consistent and derive accurate column densities we need to use the volume density derived from CH 3 OH data. We hence decided to use the Markov Chain Monte Carlo method (MCMC) together with the RADEX non-LTE radiative transfer code within the CASSIS software (Vastel et al. 2015) in order to compute the methanol column density, the excitation temperature, and the gas volume density towards all positions in our dataset. To assess the goodness of the fit, the reduced χ 2 was computed for each model. The reduced χ 2 in our sample ranges from 1.5 to 4.5, with most of the models having a reduced χ 2 ∼ 2. The results of the MCMC+RADEX analysis with CASSIS on methanol are reported in Table 3. We used the collisional rates for methanol reported in Rabli & Flower (2010) at 10 K.
We did not use the MCMC method for the C 18 O because for most of the cores we observed only one transition (with the exception of the three cores in TMC-1), and the observed lines of C 18 O are not very sensitive to variations in the volume densities. Furthermore, C 18 O often has more velocity components than CH 3 OH, and we are interested in deriving a column density for C 18 O from the velocity component that corresponds to the one observed with CH 3 OH. We used two approaches to derive the column densities of C 18 O in our dataset. In the offsets with a good match between the rest velocity and line width of the lines of methanol and C 18 O, we fit the lines of C 18 O and used RADEX to derive the column density, assuming the T kin and n H2 derived for CH 3 OH. In these cases, the errors derived on the v LSR and δv from the line fit have been propagated to derive the error on the column density of C 18 O, see for example the values for TMC-1C and B213-C6 reported in Table 3. For the offsets with a bad match between the rest velocity and line width of the lines of methanol and C 18 O (see for example the spectra of TMC-1 NH 3 in Figure A.3), we have used the GUI of CASSIS and performed a single-line fit with RADEX of the portion of the spectrum of C 18 O that is comparable with methanol in velocity and line-width. In order to do this we have fixed the v LSR and adjusted the line-width in order to be similar to the line-width of methanol, and at the same time model well the line shape of C 18 O. The error on the column densities derived in this manner, and reported in Table 3, has been estimated to be 15%. While in the case of CH 3 OH our analysis shows that the molecules are not in LTE (as the resulting T ex is lower than the kinetic temperature), this is not the case for the C 18 O 1-0 transition, where T ex ∼ T kin .
The visual extinction tabulated in Table 3 has been calculated from the H 2 column density maps reported in Rodríguez- Baras et al. (2021) (Bohlin et al. 1978).
We do not report on the upper limits for the column density of methanol towards positions where the line was not detected because it is not possible to derive them in the same fashion as we did for the methanol column densities reported in Table 3, and we would have to assume a value for the H 2 volume density. The uncertainty induced by this assumption would make the comparison of the upper limits on N(CH 3 OH) with the N(CH 3 OH) computed with the MCMC+RADEX method not meaningful.  Fixed to a value that matches best the δv of CH3OH and the line shape of C 18 O. vLSR and δv are reported without uncertainty because the statistical uncertainty is much lower than the spectral resolution. Figure 2 shows the column density ratio of methanol and C 18 O in all observed positions in our sample plotted against the visual extinction and the dust temperature derived from Herschel and P lanck data (Palmeirim et al. 2013;Rodríguez-Baras et al. 2021), as well as against the volume density derived from the methanol MCMC+RADEX analysis, presented here in Section 4. Figure 2 shows no correlation among N(CH 3 OH)/N(C 18 O) and A V , T dust , or n H2 . This also holds if we consider B213 and TMC-1 separately or if we consider solely the column density of methanol instead of the column density ratio with respect to C 18 O. In Figure 2, the points belonging to B213 and TMC-1 are shown as full circles and as crosses, respectively. In the case of multiple velocity components reported for methanol and C 18 O in Table 3, we assumed that most of the dust emission comes from the fiber with a larger methanol column density (supposedly the densest). As a consequence, only the velocity component with the largest methanol column density is plotted against A V and T dust in Figure 2. In particular, in the case of TMC-1 CP and NH 3 , we plot the velocity component at 5.7 km/s. As we aim at studying methanol in different positions across different cores, comparing solely the column density of methanol would give us a biased picture because of the different densities of the positions involved. On the other hand, comparing methanol abundances with respect to molecular hydrogen is also not an optimal solution given the presence of different velocity components in several spectra, corresponding to different fibers on the line of sight. We decided to use the column density ratio of methanol and carbon monoxide because they are chemically related, as the depletion of carbon monoxide on the surface of dust grains is necessary to the formation of methanol. Using the rare isotopologue C 18 O we can limit, if not exclude, the risk of large optical depth. A possible reason for the lack of definite trends in Figure 2 might be that methanol emits from the outer layers of starless cores (e.g. in L1544 as shown in Bizzocchi et al. 2014;Vastel et al. 2014). Furthermore, the emission of methanol is not homogeneous around pre-stellar and starless cores, and its distribution has been shown to depend on the density structure around the core, and hence the illumination onto the core due to the interstellar radiation field (Spezzano et al. 2016(Spezzano et al. , 2020. In Spezzano et al. (2016), for example, the emission maps of methanol and cyclopropenylidene towards the inner 2.5"×2.5" of the pre-stellar core L1544 show that both molecules have asymmetric distribution around the core, with methanol peaking towards the North-East and cyclopropenylidene towards the South-West. It is important to note that the visual extinction at the methanol and the c-C 3 H 2 emission peak in L1544 is the same, within error-bars. This shows that the H 2 column density of the medium computed on the line-of-sight is not a good enough indicator to discriminate among the regions where the illumination is effective enough to keep more Carbon in its atomic form and hence available to form hydrocarbons like c-C 3 H 2 , and the regions where there is enough shielding to allow carbon to be locked in to carbon monoxide, and eventually form methanol. In addition to the density structure around the core, also larger scale environmental effects have shown to influence the distribution of methanol around starless and pre-stellar cores. In Spezzano et al. (2020), methanol was mapped towards 4 pre-stellar and 2 starless cores, and it was shown that the methanol peak was influenced by the presence of nearby massive stars in the two cores mapped in Ophiucus (OphD and HMM-1). The GEMS dataset offers the opportunity to study the emission of methanol towards 12 starless core, a larger sample with respect to previous studies. Furthermore, the GEMS dataset is not limited to the cores, giving us information on the transition from core to cloud. In order to spot trends that might be related to the density structure surrounding the core, in Figures A.12 to A.23 we plot the column density ratios ordered by offset position across the observed cut within the core for each source, while the offsets are shown on top of the N(H 2 ) column density maps of the cores in the upper panel of the Figure to facilitate the comparison. The individual results are discussed in detail in Appendix A. Generally, we can recognise two different behaviours, shown in Figure 3: the cores where N(CH 3 OH)/N(C 18 O) peaks at the dust peak, and the cores where the N(CH 3 OH)/N(C 18 O) ratios peak has a slight offset with respect to the dust peak (∼10000 AU). This trend is particularly clear in the cores where more than 5 offsets have been observed. B213-C1 and C5 for example show a clear peak of the column density ratio towards the dust peak (see Figures 3, A.15, and A.17), while in B213-C10 and C16 the peak is slightly shifted from the dust peak (∼10000 au), (see Figures 3, A.20 and A.22). One reason for this behaviour might be the external irradiation on the cores. B213-C1 and C5 are in fact located towards the northern part of B213, that is known to be actively forming stars with ∼40 Class 0/I/II sources present in the vicinity of B213-C1, C5, C2, C6 and C7. B213-C10 and C16, on the other side, are located in the central/southern part of B213, which only contains ∼15 Class 0/I/II sources around B213-C10, C12 and C16, and 5 around B213-C17. In Figure A.25 the positions of the young stellar objects in B213 are marked as red triangles on the H 2 column density map of the filament. Low-mass protostars can accelerate energetic particles along the outflow, and consequently increase the local cosmic rays flux (Padovani et al. 2016;Fitz Axen et al. 2021). This will in turn increase the destruction of CO, and consequently decrease the production of methanol towards the outer layers of starless cores in the North of B213.

Results
The line-width of the methanol lines in the cores located in the North of B213 are ∼25% broader than the lines of methanol observed towards the central and southern cores of B213, hinting at more prominent effects due to stellar activity in the northern part for B213. This effect has already been observed with ammonia in Seo et al. (2015). An additional indication of the environmental effects is the fact that the dust temperature in the southern cores of B213 is lower with respect to the northern cores, at comparable visual extinctions, see Figure A.24. In addition to the irradiation, there are also differences in the velocity structure. In B213-C1, C2, and C7, the profiles of C 18 O and CH 3 OH are very different. In B213-C6, C10 and C16 we have however a better agreement between the profiles of C 18 O and CH 3 OH.
The N(CH 3 OH)/N(C 18 O) column density ratios towards the three cores in TMC-1 can be only calculated towards few offsets close to the dust peak, because methanol was not detected at larger radii, in contrast with B213, likely because the cores in B213 are denser. Towards TMC-1 NH 3 the two different velocities show peaks of the N(CH 3 OH)/N(C 18 O) ratio at different offsets, maybe hinting at the fact that the gas they are tracing does not match in chemical age. Towards TMC-1 CP, if we exclude the offset 4 for the component at 5.6 km/s because of the large error bar, both velocity components show the same trend, an increase of the N(CH 3 OH)/N(C 18 O) ratio towards the dust peak. Towards TMC-1 C instead, the N(CH 3 OH)/N(C 18 O) ratio is constant within error-bars.

Physical structure for spherical 1D static chemical model
For a proper chemical modeling of the region under study, it is necessary to have a good knowledge of the density and temperature across the region. However, sometimes our knowledge is limited to a few points, and geometrical considerations and simplifications are mandatory. Temperature parameterizations of prestellar cores can be found in, for example, Crapsi et al. (2007), in which they choose the spherically symmetric profile where T out and T in are the temperatures at the outermost and innermost observed positions, respectively. The dependence of the temperature with the radius is controlled by r T 0 , the flat radius of the temperature, and α, an asymptotic power index. A similar functional dependence with the radius is found in Plummer-like density profiles. These profiles assume spherical symmetry, and they have become popular density profiles to describe dense cores (Tafalla et al. 2002;Priestley et al. 2018): where n H is the total atomic hydrogen number density across the profile, n 0 is the central density, r 0 is the flat radius, and α is the asymptotic power index. These profiles allow us to obtain density and temperature profiles across the B213 cuts in Table 3. To have more data to derive physical profiles, we may fit data of different cuts with similar temperatures simultaneously. Therefore, we merge and average the data from the cut B213-C2 with B213-C6 (B213-C2-C6), B213-C10 with B213-C12 (B213-C10-C12), and consider B213-C16 as independent from the rest. The temperature profiles of B213-C2-C6, B213-C10-C12, and B213-C16 can be determined by finding the set of parameters that best fit the experimental data from Table 3. The results are shown in Table 4.
To perform the fitting of the parameters to the extinction values, we have subtracted the background extinction, that is, the extinction at the outermost positions, to each cut. The results of the fitting process are shown in Table 4.

Chemical modelling approaches
Two approaches were used to reproduce our observations with chemical models. In the first, we used the physical structure of some of the B213 cores in our sample derived in Section 6.1, and ran a spherical 1D static chemical model to compute the radial abundance profiles of methanol and CO. In the second approach, applied to the low-density parts of the filament beyond the starless cores, we ran a 0D chemical model assuming a grid of temperature and A V to compute the abundances of methanol.
To model the emission of methanol in cold starless cores, it is essential to utilize a so-called three-phase astrochemical model. By three phases, we mean gas phase, surface of icy mantles of interstellar grains, and bulk ice of grain mantles. The usage of a three-phase model, a more advanced approach in comparison to two-phase (gas phase + ice without distinction between surface and bulk), is justified by the importance of the details of chemical processes on interstellar grains for the formation of methanol and its abundance in the gas phase. Since Geppert et al. 2006, the consensus is that observed abundances of interstellar methanol cannot be formed in gas-phase chemical reactions. In contrast, formation of methanol during hydrogenation of CO molecules on interstellar grains is shown to be very efficient (e.g. Watanabe & Kouchi 2002). The delivery of methanol formed on cold (∼10 K) grains to the gas phase is most likely due to the process of reactive desorption (Garrod et al. 2007;Minissale et al. 2016). The efficiency of reactive desorption depends on the details of a particular chemically reacting system and the composition of the environment, i.e., the underlying ice surface (Minissale et al. 2016;Chuang et al. 2018). Thus, in order to model properly a complicated relation between CO and CH 3 OH in prestellar clouds, we need to utilize the reasonably detailed approach to grain chemistry and gas-grain interaction. To model the emission of methanol, we utilized the MONACO model previously applied to studies of chemistry in a number of prestellar cores (Vasyunin & Herbst 2013;Vasyunin et al. 2017;Nagy et al. 2019;Lattanzi et al. 2020;Harju et al. 2020;Scibelli et al. 2021;Jimenez-Serra et al. 2021). MONACO is a rate equations-based numerical code capable of simulation chemistry in interstellar medium under a three-phase approach. As described in details in Vasyunin et al. (2017), the code includes treatment of chemical reactions in the gas phase, as well as diffusive chemistry on surfaces of icy mantles on interstellar grains, and in the bulk of icy mantles. The code includes treatment of thermal desorption, cosmic ray-induced desorption , photodesorption with photodesorption yield equal to 10 −3 particles per incident photon (Öberg et al. 2008) and reactive desorption following parametrization by Minissale et al. (2016). The efficiency of reactive desorption is further adjusted to the fraction of ice surface covered with non-water species (Vasyunin et al. 2017). We assume low diffusion-to-desorption energy ratio E dif f /E des = 0.3 for thermal diffusion of species on ice surface. In addition, quantum tunneling as a source of mobility for atomic and molecular hydrogen is assumed through a rectangular barrier for diffusion with thickness of 1.2 angstroms (Vasyunin et al. 2017). For diffusion of bulk species, E dif f /E des is assumed twice that of for surface, i.e., 0.6 following Garrod (2013). The results of the spherical 1D static chemical models are presented in Figures 4 and 5. More details on the physical profiles used for the cores can be found in Section 6.1. To convert abundances of species per unit volume to directly observed values of column densities, we assumed that modeled starless cores are spherically symmetric with radial density profiles derived as described in Section 6.1. Convolution with IRAM 30m beam has also been taken into account. The procedure is described in more details in Jiménez-Serra et al. (2016).
To explore the abundances of CO and CH 3 OH in the molecular gas that surrounds starless cores in the filament, we utilized the same chemical model as for the cores. However in contrast to starless cores, the surrounding gas does not exhibit a well-shaped spatial structure. Thus, we utilized simple 0D chemical models with constant values of temperature and density. Since both parameters can vary in the surrounding gas, we calculated grid of 0D models on a parameter space that we believe reasonable covers the ranges of temperature and density in the filament gas. The results of such gridded modeling allows us to explore the behavior of CO and CH 3 OH abundances beyond the starless cores in B213 filament. The results of the 0D models are presented in Figure 6. Figure 4 compares the N(CH 3 OH) derived from the spherical 1D static chemical model used with the physical structures derived for B213 C2-C6, C10-12 and C16 at 10 5 , 10 6 and 10 7 years, with the column densities derived in this paper for B213-C6, C10 and C16 and reported in Table 3. The models at 10 5 yr reproduces fairly well (within a factor of a few) the column densities of methanol. The model using the physical structure derived for the cores C2 and C5 (left panel in Figure 4) reproduces very well also the radial profile of observed column densities of methanol towards B213-C5. We decided to plot only the N(CH 3 OH) computed for C5 and not for C2 because the latter has larger error bars. The radial profiles observed for C10 and C16 (central and right panels in Figure 4) are instead not very well reproduced by the models. While it is important to keep in mind the intrinsic large uncertainties of the column densities computed by chemical models (Vasyunin et al. 2004), we observe that the changes in N(CH 3 OH) as a function of radii cannot be reproduced by the chemical models for C10 and C16, while it can reproduced fairly well for C6. This is shown even if we look at the radial profiles of methanol and C 18 O in Figures A.12 to A.23, with respect to the one predicted by the models in Figure 5. It is important to note that the physical structures used to run the spherical 1D static chemical model have been derived from Herschel observations, which are not sensitive to the central and densest part of the cores because of the large beam (∼40"). The volume density in starless cores could in fact have a steep increase towards the inner 40", especially in the more evolved pre-stellar cores (see for example Figure 5 in Crapsi et al. 2007). Figure 5 shows the radial distribution of the N(CH 3 OH)/N(C 18 O) column density ratio, as well as the column densities of CH 3 OH and C 18 O (calculated from the main isotopologue assuming a 16 O/ 18 O ratio of 500, Solar value in Penzias 1981) in order to facilitate the comparison with our observational results shown in Figures A.12 to A.23. Independently from the physical structure used, the model predicts a drop of CH 3 OH column density of about one order of magnitude from 10 5 to 10 6 years, and an additional, although less significant, drop from 10 6 to 10 7 years. The column density of C 18 O is underpredicted by the model. However, the model is only considering the core, and not also the cloud where the core is embedded. To check if the difference in C 18 O column densities between model and observation is only due to the lack of the contribution from the cloud, we have calculated the C 18 O column density at offsets outside of the cores, for example offset 7 in B213-C2, and then added it to the column densities predicted by the model, and we saw that the contribution of the cloud can account for the missing CO.

Modelling results
In some of our cores the observed cuts cover over 40000 AU, hence merging into the cloud, while the physical models of the cores are derived assuming a Bonnor-Ebert sphere with a radius of 10000 AU. We therefore computed a grid of models in the temperature range of 8-14 K and A V range of 2-8 mag. Please note in this case the A V is the local A V , so it needs to be multiplied by a factor of two to be compared with the A V tabulated in Table 3. The column density of methanol resulting from the grid models at 10 5 , 10 6 and 10 7 yr are shown in Figure 6. Also in this case, the models cannot reproduce the variety of profiles that we observe. They predict a substantial drop of the methanol in gas-phase between 10 5 and 10 6 yr because of freeze-out, as well as a defined peak of the methanol column density, independent from the dust temperature, at A V ∼ 4 mag (corresponding to a total visual extinction of ∼8 mag along the line of sight). This is in contrast with what we observe, i.e. that the cores where the dust temperature is higher, tend to have the methanol peaking at larger A V . It is important however to note that protostellar feedback (e.g. higher cosmic-ray ionization rate and/or shocks) is not included in the models.

Conclusions
Our observations of methanol with the GEMS dataset show that methanol is present and its emission is extended in cold starless cores and in their surrounding clouds. In our dataset in fact, methanol has been detected at visual extinction ranging from 4 to 27 mag. Furthermore, the different behaviours in our sample of starless cores in Taurus show that local environmental differences are significative and have an impact on the distribution of methanol.
Previous work on emission maps of methanol towards 4 starless and 2 pre-stellar cores suggested that the asymmetric distribution of methanol is linked to the large scale density structure around the cores as well as different amount of illumination from the external radiation field (Spezzano et al. 2016(Spezzano et al. , 2020. With this work we can confirm that environmental effects have a strong impact on the methanol distribution. However, we show that the column density and temperature of the medium computed on the line-of-sight are not good enough indicators to identify the effects of the environment on the methanol, and that the level of star formation activity needs to be taken into account. It is clear from our data that low-mass protostars in the surrounding environment have an impact on the distribution of methanol towards starless cores. Low-mass protostars in fact can accelerate energetic particles and increase the local cosmic-ray flux and ionisation rate. This can locally increase the abundance of He + , the main destruction partner of CO in dark clouds, thus reducing its abundance in the gas and solid form (the latter being the first step toward CH 3 OH production). Increased fluxes of cosmic rays can also increase the CO desorption rate (e.g. , again reducing the amount of available solid CO for its transformation into methanol via successive hydrogenation. The higher CO desorption rate could compensate the gas phase destruction of CO via He + ; indeed we do not see any local reduction of the C 18 O column density around the cores in the North. Dedicated observations targeting molecular ions are needed to map the variation of ionization rate across low-mass star forming regions in order to test the effect on the chemistry with the chemical models. Another consequence of the presence of low-mass stars nearby is that they contribute to the clumpiness of the medium, escavated by the cones of the outflows, which could lead to a more efficient penetration of the interstellar radiation field in the molecular cloud. In these conditions, only the densest regions (i.e. the centers of the starless cores) will be screened enough from UV photons to allow CO to accumulate on the surface of dust grains. The higher dust temperature measured by Herschel in the northern region of B213 could indeed hint to this scenario, but higher resolution mapping of the dust continuum emission and gas tracers are needed to make quantitative conclusions.
The comparison with state-of-the-art chemical models shows that overall the agreement of the modelled column densities with the observed ones is quite good, and that reactive desorption is a very efficient method to release the methanol formed on the surface of dust grains into the gas-phase, reproducing the observations. The discrepancies among models and observations in the reproduction of the radial profiles might be due to incomplete formation/destruction/desorption mechanisms in the chemical models, to the fact that the protostellar feedback is not included in the models, as well as to the intrinsic limit of using spherical models to reproduce an asymmetric emission. Additional experimental and theoretical work on, for example, the chemical desorption and cosmic ray sputtering on different ice mixtures would be helpful for a better match between observations and models. Collisional rates for methanol with molecular hydrogen at temperatures lower than 10 K would also be beneficial. We have observed CH 3 OH and C 18 O (both lines) towards four positions in TMC-1 C. The J Ka,Kc = 2 1,2 -1 1,1 (E 2 ), 2 0,2 -1 0,1 (A + ), 2 0,2 -1 0,1 (E 1 ), and 2 0,2 -1 0,1 (E 1 -E 2 ) transitions of methanol have been observed towards all four positions, while the 1 0,1 -0 0,0 (A + ) transition has been observed only towards the offsets 3 and 4. All lines have been detected with a signal to noise larger than 3 towards all positions. The spectra of the 2 1,2 -1 1,1 (E 2 ) transition of methanol overlaid with the J = 1-0 transition of C 18 O are shown in Figure A.1, and they show a very good match both in v LSR and in line-width. Figure A.12 shows in the upper panel the observed offsets within TMC-1 C on the H 2 column density map derived from Herschel and P lanck data (Rodríguez-Baras et al. 2021) data, and the variation of the CH 3 OH and C 18 O column density ratio, as well as the single column densities N(CH 3 OH) and N(C 18 O) in the observed cut across TMC-1 C in the lower panel. Both N(CH 3 OH) and N(C 18 O) increase when going towards dust peak of TMC-1 C. The column density ratio is constant, within error-bars, across the observed cut.
Appendix A.2: TMC-1 CP We have observed CH 3 OH and C 18 O (both lines) towards four positions in TMC-1 CP. The J Ka,Kc = 2 1,2 -1 1,1 (E 2 ), 2 0,2 -1 0,1 (A + ), 2 0,2 -1 0,1 (E 1 ), and 2 0,2 -1 0,1 (E 1 -E 2 ) transitions of methanol have been observed towards all four positions, while the 1 0,1 -0 0,0 (A + ) transition has been observed only towards the offset 4. All lines have been detected with a signal to noise larger than 3 towards all positions. The spectra of the 2 1,2 -1 1,1 (E 2 ) transition of methanol overlaid with the J = 1-0 transition of C 18 O are shown in Figure A.2. TMC-1 CP is known to have multiple velocity components (Fuente 2019 and references therein). From the spectra shown in Figure A.2 we see that methanol is tracing two of the many velocity components seen in C 18 O. Two velocity components are present towards TMC-1 CP, one at 5.6 km/s and one at 6.0 km/s, with the latter having column densities of both CH 3 OH and C 18 O larger by about a factor of 2 with respect to the lower velocity component. Figure A.13 shows the observed offsets on the H 2 column density map in the upper panel, and the variation of the CH 3 OH and C 18 O column density ratio, as well as the single column densities N(CH 3 OH) and N(C 18 O) in the observed cut across the core for both velocity components, in the lower panel. In the lower velocity component the column density of both methanol and C 18 O increase towards the dust peak. Excluding the offset 4, that has a larger error in comparison with the other positions, also the column density ratio increases towards the dust peak.
Appendix A.3: TMC-1 NH 3 We have observed CH 3 OH and C 18 O (both lines) towards four positions in TMC-1 NH 3 . The J Ka,Kc = 2 1,2 -1 1,1 (E 2 ), 2 0,2 -1 0,1 (A + ), 2 0,2 -1 0,1 (E 1 ), and 2 0,2 -1 0,1 (E 1 -E 2 ) transitions of methanol have been observed towards all four positions, while the 1 0,1 -0 0,0 (A + ) transition has been observed only towards the offset 4. All lines have been detected with a signal to noise larger than 3 towards all positions. The spectra of the 2 1,2 -1 1,1 (E 2 ) transition of methanol overlaid with the J = 1-0 transition of C 18 O are shown in Figure A.3. The velocity structure that we can see in the C 18 O lines is quite complex, and the mismatch with the peaks of the methanol lines suggests the presence of more velocity components that we can actually resolve, at least in C 18 O. TMC-1 NH 3 is the ammonia peak in TMC-1 and supposed to be the most evolved core, with respect to TMC-1 CP for example (Hirahara et al. 1992). Two velocity components are present for methanol, while a more complex structure is present in C 18 O, see Figure A.3. Figure A.14 shows in the upper panel the observed offsets within TMC-1 NH 3 on the H 2 column density map derived from Herschel and P lanck data (Rodríguez- Baras et al. 2021), and the variation of the CH 3 OH and C 18 O column density ratio, as well as the single column densities N(CH 3 OH) and N(C 18 O) in the observed cut across the core for both velocity components, in the lower panel. In the component at 5.7 km/s, the column densities of both molecules have a peak at the offset 3, and then decreases moving towards the dust peak. The column density ratio shows the same behaviour. In the component at 6.0 km/s instead, while the column density of C 18 O decreases towards the dust peak, the methanol column density has a peak at the offset 2, and so does the column density ratio.

Appendix A.4: B213 C1
We have observed all lines of CH 3 OH reported in Table 1, with the exception of the 1 0,1 -0 0,0 (A + ) transition, and the 1-0 transition of C 18 O towards nine positions in B213-C1. The observed lines have been detected with a signal to noise larger than 3 towards all positions. The spectra of the 2 1,2 -1 1,1 (E 2 ) transition of methanol overlaid with the J = 1-0 transition of C 18 O are shown in Figure 1. The velocity structure is quite complex, with both molecules showing multiple velocity components. A shift in the v LSR of the main velocity component across the observed cut, is also clearly seen. The upper panel in Figure A.15 shows the H 2 column density map of B213-C1 as well as the offsets of our observations. The observed stripe is centred at the denser peak, labeled as offset 1, but also passes through a less dense tail towards the East (offsets 6 and 7). In the lower panel of Figure A.15 are shown both the column density ratio of CH 3 OH/C 18 O and the single column densities across the observed stripe, which all decrease rather sharply towards the East after the offset 3 and peak towards the dust peak.   (Palmeirim et al. 2013). The young stellar objects reported in Rebull et al. (2010) are shown as triangles.