The Perseus ALMA Chemistry Survey (PEACHES). I. The Complex Organic Molecules in Perseus Embedded Protostars

To date, about two dozen low-mass embedded protostars exhibit rich spectra with complex organic molecule (COM) lines. These protostars seem to possess different enrichment in COMs. However, the statistics of COM abundance in low-mass protostars are limited by the scarcity of observations. This study introduces the Perseus ALMA Chemistry Survey (PEACHES), which aims at unbiasedly characterizing the chemistry of COMs toward the embedded (Class 0/I) protostars in the Perseus molecular cloud. Among the 50 embedded protostars surveyed, 58% of them have emission from COMs. A 56%, 32%, and 40% of protostars have CH$_3$OH, CH$_3$OCHO, and N-bearing COMs, respectively. The detectability of COMs depends on neither the averaged continuum brightness temperature, a proxy of the H$_2$ column density, nor the bolometric luminosity and the bolometric temperature. For the protostars with detected COMs, CH$_3$OH has a tight correlation with CH$_3$CN, spanning more than two orders of magnitude in column densities normalized by the continuum brightness temperature, suggesting a chemical relation between CH$_3$OH and CH$_3$CN and large chemical diversity among the PEACHES samples at the same time. A similar trend with more scatter is also found between all identified COMs, hinting at common chemistry for the sources with COMs. The correlation between COMs is insensitive to the protostellar properties, such as the bolometric luminosity and the bolometric temperature. The abundance of larger COMs (CH$_3$OCHO and CH$_3$OCH$_3$) relative to that of smaller COMs (CH$_3$OH and CH$_3$CN) increases with the inferred gas column density, hinting at an efficient production of complex species in denser envelopes.


INTRODUCTION
Planet formation may start during the embedded phase of star formation ). In the scenario where planets form from the embedded disks, resulting in substructures, the chemistry of embedded disks may play a significant role for the chem-ical composition of the forming planets. In the past two decades, observations show the emission of complex molecules toward the center of several protostellar cores. From the astronomical point of view, complex molecules are usually defined as species that contain 6 or more atoms (Herbst & van Dishoeck 2009). All detected complex molecules are organic. It can be saturated (e.g., Cazaux et al. 2003;Bottinelli et al. 2007;Jørgensen et al. 2016;Lee et al. 2017;Ceccarelli et al. 2017, and the ALMA PILS Survey) or unsaturated (e.g., carbon-chain molecules: Sakai & Yamamoto 2013;Sakai et al. 2014;Law et al. 2018). The saturated organic molecules, often called complex organic molecules (COMs) or interstellar COMs (iCOMs), have single covalent bonds between carbon atoms, making them rich in hydrogen, while the unsaturated organic molecules contain double or triple bonds between carbon atoms, making them poor in hydrogen. Most of the COMs appear at the inner warm envelope (e.g., Bottinelli et al. 2004a) and/or the surface of the embedded disk where the temperature is warm ( 100 K; e.g., Lee et al. 2017), while COMs in some sources, such as L1157, are linked to molecular shocks (e.g., Bachiller et al. 2001;Lefloch et al. 2017;Codella et al. 2020). The connection between a disk-like structure and the emission of COMs has been elusive (Belloche et al. 2020). However, if the disk formation can inherit COMs from the embedded envelope (e.g., Drozdovskaya et al. 2014), the chemistry of COMs at the embedded phase of star formation may implicate future developments of organics on the planets.
For low-mass protostars, gas-phase COMs typically emit from the warm inner protostellar envelope (T 100 K), which corresponds to 100 au for typical lowmass protostars (Bottinelli et al. 2004a;Ceccarelli et al. 2007). Thus, COMs serve as a tracer of the inner envelope where a disk may be forming (Aikawa 2013), measuring the chemistry and dynamics of embedded disks Lee et al. 2017). Other processes, such as jets and outflows, can also produce the emission of COMs by sputtering and/or shock chemistry (e.g., Arce et al. 2008;Lefloch et al. 2017). Furthermore, accretion shocks occurred when the gas falls onto the disk may enhance the desorption of COMs Miura et al. 2017;Csengeri et al. 2018). However, Belloche et al. (2020) showed no clear correlation between the COM emission and the occurrence of disk-like structures possibly due to specific conditions required for efficient desorption via accretion shocks (Miura et al. 2017) or insufficient sensitivity. Outbursting protostars, such as FU Orionis objects, can temporarily increase the COM abundance with elevated temperature (van 't Hoff et al. 2018b;Lee et al. 2019). In the more evolved stage, COMs have only been detected in a few disks (Öberg et al. 2015;Walsh et al. 2016;Favre et al. 2018;Podio et al. 2020).
While recent observations show abundant COMs in several embedded protostars, the probability of the existence of COMs and their relationship to the star formation process remain yet to be understood. Several protostars are rich in long carbon-chain molecules but have little COMs, such as L1527 (Sakai et al. 2010) and IRAS 15398−3359 (Sakai et al. 2009;Oya et al. 2014;Okoda et al. 2018). In comparison, interferometric observa-tions detected both COMs and carbon-chain molecules in several protostars, such as IRAS 16293−2422 (Jaber Al-Edhari et al. 2017;Murillo et al. 2018), L483 (Oya et al. 2017), and B335 (Imai et al. 2016). Single-dish surveys showed no correlation between COMs and long carbon-chain molecules (Graninger et al. 2016;Higuchi et al. 2018). Interestingly, carbon-chain molecules tend to distribute at a larger spatial scale than that of COMs (Imai et al. 2016;Oya et al. 2017;Bouvier et al. 2020). Aikawa et al. (2020) demonstrated that the conditions of prestellar cores can result in larger spatial extent of carbon-chain molecules coexisting with COMs or deficient of carbon-chain molecules in COM-rich sources. Moreover, dust opacity may obscure the emission of COMs in submillimeter wavelengths, introducing additional uncertainty to our understanding of COMs (De Simone et al. 2020). Therefore, the chemical pathways of complex molecules at the embedded protostars remain ill-constrained, highlighting the need of unbiased chemistry surveys to characterize the statistics of COMs (e.g., Caux et al. 2011;Lefloch et al. 2018).
The Perseus molecular cloud is one of the most active nearby star-forming regions, which extends ∼10 pc on the sky. Many studies assume a distance of 235 pc for Perseus based on the maser observations toward NGC 1333 SVS 13 (Hirota et al. 2008). Recently, Ortiz-León et al. (2018) estimated a distance of 321±10 pc for IC 348 using VLBA and a distance of 293±22 pc for NGC 1333 using Gaia parallaxes. Zucker et al. (2020) combined the Gaia parallaxes and photometric data with a Bayesian framework to revise the distances toward several sightlines of Perseus, resulting in distances ranging from 234 pc to 331 pc.
In Perseus molecular cloud, infrared and submillimeter surveys reveal more than 400 young stellar objects as well as ∼100 dense cores, which contains ∼50 Class 0 and I protostars (Hatchell et al. 2005;Jørgensen et al. 2008;Dunham et al. 2013; hereafter the embedded protostars). The Perseus molecular cloud contains starforming regions in a wide range of environments. The majority of protostars in Perseus are associated with the two clusters, NGC 1333 and IC 348. NGC 1333, which has many active outflows (Knee & Sandell 2000;Plunkett et al. 2013), was thought to be younger than IC 348 based on the higher extinction and abundance of protostars; however, Luhman et al. (2016) showed a similar age for two regions based on the evolutionary modeling of young stars and brown dwarfs. The embedded protostars associated with IC 348 lie at the southwest of the open cluster, IC 348, near a prominent outflow, HH 211 (e.g., Lee et al. 2009). Most of other embedded protostars at Perseus are related to L1448, L1455, Barnard 1 (B1), and Barnard 5 (B5). L1448 has active outflows that may regulate the ongoing star formation (Curtis et al. 2010); The low abundance of Class 0 protostars in L1455 suggests a more evolved protostellar population than that in other regions (Hatchell et al. 2007); B1 exhibits rich spectra of deuterated species Marcelino et al. 2005). Furthermore, the protostellar multiplicity in Perseus has been extensively studied by the VLA Nascent Disk and Multiplicity Survey (VANDAM) of Perseus protostars, showing a multiplicity fraction of 0.57±0.09 and 0.23±0.08 for Class 0 and I protostars, respectively Tychoniec et al. 2018;Segura-Cox et al. 2018;Tobin et al. 2018). Murillo et al. (2016) further characterized the spectral energy distributions of all sources identified in the VANDAM survey to derive protostellar properties, such as bolometric luminosities (L bol ), bolometric temperatures (T bol ), and evolutionary stage classification. Thus, the Perseus molecular cloud provides an ideal test bed for chemistry in embedded protostars by surveying the protostars within each region and across the entire cloud. Higuchi et al. (2018) presented a pilot survey of the chemistry in the Perseus embedded protostars with Nobeyama 45m telescope and IRAM 30m telescope, which surveyed all Class 0/I protostars (Hatchell et al. 2007) that have L bol greater than 1 L (0.7 L for protostars in B1 and B5) and the envelope mass greater than 1 M . This pilot survey probes the molecules such as C 2 H, c-C 3 H 2 , and CH 3 OH. The majority of the sources have emission of both CH 3 OH and carbonchain molecules. They suggested a possible correlation between the location of sources within the clouds and the ratio of CH 3 OH/C 2 H indicating the environmental effect on the chemistry. However, single-dish observations of CH 3 OH are contaminated by the contribution of the photodissociation region of the molecular clouds hosting the protostar (Bouvier et al. 2020). In addition, chemical surveys with single-dish telescopes are often less sensitive to detect the COMs in the vicinity of the protostars. To trace those COMs , high spatial resolution from interferometric observations is crucial. Based on the NOEMA observations of 26 embedded protostars in several molecular clouds, Belloche et al. (2020) detected CH 3 OH and other COMs in about a half of them. Interestingly, CH 3 OH and CH 3 CN show strong correlation, while their chemical link remains unclear. They also investigated the origins of COMs and revealed apparent chemical difference among binary systems in their survey. Internal luminosity is the most impacting parameter for the COM chemical composition, while the existence of a disk-like structure has no obvious impact to the COM emission. A combined analysis of a few protostars from four different star-forming regions (Perseus, Serpens, Ophiuchus, and Orion) showed no significant correlation between the occurrence of COMs and the L bol (van Gelder et al. 2020). Furthermore, the abundance of a few O-bearing COMs to CH 3 OH are similar for the sources in different regions. They suggested an inheritance scenario where the evolution in the prestellar phase dominates the chemistry. Despite the growing sample of COM detections in embedded protostars, the statistics of COM abundance remain unconstrained due to the biases from the source selection, limited resolution, and small sample sizes. To unbiasedly understand the chemistry, we conducted the Perseus ALMA Chemistry Survey (PEACHES) that probes the complex chemistry toward Perseus embedded protostars. The source selection followed the same criteria as Higuchi et al. (2018) with a few pointing modifications using the results from the VANDAM of Perseus protostars . The PEACHES survey also probed the outflows of embedded protostars (e.g., NGC1333 IRAS 4C: Zhang et al. 2018).
Section 2 describes the details of our ALMA observations. Section 3 presents the identification of protostellar sources and the methodology of spectral extraction, line identification and modeling. Section 4 shows the detection statistics of COMs and their correlations. Section 5 discusses the implications of the observed COM abundance on the current understanding of chemistry of COMs at protostellar sources. Finally, Section 6 summarizes the findings of this study.

OBSERVATIONS
The PEACHES observations were conducted in two ALMA projects (2016.1.01501.S and 2017.1.01462.S; PI: N. Sakai), surveying 37 fields toward the Perseus molecular cloud. Each project surveyed different regions in the Perseus molecular cloud. Table 1 lists the basic information for each project, while Table 2 shows basic information of the sources observed in the PEACHES survey. The ALMA correlator was configured to have 13 spectral windows at Band 6, which have 12 narrow spectral windows with 480 channels and a wide spw with 980 channels. The narrow spectral windows were tuned to observe specific molecular species, such as SiO, C 2 H, CS, CH 3 OH and CH 3 CN, with a spectral resolution of 122 kHz (∼0.15 km s −1 ), while the wide spectral window aimed to observe the continuum with a resolution of 0.976 MHz (∼1.2 km s −1 ). In later imaging processes, we combined the two spectral windows potentially affected by the broad SiO emission, resulting in a total of 12 spectral windows. The frequency setups for the two ALMA projects are largely identical except for the wide continuum window, which shifts by ∼380 MHz. Table 3 lists the frequency ranges for each spectral window and the corresponding channel widths. We used CASA (McMullin et al. 2007) for standard calibration and imaging of the continuum and spectral lines with tclean. Self calibration is not applied in this study, since some of the samples have weak continuum emission. Because of the rich spectra, we manually flagged the lines for the continuum imaging. The images were cleaned down to 0.022Jy beam −1 , except for the continuum spectral windows, whose images were cleaned down to 0.008Jy beam −1 due to their lower spectral resolution. The line imaging used the "multiscale" deconvolver with a robust parameter of 0.5, because the targeted emission traces different spatial scales (e.g., SiO for outflows and COMs for the inner envelope/disk). Lastly, we applied the primary beam correction to the image cubes. The synthesized beam is about 0. 6×0. 4 averaged across all spectral windows. In this study, we assume a distance of 300 pc for the entire Perseus cloud, resulting in a synthesized beam of 180 au×120 au.

Identification of Young Stellar Objects
In the 37 fields of view, 51 continuum peaks are identified. We used the CASA task imfit to iteratively fit a 2D Gaussian profile for continuum emission down to 5σ of the residual image within the central 70% of the size of the primary beam (20 ). We measured the noise of the continuum emission (σ) from the vicinity of the continuum emission. For the field centered on Per-emb 16, which also covers Per-emb 28, the fitting used a threshold of 4σ and extended the mask to the entire primary beam, because Per-emb 28 is detected toward the edge of the primary beam, where the noise is higher. The continuum of the multiple systems were manually fitted to the individual continuum peaks. The two protostars, L1448 IRS 2E and SVS 3, become non-detections in our observations due to their low brightness. SVS 13C locates at the edge of the primary beam, resulting in a noisy continuum. Thus, we excluded SVS 13C from this study, making a total of 50 protostars. Table 4 lists the properties ot the fitted continuum emission, while Figure 1 shows the continuum emission along with the fitted shapes. The continuum emission peak appears as a compact circular or elliptical shape. Some sources show extended continuum emission resembling the shape of outflow cavities, such as Per-emb 22 A and B. The source velocities collected from the literatures agree with our observations. Three sources, EDJ2009-237, Per-emb 60, and EDJ2009-172, show no spectral line and no reliable measurement of source velocity in literature; therefore, we excluded them from spectral extraction as well as the line identification and modeling. However, these three sources still contribute to the total number of sources for characterizing detection statistics.
The observations resolved or marginally resolved 94% (47 of 50) of the continuum sources. Our sample includes single sources as well as resolved and unresolved multiple systems. According to the VANDAM survey , nine close binary systems remain unresolved in the PEACHES survey (Table 4).
The mean continuum brightness temperature (T b,cont ) traces the optically thin dust mass; however, the dust emission from the center of protostars may be marginally optically thick at ALMA Band 6 (Ko et al. 2020), making the derived dust mass a lower limit. The T b,cont of PEACHES sample correlates with the dust mass of the disk-like structure of Perseus protostars derived from the 1 mm ALMA and 9 mm VLA observations ) except for that of Per-emb 54 ( Figure 2). The 1 mm masses are systematically lower than the 9 mm masses because the inner disk-like structure becomes optically thick at 1 mm. Nonetheless, the T b,cont increases linearly with the 1 mm masses with a similar slope as that with the 9 mm masses, suggesting that the T b,cont traces the central mass without significant impact of optical depth. The 1 mm mass of Per-emb 54 follows the linear trend, indicating that the deviation of its 9 mm mass may be due to excessive free-free emission.
Seven sources in our sample have T b,cont > 10 K (Table 4). Per-emb 12 A has a T b,cont of 21.9 K, the brightest among our sample. In fact, continuum opacity limits the detectability of COMs toward 4A1, where COMs show  . Toward the other six sources with T b,cont > 10 K, Per-emb 33 A, Per-emb 12 B, Per-emb 13, Per-emb 5, B1-b S, and Per-emb 11 A (see Section 3.2 and figures mentioned there), our observations detected several molecular lines from the central region, suggesting that the continuum opacity has less impact compared to that for Per-emb 12 A. Previous observations also suggest less impact from the continuum opacity for Per-emb 13 (Belloche et al. 2020) and B1-b S (Marcelino et al. 2018). Most of the sources have the mean brightness temperature lower than 10 K, while spatially resolved embedded disks tend to have a warm temperature, > 20 K within a radius of 100-150 au (van 't Hoff et al. 2018a(van 't Hoff et al. , 2020a. Thus the continuum source is likely to be unresolved and have a insignificant impact on the COMs emission at the angular resolution of this study. Note that some of the sources such as Per-emb 18 and Per-emb 21 show systematical shift of the peaks in multiple species. This could be understood by the high dust opacity at the continuum peak. However, correlations between the column density of COMs (Section 4.2.2) would be less affected since the dust opac-ity effect could be similar for all the species. We take the T b,cont as a tracer of the averaged gas column densities, judging from a reasonable correction in Figure 2. Constraining the actual column densities requires higher resolution observations at multiple frequency bands and insights on the dust properties.

Spectral Extraction
The ALMA image cubes were post-processed to extract 1D spectra for identifying the emission of complex molecules and further analyses. Given the compact size of COM emission ( 100 au) and the spatial resolution of ∼0. 5 (∼150 au at 300 pc), we focused on the spectra toward the continuum sources to search for the COMs in the inner envelope. Figure 3 shows representative sample of the COMs emission, while the maps for the rest of the sample are shown in Appendix A. In most of the cases, COMs emission concentrates around the protostars at 300 au scale. Three steps of post-processing reduced the image cubes to 1D spectra, which are summarized below.
• Extracting spectra: We used the CASA task specflux to extract the mean flux density within the ellipse which has the same major and minor axes as well as the position angle as the fitted continuum sources (see also the dashed ellipses in Figure 1).
• Baseline subtraction: The continuum has been removed before the imaging process; however, the a The "binary" label indicates that the continuum source is an unsolved binary in our observations, according to Tobin et al. (2016). b The deconvolved continuum source is point source no larger than the listed size. c The continuum cannot be deconvolved because the emission may be resolved in only one direction. extracted spectra sometimes still show imperfect baselines due to rich spectra of COMs and broad features. Thus, we manually selected the frequency ranges for baseline subtraction for each spectral window and each field.
• Velocity correction: Finally, the frequencies of the extracted spectra were corrected according to the source velocities. We collected the source velocities from the literature as well as from the strong emission lines in our spectra, such as SO and CS. Table 2 lists the adopted source velocities and the corresponding references.

Line Identification
Line identification starts with manual identification and verification for a few sources with rich spectra, including Per-emb 12 B and B1-b S. Then, we tested the list of identified species to the rest of sources for identification. We used splatalogue 1 to identity the molecular species, where the molecular data come from the Cologne Database of Molecular Spectroscopy (CDMS; Müller et al. 2001Müller et al. , 2005Endres et al. 2016) and the Jet Propulsion Laboratory Millimeter and Submillimeter Spectral Line Catalog (JPL; Pickett et al. 1998). Appendix B and Table B1 summarize the catalogs relevant for each species. Any tentatively identified species requires verification using synthetic spectra modeled with the xclass (Möller et al. 2017), which performs LTE radiative transfer calculations using the molecular data from CDMS and JPL. We describe details about the xclass modeling in Section 3.4. An identification needs to satisfy the following criteria.
• The predicted strengths in synthetic spectra agree with the observations, assuming a fiducial column density that produces the observed intensity, an excitation temperature of 100 K, and a source size of 0. 5. Identification of a species requires an unblended line detected at S/N > 3.
• For each species, at least one spectral line is not blended with other emission. • Identified molecules need to be already found toward young stellar objects as summarized in McGuire (2018).
The identified species and transitions that are detected in at least one of the PEACHES sample are listed in Appendix C. Only identifiable transitions are listed. The xclass modeling includes all the transitions for the identified species in our frequency coverage so that we can test the relative strengths of different transitions in the identification process. Not all emission is identified; we reserve the analysis of the unidentified lines in future studies.

Modeling the Spectra of COMs
In addition to verifying the identification, we use the xclass to model the observed spectra to constrain the column densities of COMs. The xclass solves radiative transfer equation for an isothermal source in one dimension, called detection equation (Stahler & Palla 2004). The source brightness distribution follows a 2D Gaussian distribution, described by its full width half maximum (F W HM COM ), the column density within F W HM COM (N COM ), and the excitation temperature (T ex ). The model has two additional parameters, the line width (∆ν) and velocity offset (v off ), to calculate the spectra. For each species, the xclass solves the radiative transfer equations independently. For a few species whose emission blends at some frequencies, we fit them together in pairs, such as CH 3 OH & CH 3 OCHO, HDCO & 13 CH 3 OH, CH 3 CHO & CH 2 DOH, 34 SO & C 2 H 5 OH, and CH 3 OCH 3 & CH 2 DCN. We detect several simple molecules, such as CS, H 13 CN, HDCO, SO, and SO 2 . However, these lines often exhibit double-peaked or complex line profiles (see Figure 5 for an example) due to the self-absorption or contamination from the envelope and outflows, which will be analyzed in future studies. While our modeling includes these species to account for the weaker emission in the spectral coverage, we excluded these species from the following analysis because of the limitation of the simple LTE model.
To test the assumption of LTE, we ran 1D non-LTE models to estimate the discrepancy between the kinetic temperature (T k ) and the excitation temperature (T ex ) for the emission of CH 3 OH, which often is the most abundant COMs (e.g., Belloche et al. 2020). We followed a similar procedure demonstrated by Jørgensen et al. (2016) using RADEX, an 1D non-LTE radiative transfer package (van der Tak et al. 2007). The collision rates were taken from the Leiden Atomic and Molecular Database (Schöier et al. 2005) based on calculations by (Rabli & Flower 2010  observations, the T k deviates from the T ex by more than 10% at n(H 2 ) < 10 10 cm −3 , or the CH 3 OH column density 10 16 cm −2 ( Figure 4). The disagreement increases toward lower gas density where the negative T ex , indicative of a masing condition, leads to the spikes in Figure 4. Assuming optically thin continuum emission, we estimated a gas density ranging from ∼ 4 × 10 9 − 10 12 cm −3 for the continuum sources in our sample (Appendix D). In this range, T k agrees with T ex within 10-30%. Our calculations for the CH 3 OH line at 243915.79 MHz show a maximum difference smaller than 0.1% for the combinations of parameters. The 243915.79 MHz CH 3 OH line has ∆K = 0 so that it is not a maser line, which requires ∆K = 0. Therefore, LTE is a valid assumption for modeling the spectra of COMs.  . The difference between the excitation temperature (Tex) and the kinetic temperature (T k ) normalized to T k for the CH3OH 211 → 101 transition at 261805.68 MHz as a function of the number density of H2 at different T k . The calculations were performed using RADEX. The results with CH3OH column densities of 10 16 , 10 17 , and 10 18 cm −2 are shown in dotted, dashed, and solid lines, respectively, while the results with T k of 100, 200, and 300 K are shown in decreasing transparency. The legends indicate the combination of NCH3OH and T k as (NCH3OH, T k ).
Our observations cover one or a few transitions for each species, making the determination of excitation temperature uncertain. Thus, instead of choosing an excitation temperature, we optimized the model with five different T ex , 100, 150, 200, 250, and 300 K. Then, we took the mean and the range of fitted column density as the best fitted value and the corresponding uncertainty, respectively. For the COMs with A-and E-species due to the nuclear spins, such as CH 3 OH and CH 3 CN, we assumed that the A-and E-species have the same excitation temperature so that they can be modeled together. These T ex covers the typical range of T ex for the COMs detected toward embedded protostars (e.g., Bergner et al. 2019). We derived a similar range of T ex for the species that have multiple transitions detected, such as CH 3 OH and CH 3 OCHO (Appendix E), suggesting that the range of 100-300 K represents the typical T ex of COMs in the PEACHES sample. Therefore, we can take the range of the derived column densities as an estimate of the true column density without fine tuning the model for each source. To have a robust fitting result, we fixed F W HM COM as 0. 5 and assumed v off = 0. We also assumed a negligible dust optical depth, making the fitted column density a lower limit if the dust opacity turns out to be substantial. The synthetic spectra were calculated for the size of continuum emission (Table 2) to include the effect of beam dilution. We allowed the line width varying between 1.2 km s −1 and 3.5 km s −1 for better fitting quality and chose a range of the column density for each molecule according to the strength of the emission.
The fitting function in xclass includes several optimization algorithms that can be used in sequence to reduce biases. We configured the algorithm chain that starts with the genetic algorithm followed by the Levenberg-Marquardt χ 2 minimization. The genetic algorithm iteratively searches for the best-fitting parameters with generations that evolve like a natural selection, where the better fitting models get less modification over generations. We set up the genetic algorithm to search for the top two best-fitting models with 30 generations. Then, the Levenberg-Marquardt χ 2 minimization applies to the two best-fitting models for 20 iterations to find the best-fitting models. The genetic algorithm aims to find global minimums and the Levenberg-Marquardt χ 2 minimization further finds the best-fitting models in the global minimums. The two best-fitting models found by the genetic algorithm often very similar, suggesting that there is only one minimum. In the rare cases of two distinct global minimums, we picked the model with the lower χ 2 values from the two best-fitting models constrained by the Levenberg-Marquardt χ 2 minimization. We allowed the column density ranging by a few orders of magnitude to ensure that a global minimum can be found.
From the five fitted column densities assuming five different T ex , if a molecule is detected according to the criteria listed in Section 3.3, we took the mean column densities as the best-fitting column density, while the range of the column densities indicates the uncertainty. If a species is considered non-detection (Section 3.3), we used the synthetic spectra to derive the upper limit of the column density corresponding to the noise of the spectra. The fitted column density assumes a FWHM of 0. 5 (150 au). Figure 5 shows an example of the fitted spectra for Per-emb 44 (SVS 13 A), assuming T ex = 150 K. The fitted spectra may underestimate or overestimate the strengths of some transitions because the assumed excitation temperature may not be the true excitation temperature.
To benchmark this hybrid optimization process, we compared the fitted column densities with the fitting using the Monte Carlo Markov Chain (MCMC) algorithm on the most chemically rich source, Per-emb 44 (SVS 13 A). The MCMC algorithm uses the affine-invariant MCMC package emcee (Foreman-Mackey et al. 2013) to sample the parameter space. In this MCMC optimization, we fitted both the excitation temperature and the column density. For a single-species fitting, we used 100 walkers with 100 iterations after 30 burnin iterations; for a combined fitting of two species, we used 500 walkers to better sample the parameter space. The fitting column densities from the hybrid (genetic and Levenberg-Marquardt χ 2 minimizations) method are consistent with that fitted with the MCMC optimization ( Figure 6). Also, our hybrid method yields comparable uncertainties as the ones from the MCMC method, validating our modeling approach.

Detection Statistics
We summarize the detection statistics in Figure 7, which includes COMs, carbon-chain molecules, and simple molecules, such as CS, H 13 CN, SO, 34 SO, and SO 2 . To derive the detection fraction, we included the three sources that are excluded from modeling due to no reliable source velocity (Section 3.1), making a total of 50 sources. The COMs discussed here are derived from the spectra taken from the area surrounding the continuum peak as explained in Section 3.1. We focus on the chemistry of COMs toward the disk-forming regions in this study; therefore, this study excludes any molecules that are only detected outside the region of continuum emission. A comparison of chemical composition in protostellar envelopes (100-1000 au) requires observations with larger maximum recoverable scale (θ MRS ).
The PEACHES sample show a diverse chemistry of COMs from no molecule detected (e.g., B1-b N and L1455 IRS 2) to rich spectra of COMs (e.g., Per-emb 13). Most of the PEACHES sample have simple molecules, such as SO, CS, H 13 CN, and HDCO, and ∼ 60% of sources have SO 2 and 34 SO. Emission of C 2 H can be easily identified from the spectra. However, the C 2 H toward the continuum sources often shows irregular line profiles together with velocity offsets and absorption (Appendix F, Figure 22). Warm environments, such as the outflow cavity wall, easily enhance the abundance of C 2 H because of the elevated abundance of C + (e.g., Zhang et al. 2018;Imai et al. 2019). Thus, C 2 H emission tends to extend along with the outflow cavities, making the 1D spectra toward the continuum source unrepresentative of the abundance of C 2 H. Future studies will investigate the comparison between unsaturated species and COMs using follow-up observations with a higher dynamic range of the spatial distribution (from ∼100 au to several 1000 au scale). The previous observations by single-dish telescopes has beam size of ∼3000 au (e.g., Higuchi et al. 2018). Thus, observations that fill up the spatial scale from several 100 au to 1000 au are necessary.
Several sources have broad SiO emission, significantly contaminating the emission of CH 3 CH 2 CN and CH 3 CHO. In comparison, emission of COMs has relatively narrow line width of a few km s −1 . In the later quantitative discussion, we exclude the spectral windows contaminated by the SiO emission. For a few sources, such as Per-emb 26, we can still distinguish the emission of CH 3 CH 2 CN and CH 3 CHO from the broad SiO emission and include these identification in the detection statistics. Figure 7 shows the number of COMs detected toward the PEACHES sample. Twenty-nine (58%) sources have COMs detected. CH 3 OH is detected in 28 sources (56%); CH 3 OCHO is detected in 14 sources (28%); and N-bearing COMs are detected in 20 sources (40%). Comparing to the COMs in the CALYPSO survey (Belloche et al. 2020), the fraction of sources that have methanol, ∼50%, is similar to that for the PEACHES sample. Also, 30% of the CALYPSO sources have at least three species of COMs, while 28% of the PEACHES sample have at least three species of COMs. In a smaller sample selected from several molecular clouds, van Gelder et al. (2020) found 3 of 7 protostars show emission of COMs.
To investigate the impact of protostellar properties on the appearance of COMs, we divided the PEACHES sample into four groups according to the number of detected COMs from no COM, at least one, five, and ten species of COMs (Figure 8). The compared protostellar properties include the bolometric luminosity (L bol ), the bolometric temperature (T bol ), and the mean continuum brightness temperature (T b,cont ), where the T bol and L bol are collected from Tobin et al. (2016) and Murillo et al. (2016). The T bol traces the evolutionary stage of protostars with larger T bol for more evolved protostars (Myers & Ladd 1993;Chen et al. 1995;Evans et al. 2009); and the L bol serves as a proxy of the size of the warm inner envelope as the accretion luminosity t-HCOOH dominates the L bol for Class 0/I protostars (Dunham et al. 2014). The number of detected COMs shows no obvious trend with L bol , T bol , and T b,cont . If the COMs mostly come from thermal desorption, the region with T > T desorption may be smaller for low luminosity sources, resulting in fainter emission of COMs and reducing their detectability. In contrast to this expectation, the COMs detection shows no obvious trend with L bol , except that the group with more COMs has a higher minimum L bol . However, the range of L bol in each group is much larger than the change in the minimum L bol ; thus, we cannot establish a significant correlation between the appearance of COMs and L bol .
The high T bol value for the sources without COM detection is driven by the high T bol of EDJ2009-172 (T bol =1100 K), while the second highest T bol for the group is L1455 IRS 2 with 740 K. The lower maximum T bol for the sources with more COMs detected hints that evolved sources may have less COMs in gas-phase. In Ophiuchus molecular cloud, Artur de la Villarmois et al.
(2019) also find no detection of CH 3 OH toward a sample of Class I protostars. If we take T bol of 70 K as the boundary between Class 0 and I, our survey shows no obvious difference for the Class 0 and I sources; however, the evolutionary stage classified with T bol may not represent the true evolutionary state of the protostar due to extinction and the effect of inclination (Launhardt et al. 2013;Dunham et al. 2014). The detectability of COMs may depend on other factors, such as the chemistry in the envelope scale, instead of the current physical condition. Thus, the chemo-physical history of protostars may play a dominant role in the chemistry of COMs (e.g., Garrod et al. 2008;Aikawa 2013;Aikawa et al. 2020).
Several COM-rich protostars in Perseus are previously known (e.g., Bottinelli et al. 2004bBottinelli et al. , 2007Sakai et al. 2006;Marcelino et al. 2018;Bianchi et al. 2019). Our survey detects more COMs in additional to the known CH 3 OH toward Per-emb 35 A and Per-emb 11 A (Figure 8). Per-emb 35 A is in a binary system (NGC 1333 IRAS 1) along with Per-emb 35 B. Source A has emission of CH 3 OH, CH 3 CN, CH 3 OCHO, and t-HCOOH, while Source B only has weak emission of CH 3 OH, showing an apparent chemical differentiation. Per-emb 11 A is the brightest source in a wide triple system, where Per-emb 11 B is ∼3 away toward southwest and Per-emb 11 C is ∼9. 5 away toward northeast ). Per-emb 11 A has rich spectra of COMs, including CH 3 OCHO, CH 3 OCH 3 , CH 3 COCH 3 , CH 3 CHO, and CH 3 CN.
Compact emission of COMs is consistent with an origin of thermal desorption at the inner envelope, so called hot corinos (Ceccarelli 2004). However, we cannot rule out other origins, such as the local enhancement due to accretion shocks , because our observations only marginally resolved most emission of COMs. Only L1448 IRS 3A shows extended CH 3 OH emission, and B1-b N shows no CH 3 OH emission toward the continuum but peaks at ∼1 away from the continuum, hinting the existence of CH 3 OH toward the continuum source blocked by opaque continuum (Marcelino et al. 2018). The enrichment of COMs due to outflows via sputtering (Arce et al. 2008) or photodesorption in the cavity walls (Drozdovskaya et al. 2015) dominates at a larger scale (>1000 au) and can be time-dependent, therefore playing a minor role for the COMs detected toward the continuum sources. Excluding L1448 IRS 3A, 28 sources (56%) are likely to have a hot corinos type chemistry. However, studies on the origin of COMs require higher spatial resolutions and are beyond the scope of this study.

Column Densities
From the method described in Section 3.4, we constrained the column densities of COMs along with their uncertainties (Table 5). Our modeling pipeline successfully reproduced the spectra of most PEACHES sample, except for Per-emb 17 due to its double-peaked line profile, which we further discuss in Section 4.2.1. Using the fitted column densities, we discuss the correlation between COMs in Section 4.2.2.

The Double-peaked Features of Per-emb 17
Per-emb 17 is a close binary system with a separation of 0. 278±0. 014 (83.3±4.0 au with a distance of  300 pc) ). Our observations show clear double-peaked features for most emission (Figure 9). The position-velocity diagram of the CH 3 OH have a bar-shape morphology with asymmetric bright spots, which may come from an unresolved disk (e.g., Lee et al. 2017;Yang et al. 2020). However, the nature of this double peaked features remains unclear due to the insufficient spatial resolution in our observations. The double-peaked feature challenges our standard modeling pipeline that assumes a simple Gaussian line profile. Thus, we modelled the spectra of Per-emb 17 separately from the entire sample. In xclass, we set up two components of molecular gas using the same assumptions described in Section 3.4. Then, we configured the fitting to allow both components varying their velocity offset from the rest frequency by ±5 km s −1 . The two velocity com-ponents for the emission of COMs appear at 1.96 km s −1 and −3.15 km s −1 on average, while the simple organics appear at slightly smaller velocity offsets. We reserved detailed analyses about the kinematics of Per-emb 17 for future studies. The fitting results become less robust with the additional parameters of velocity offsets due to the blending of the double-peaked lines. Thus, we performed a two-components fitting with fixed velocity offsets at 1.96 km s −1 and −3.15 km s −1 , which reproduces the observations (Figure 9). To minimize the uncertainty due to substantial line-blending, we estimated the total column density from the sum of the two components instead of reporting them individually. Then, we applied the same procedure to derive the mean column densities and upper limits as discussed in Section 3.4. Line-blending still prevents accurate determinations of  column densities for a few species. In particular, the emission of CH 3 OCH 3 at ∼259310 MHz may be underestimated or severely contaminated by the emission of CH 2 DCN, which has another transition at 260523 MHz that may be contaminated by SiO emission.

Correlations of COMs
The chemical evolution of protostars may leave certain patterns in the abundance of molecules as the dynamical evolution determines the density and temperature structures, which regulate chemical reactions. Thus, the abundance of COMs and their correlations provide critical information to study the chemical evolution at embedded protostars. As described in Section 3.4, we fitted the column density and line width with different excitation temperatures, resulting in a range of column density as its uncertainty. Thus, we can investigate the correlation between the column density of each species to probe the underlying chemistry.
To quantify the goodness of correlation, we calculated the Pearson's correlation coefficient (r), which tests the linearity of two variables. We used the bootstrap method to sample the fitted column densities to calculate the Pearson's r. The bootstrap method is an iterative process that resamples the column densities and calculates the Pearson's r for the drawn sample for each iteration. We took the correlation between CH 3 OH and CH 3 CN as an example to demonstrate this bootstrap method (Figure 11). For detections, we assumed an asymmetric normal distribution centers on the best-fitted column density to account for the asymmetric uncertainty in logarithmic scale. The distribution on either side of the best-fitted value has an equal probability, where the width of each "half" normal distribution follows the corresponding uncertainty. We ran the bootstrap process for 1000 iterations to characterize the distribution of the Pearson's r. Including the non-detections for bootstrapping the correlation coefficient requires an assumption of their underlying distribution that is consistent with zero. However, assigning a distribution of the column densities of two species requires an assumption of their covariance, which determines the correlation coefficient. Thus, we only considered detection for the calculations of the correlation coefficient. Figure 12 shows the correlations of several COMs selected from their detection rates as well as the ratios between species, which are discussed in Section 5.2. The column density of CH 3 OH best correlates with that of CH 3 CN (see also Figure 11). Belloche et al. (2020) also found the tight correlation between these two molecules from the CA-LYPSO survey, which has a smaller sample, 26 proto-  stars, in several molecular clouds. The column densities of CH 3 OCH 3 and CH 3 OCHO also show a tight correlation. Jaber et al. (2014) also show a tight correlation between CH 3 OCH 3 and CH 3 OCHO from different ISM sources, ranging from galactic center clouds, hot cores, protostars, cold clouds to comets. Moreover, CH 3 OH correlates with CH 3 OCHO and CH 3 OCH 3 , which are the daughter species of CH 3 OH, with a weaker correlation than that with CH 3 CN. The correlations be-tween CH 3 CN and CH 3 OCHO or CH 3 OCH 3 also show the same behavior. The decreasing correlation between CH 3 OH or CH 3 CN and the daughter species of CH 3 OH, CH 3 OCHO and CH 3 OCH 3 , is consistent with the scenario where CH 3 OH and CH 3 CN turns into CH 3 OCHO and CH 3 OCH 3 along with other COMs as the chemistry evolves (e.g., Garrod & Herbst 2006  tions are driven by a few sources due to the low detection rate of other COMs (Figure 13). To directly compare the abundance of COMs, we normalized the column densities of COMs by the T b,cont , which is a proxy of gas column density. The normalized column densities of COMs show similar correlations as that of the column densities ( Figure 14, red markers). We further normalized the ratio of column densities to T b,cont with L bol and T bol to test whether the correlation seen in Figure 14 (red markers) is dominated by the pro-tostellar properties. After the normalization with L bol and T bol (Figure 14, blue and black markers), the correlation between CH 3 OH and CH 3 OCH 3 weakens significantly, while CH 3 CN also shows a considerable lower correlation with CH 3 OCH 3 . The correlations between CH 3 OCHO and CH 3 OCH 3 remain similar with different normalizations. Per-emb 17 is an outlier with low and uncertain abundance of CH 3 OCH 3 (Section 4.2.1), reducing the correlation strength. However, the correlation analyses excluding the abundance of Per-emb 17 show similar trends. In summary, CH 3 CN has the best correlation with CH 3 OH, followed by that of CH 3 OCHO and CH 3 OCH 3 . For the protostars with compact emission of COMs, so called hot corinos, the abundance of COMs correlates well between species. With the normalization of T b,cont , the Pearson's r for the correlations between the four most detected species (CH 3 OH, CH 3 CN, CH 3 OCHO, and CH 3 OCH 3 ) range from 0.76 to 0.93 ( Figure 12). The Pearson's r for the same four species with all COM species has a median value of 0.86 with a range from 0.57 to 0.95 (Figure 13). The correlations remain unchanged with the normalizations of L bol and T bol , suggesting that the correlation, which represents the chemistry of COMs, may be independent of the evolutionary stage of protostars as also suggested by Bianchi et al. (2019). The source size and beam dilution have limited impact to the derived column densities (see the discussion in Appendix G). The limited number of detections of the COMs except for the four most detected COMs hinders further quantification of their correlation strengths. While the PEACHES sample show similar abundance ratios between COMs, the absolute column densities of COMs vary by 1-3 orders of magnitude, suggesting a diverse environment in the PEACHES survey.
To understand chemical similarity of COMs seen in embedded protostars in different environments, we further test the effect of protostellar evolution with the abundance ratios of the four most abundant COMs ( Figure 15). Per-emb 17 has a significantly lower column density of CH 3 OCH 3 compared to that of the PEACHES sample, resulting in apparent outliers in Figure 15 (opened circles). The ratios of CH 3 CN/CH 3 OH and CH 3 OCH 3 /CH 3 OCHO, which are molecules with a similar complexity, have little variation as the functions of L bol , T bol , and T b,cont . In contrast, the ratios of CH 3 OCHO or CH 3 OCH 3 to CH 3 OH which are the molecules with more complexity over the molecules with less complexity, increases with T b,cont , whereas the ratios tentatively decrease with L bol and T bol . Because of its high abundance, CH 3 OH may be more optically thick at a higher gas column density (high T b,cont ), resulting in an underestimation of CH 3 OH column density, hence the elevated ratios. The emission of CH 3 CN is likely to be optically thin given its low abundance compared to that of CH 3 OH. In fact, the highest fitted column density of CH 3 CN, ∼ 10 16 cm −2 , and the lowest modeled excitation temperature, 100 K, results in an optical depth of ∼0.4. The CH 3 OH becomes optically thick when the column density exceeds 2 × 10 17 and 5 × 10 17 cm −2 at 100 K and 150 K, respectively. If the temperature is 200 K or higher, CH 3 OH remains optically thin when N≤ 10 18 cm −2 , while our modeling estimates the highest column density of 1.1 × 10 18 cm −2 in Per-emb 27. If the optical depth of CH 3 OH plays a major role in the trend between the CH 3 CN/CH 3 OH and CH 3 OCH 3 /CH 3 OCHO to T b,cont , we would ex-  Figure 12. Corner plot of the correlations of the column densities between CH3OH, CH3CN, CH3OCHO, and CH3OCH3. The color code follows that in Figure 11. The annotated texts indicate the Pearson's r (r d ) and the logarithmic ratio of the two molecules (Ny/Nx) for the detection-only sample. The four most detected COMs are shown in this figure, while other COMs are shown in Figure 13.
pect to see a similar trend in the CH 3 CN/CH 3 OH to T b,cont , which shows a flat relation instead. Therefore, the optical depth effect is unlikely to dominate the observed trend between the CH 3 CN/CH 3 OH and CH 3 OCH 3 /CH 3 OCHO to T b,cont . Chemical evolution is an alternative scenario for such trend. Under the scheme of grain-surface chemistry, CH 3 OH primarily forms from the hydrogenation of CH 3 O or CH 2 OH on grains, while CH 3 OCHO and CH 3 OCH 3 primarily form from the radicals associated with H 2 CO and CH 3 OH destruction, which requires T= 20 − 40 K so that the radicals become mobile (Garrod & Herbst 2006). Thus, the elevated ratios of CH 3 OCHO and CH 3 OCH 3 to CH 3 OH hint a higher abundance of HCO and CH 3 , a longer time of warm condition at high T b,cont sources for more efficient formation of CH 3 OCHO, or more efficient formation of CH 3 OCHO and CH 3 OCH 3 in the gas phase after the evaporation/desorption of parent species, such as CH 3 OH and CH 3 O, as suggested by Balucani et al. (2015). Higuchi et al. (2018) found a tentative trend where the sources closer to the edge of the cloud or the isolated sources show a higher ratio of C 2 H to CH 3 OH. Although our ALMA observations cannot fully sample the emission of C 2 H, we can test if a similar trend exists between the abundance of COMs and the location of the sources in the cloud. We follow the approach described in Higuchi et al. (2018) to calculate the minimum distance (D min ) from the source to the edge of the cloud that is arbitrarily defined as the 10σ level in the Planck the C 2 H emission to confirm that will be presented in future papers.

The Abundance Ratios of O-bearing and N-bearing COMs
The abundance ratios of COMs reflect the chemistry of COMs. Figure 17 compares the abundance ratios of COMs toward the PEACHES sample to the ratios from the CALYPSO survey (Belloche et al. 2020), the observations of individual protostars, and the model pre-dictions from Garrod (2013). Belloche et al. (2020) divide up the CALYPSO sample into three groups, where Group 1 has low abundance of O-bearing COMs to CH 3 OH; Group 2 has high abundance ratio of O-bearing COMs to CH 3 OH compared to Group 1; and Group 3 is similar to Group 2 but with a higher abundance of CH 3 CN and CH 3 CH 2 CN relative to CH 3 OH and a lower abundance of O-bearing COMs. The abundance ratios of COMs toward the PEACHES sample generally agree with the ratios from the CALYPSO survey, which  is expected due to the significantly overlapped sample. Three out of four sources in Group 2 are Perseus protostars. If we compared to one of the archetype hot corinos, IRAS 16293−2422 B, the ratios of COMs in the PEACHES sample are 2-18 times higher (Jørgensen et al. 2016(Jørgensen et al. , 2018Calcutt et al. 2018). However, our sample has only a few sources with a comparable column density of CH 3 OH, 10 18 cm −2 , as that in IRAS 16293−2422, making the disagreement less robust. Also, our observations may underestimate the column densities of CH 3 OH due to the dust optical depth. Constraining the true difference of the ratios requires future analyses on the effect of dust opacity as demonstrated by De Simone et al. (2020).
We also compare the ratios of COMs to CH 3 OH with the numerical results of chemo-physical models. Garrod (2013) modeled three different timescale for the warmup phase in their model, where the formation of COMs become efficient (Garrod et al. 2008). The ratios of COMs toward the PEACHES sample differ from the ranges of the peak abundance ratios in the warm-up models by Garrod (2013). The models underestimate the abundance ratios for most of the O-bearing COMs, while the ratios of CH 2 OHCHO and N-bearing COMs agree with the ratios derived from the PEACHES sample ( Figure 17). Since we may underestimate the column density of CH 3 OH due to the optical depth, the abundance ratio to CH 3 OH could be overestimated. If the ratios derived from the PEACHES sample reduce by a factor of 30, they would be consistent with the ratios in Garrod (2013) for the five COMs in the left but the model would then overestimate the abundance ratios of CH 2 OHCHO, NH 2 CHO, and CH 3 CH 2 CN. This disagreement may point to a dominant gas-phase formation of (some) COMs other than CH 3 OH as suggested by other studies (e.g., Balucani et al. 2015;Skouteris et al. 2018;Vazart et al. 2020). The well correlated abundance of N-bearing COMs and O-bearing COMs (Figure 12 and 13) indicates no "N-/O-bearing COMs differentiation," which has long been discussed for Orion-KL, a massive star-forming region (e.g., Guélin et al. 2008;Friedel & Snyder 2008;Feng et al. 2015;Tercero et al. 2018). The N-bearing COMs and O-bearing COMs segregate toward the hot core and the compact ridge of Orion-KL, respectively, which has been interpreted as a result of different temperatures at the prestellar phase or different chemical evolutionary stage (Caselli et al. 1993;Neill et al. 2011;Laas et al. 2011;Garrod 2013). If the Perseus protostars have different initial conditions at the prestellar stage or are in different chemical stages that would result in a similar differentiation in Orion-KL, we would expect to detect a variety of N-bearing COMs to O-bearing COMs ratios. However, the well correlated abundance of Nbearing COMs and O-bearing COMs suggests that the Perseus protostars may form from a similar initial condition or be in a similar chemical stage. Further studies of region-to-region comparisons await.

Origin of the CH 3 OH/CH 3 CN Correlation
The tight correlation between CH 3 OH and CH 3 CN is a striking finding of the PEACHES survey, while studies, such as Bergner et al. (2017) and Belloche et al. (2020), have shown a similar trend with fewer detections or larger scatter. Belloche et al. (2020) argued that this correlation between CH 3 OH and CH 3 CN may not be due to chemistry because of the somewhat unrelated formation pathways of two molecules. If the gas-phase chemistry is negligible for the production of CH 3 OH and CH 3 CN, the tight correlation between CH 3 OH and CH 3 CN suggests a similar abundance ratio of two molecules on the icy grains prior to being thermally desorbed, hinting a common chemistry of COMs in the warm protostellar envelope. The tight correlation may also reflect a uniform elemental abundance of O and N in the Perseus molecular cloud before the star formation.
The compact emission of CH 3 OH has a well known icy origin, where CH 3 OH forms via subsequent CO hydrogenation on the surface of dust grains (Tielens & Whittet 1997;Watanabe et al. 2003;Rimola et al. 2014). At cold temperature, the non-thermal desorption and/or the excessive energy from the formation of CH 3 OH that kicks it off from dust grains are thought to be responsible for the presence of gaseous CH 3 OH (e.g., Garrod et al. 2007;Vasyunin & Herbst 2013;Soma et al. 2015, ; but see the discussion in Pantaleone et al. (2020) on the chemical desorption probability). CH 3 CN can form in gas-phase via radiative association (HCN + CH + 3 → CH 3 CNH + + hν) followed by the dissociative recombination of CH 3 CNH + or on the grain surface. The grain-surface reactions to form CH 3 CN can occur at cold temperature by the hydrogenation of CCN, which can form via C+CN or C 2 +N. At warm temperature, CH 3 CN can form directly from reactions between radicals, CH 3 and CN. For embedded protostars, the grainsurface reactions are thought to dominate the production of CH 3 CN; however, it remains unclear whether the cold phase or the warm phase grain-surface pathway regulates the abundance of CH 3 CN.
If we assume an icy grain origin for both CH 3 OH and CH 3 CN, we can test the efficiency of the cold and warm phase pathways by comparing the ratio of CH 3 CN/CH 3 OH toward the prestellar and protostellar sources. In prestellar cores, [CH 3 CN/CH 3 OH] are 0.28, 0.0192±0.0075, and 0.021±0.008 for TMC 1-CP (cyanoployyne peak) (Kaifu et al. 2004;Gratier et   The ratio of TMC1-CP is substantially higher than that in our survey, 0.012 +0.009 −0.005 , while the mean ratios of L1521 E and L1544 seem higher than the ratio in our survey within 1σ uncertainty. Furthermore, we can divide the PEACHES sample into Class 0 and I according to their T bol (Class 0: T bol < 70 K; Class I: T bol ≥70 K; Dunham et al. 2014). Only one Class I source has emission of both CH 3 OH and CH 3 CN along with another borderline Class I source, which has a T bol of 69 K. The mean value of CH 3 CN/CH 3 OH for Class I sources (N=2) is 6.5 +6.0 −3.0 × 10 −3 , while the mean value for Class 0 sources (N=12) is 1.4 +1.1 −1.1 × 10 −2 . Only the sources with reliable measurements of T bol were included. Although the difference is within the uncertainty, the tentative decrease of CH 3 CN/CH 3 OH from Class 0 to Class I phase would be consistent with the tentative decrease of CH 3 CN/CH 3 OH from prestellar to protostellar phase. Such difference provides constraints on the interplay between the production of CH 3 CN in the warm and cold phase. Testing the chemical pathways of CH 3 CN requires more measurements of CH 3 CN and CH 3 OH in Class I and prestellar sources.
Alternatively, the desorption process may regulate the amount of COMs in gas phase. Studies find a similar binding energy for CH 3 OH and CH 3 CN (Collings et al. 2004;McElroy et al. 2013;Wakelam et al. 2015Wakelam et al. , 2017Penteado et al. 2017). Recently, Ferrero et al. (2020) show similar range of binding energies for CH 3 OH (3770-8618 K) and CH 3 CN (4745-7652 K) simulated with amorphous solid water. Thus, both molecules may simply desorb from grains at a similar rate if the ice mantles have not completely desorbed. However, the microphysics during the desorption involve the structure of the water ice and how the COMs reside in the ice matrix, requiring further understanding of how the chemical evolution takes place on the grain surface. Also, the distribution of the binding energy can have crucial effect on the estimate of the desorbed species (Ferrero et al. 2020;Grassi et al. 2020).

CONCLUSIONS
This work presents the PEACHES survey, an unbiased chemistry survey of COMs toward 50 embedded protostars at the Perseus molecular cloud using ALMA Band 6 observations. The main conclusions are:

Perseus Class 0/I protostars commonly (58%)
show the warm COMs emission. CH 3 OH is the most detected species of COMs, while our observations show the emission of twelve O-bearing COMs and four N-bearing COMs, including isotopologues.
2. Complex organic molecules other than CH 3 OH are detected for the first time toward Per-emb 35 A and Per-emb 11 A.
3. The CH 3 OH and CH 3 CN normalized column density (N /T b,cont ) has large diversity (more than 2 orders of magnitude), suggesting the existence of chemical diversity among the sources in Perseus molecular cloud. Moreover, the detectability of those COMs is not directly related to the protostellar properties, such as T bol , L bol , and the gas mass column density.
4. The column density of CH 3 OH tightly correlates with that of CH 3 CN with a Pearson's r of 0.87, where two molecules may not have a common chemical origin. The abundance of the most frequently detected COMs, CH 3 OH, CH 3 CN, CH 3 OCHO, and CH 3 OCH 3 , correlate with other species of COMs, which are detected in fewer sources in the PEACHES sample, hinting at common chemistry of COMs in the Perseus embedded protostars. Protostellar properties, such as L bol , the evolutionary indicator, such as T bol , and the minimum distance to the cloud edge, have little impact on the number of detected COMs and the abundance of COMs. 5. CH 3 OCHO and CH 3 OCH 3 are the most abundant COMs other than CH 3 OH with abundance ratios to CH 3 OH of 0.22 +0.17 −0.10 and 0.21 +0.15 −0.09 , respectively. The abundance of COMs to CH 3 OH in the PEACHES sample are similar to that in the CALYPSO survey (Belloche et al. 2020), which includes 14 Perseus sources among 26 sample and 6 of 14 sources exhibit emission of COMs. The abundance ratios of NH 2 CHO and CH 3 CH 2 CN agree with the model predictions by Garrod (2013), which underestimates the abundance of the Obearing COMs. However, the abundance ratios of O-bearing COMs and CH 3 CN (e.g., CH 3 CN to CH 3 OCHO and CH 3 OCHO to CH 3 OCH 3 ) seem to be consistent with the ratios of abundance predicted in Garrod (2013) (Figure 17). alternatively, the disagreement may indicate that at least some COMs are the products of gas-phase chemistry (Balucani et al. 2015;Skouteris et al. 2018;Vazart et al. 2020).
6. The CH 3 CN/CH 3 OH ratio of Perseus sources seems lower than the same ratio for other prestellar cores, especially L1544. For the two Class I sources with detections of both molecules, the ratio is 6.5 +6.0 −3.0 × 10 −3 compared to 1.4 +1.1 −1.1 × 10 −2 for the Class 0 sources. This may suggest an evolution of the contributions from the gas-phase chemistry and grain-surface chemistry from prestellar to protostellar phase.
7. The ratios of more complex COMs, such as CH 3 OCHO and CH 3 OCH 3 , to the less complex COMs, such as CH 3 OH and CH 3 CN, increase with the averaged continuum brightness temperature, a proxy of the gas column density, suggesting an enhanced production of more complex COMs toward the sources with more centrally concentrated mass.   a The column densities of molecules are unconstrained due to the exclusion of the spectral window contaminated by the SiO emission.
Note-The column density is fitted with an aperture of 0. 5. a The column densities of molecules are unconstrained due to the exclusion of the spectral window contaminated by the SiO emission.
Note-The column density is fitted with an aperture of 0. 5.

APPENDIX
A. INTENSITY MAPS OF CH 3 OH, CH 3 CN, CH 3 OCHO, AND CH 3 OCH 3 Figure 18 shows the intensity maps of CH 3 OH, CH 3 CN, CH 3 OCHO, and CH 3 OCH 3 for the PEACHES sources not shown in Figure 3.

B. CATALOGS FOR MOLECULAR DATA
The spectroscopic data are taken from the Cologne Database of Molecular Spectroscopy (CDMS; Müller et al. 2001Müller et al. , 2005Endres et al. 2016) and the Jet Propulsion Laboratory (JPL; Pickett et al. 1998). Here we list only the references that cover the frequencies relevant to this study.  Table C2 continued  Table C2 continued  Table C2 continued

D. ESTIMATION OF THE COLUMN AND VOLUME GAS DENSITY
The continuum brightness temperature indirectly measures the gas column density assuming an optically thin emission. Within the extraction region for the 1D spectra, we can estimate the gas column density from the averaged continuum brightness temperature using where R g2d is the gas-to-dust mass ratio of 100, < I ν > is the averaged intensity, κ ν is the dust opacity, B ν (T d ) is the Planck function at T d , which is assumed as 30 K, and µ is the mean molecular weight of 2.37. We assume a frequency of 250 GHz. Furthermore, we can estimate the volume density assuming the emitting gas has the same spatial extent along the line of sight as that along the plane of sky. Figure 19 shows the estimated column and volume gas densities of the PEACHES sample. ] Figure 19. The derived gas column and volume densities for the PEACHES sample sorted by their gas column density.
blended transitions, making their excitation temperature unconstrained. Thus we assume excitation temperatures ranging from 100 to 300 K in our spectral modeling (Section 3.4). Fortunately, a few molecules have several transitions detected in our observations, such as CH 3 OH and CH 3 OCHO, allowing us to verify the assumption of the excitation temperature. Although the difference in the spectral coverage of the continuum spectral window results in different combinations of transitions, our spectra cover three CH 3 OH transitions, which have their upper energy ranging from ∼50 K to ∼500 K, allowing us to estimate the rotational temperature (T rot ) of CH 3 OH as a proxy of the excitation temperature (T ex ). To construct the CH 3 OH rotational diagram, we fitted the CH 3 OH emission with a Gaussian profile; then, we estimated the T rot by fitting a linear line to the rotational diagram (e.g., Goldsmith & Langer 1999). We further assumed that the E-and Aspecies of CH 3 OH have the same rotational temperature so that a single temperature component can describe the rotational diagram where the statistical weights due to the nuclear spins were considered. This derivation assumes no optical depth correction because we only have three transitions to fit a straight line, which has two parameters. The CH 3 OH emission at 243915 MHz may be affected by optical depth, especially in Per-emb 44, Per-emb 12 B, Per-emb 18, and Per-emb 29. Therefore, we excluded the emission between 243914 MHz and 342918 MHz in the Gaussian fitting to minimize the impact of the optical depth. Figure 20 shows the rotational diagram of Per-emb 22 B along with the fitted rotational temperature. The derived rotational temperature of CH 3 OH ranges from 140 K to 350 K, consistent with our assumption of T ex for the spectral modeling (Table E3).  Among the identified COMs, CH 3 OCHO has more than ten transitions detected. Thus, we can use the spectra of CH 3 OCHO to test our assumption of T ex . With many transitions detected, we can employ a more complex model to constrain the column density and the T ex of CH 3 OCHO. We used the xclass model described in Section 3.4, which includes the effect of optical depth. The model was optimized with the MCMC method (see description in Section 3.4), instead of a combination of genetic and Levenberg-Marquardt χ 2 minimization, to characterize the uncertainties of the fitted properties. Table E4 lists the sources which have sufficient detections of CH 3 OCHO and the derived properties of CH 3 OCHO. The best-fitting temperatures range from 100-300 K, consistent with our assumption of the excitation temperatures. Figure 21 shows an example of the posterior distributions of the fitting parameters toward B1-b S. The posterior distribution of the excitation temperature has a longer tail toward high temperature, an upper averaged uncertainty of 71.6 K compared to a lower averaged uncertainty of 29.3 K.
F. THE SPECTRA OF C 2 H The C 2 H spectra toward the continuum emission have irregular line profiles. Some spectra have strong selfabsorption, while some spectra only show the blueshifted emission. Due to the absorption and irregular line profile, the xclass fitting routine often fails to faithfully reproduce the observed C 2 H spectra. C 2 H can eas- ily form in the outflow cavity wall due to the abundant CH 4 sublimated from dust grains as well as C + ionized by the UV radiation. Thus, the C 2 H spectra can have broad line width and multiple components. Furthermore, the morphology of the C 2 H emission traces the outflows, making our extraction from the continuum emission non-ideal for representing the nature of the C 2 H emission. Figure 22 shows the spectra of C 2 H toward the PEACHES sample.

G. THE EFFECT OF SOURCE SIZE AND BEAM DILUTION ON COLUMN DENSITIES
Our modeling assumes a fixed 0. 5 source size to estimate the column densities because most emission of COMs is unresolved or marginally resolved in our observations. To study the impact of this assumption on the resulting correlation and scatter in the derived column densities, we tested whether the correlations shown in Section 4.2.2 still persist in two extreme cases, the column density averaged over the size of the continuum source and assuming COMs strictly emitting from the T =100 K radius. For the case of the averaged column density, we simply multiplied the modeled column density (N 0. 5 ) by Ω 0. 5 /Ω cont. , where Ω is the solid angle. Figure 23 shows the averaged column densities of the four most frequently detected COMs. The scatter in the averaged column densities is about 1-3 orders of magnitudes, similar to that in the modeled column densities assuming a size of 0. 5. If we assume that COMs

12
a The fitting uncertainty is smaller than the calibration uncertainty of 10%. Thus, a 10% uncertainty should be adopted for the fitted column densities.
b When the fitting uncertainty is smaller than the channel width, an uncertainty of the channel width, ∼0.1 km s −1 , is adopted.
emit from the T =100 K radius where water ice desorbs, we need to multiply the modeled column density by Ω 0. 5 /Ω 100 K . Using the equation for hot cores in Bisschop et al. (2007), we have Thus, Ω 0. 5 /Ω 100 K is proportional to L −1 . Figure 24 shows the modeled column densities of the four most fre-quently detected COMs normalized by L bol . The scatter of the normalized column densities is still about 1-3 orders of magnitude, suggesting that the intrinsic variation in the amount of COMs dominates the scatter in the column density. We also note that the L bol may not represent the protostellar luminosity due to the bias of the inclination of the protostars.     Figure 24. Corner plot of the correlations of the column densities normalized by L bol between CH3OH, CH3CN, CH3OCHO, and CH3OCH3. The legends and the color code follow that in Figure 11 and 12.