Complex organic molecules detected in twelve high-mass star-forming regions with Atacama Large Millimeter/submillimeter Array (ALMA)

Recent astrochemical models and experiments have explained that complex organic molecules (COMs; molecules composed of six or more atoms) are produced on the dust grain mantles in cold and dense gas in prestellar cores. However, the detailed chemical processes and the roles of physical conditions on chemistry are still far from understood. To address these questions, we investigated twelve high-mass star-forming regions using the ALMA band 6 observations. They are associated with 44/95GHz class I and 6.7 GHz class II CH$_{3}$OH masers, indicative of undergoing active accretion. We found 28 hot cores with COMs emission among 68 continuum peaks at 1.3 mm and specified 10 hot cores associated with 6.7 GHz class II CH$_{3}$OH masers. Up to 19 COMs are identified including oxygen- and nitrogen-bearing molecules and their isotopologues in cores. The derived abundances show a good agreement with those from other low- and high-mass star-forming regions, implying that the COMs chemistry is predominantly set by the ice chemistry in the prestellar core stage. One clear trend is that the COMs detection rate steeply grows with the gas column density, which can be attributed to the efficient formation of COMs in dense cores. In addition, cores associated with a 6.7 GHz class II CH$_{3}$OH maser tend to be enriched with COMs. Finally, our results suggest that the enhanced abundances of several molecules in our hot cores could be originated by the active accretion as well as different physical conditions of cores.


INTRODUCTION
Complex organic molecules (COMs; molecules composed of six or more atoms) are the building blocks of prebiotic molecules, such as amino acids. As such, one of the key questions on star formation is how COMs are synthesized and delivered to the planets and comets.
Recent astrochemical models and laboratory experiments have explained that COMs are produced on the icy surface and mantle of dust grain in cold (∼ 10 K) and dense (10 4 -10 7 cm −3 ) gas in star-forming regions, even in environments where the diffusion of large radicals cannot occur (Öberg & Bergin 2021). As a proto-giseon@khu.ac.kr jeongeun.lee@khu.ac.kr star forms, it warms up the innermost envelope, liberating the molecules from grain mantles into the gas phase when the grain surface temperature is high enough to evaporate the species. The dust grains are also heated by the shock energy generated by protostellar jets and outflows. As the protostar evolves, the temperature of the inner envelope reaches around 100 -300 K, forming a so-called hot core, or hot corino for low-mass counterpart, where the molecules are completely desorbed from dust grain mantles (Herbst & van Dishoeck 2009;Garrod 2013;Öberg et al. 2014;Öberg 2016;Garrod et al. 2022).
According to current models and laboratory studies, in cold and dense environments such as prestellar stage, radicals are produced on the icy mantle via successive atom-addition reactions or dissociation of present molecules. The production of COMs could occur in non-diffusive manners (Jin & Garrod 2020;Garrod et al. 2022), and if the mantle is warmed up (10 -100 K), radicals become mobile in the icy mantle, combine with other species and form COMs (Öberg 2016;Öberg & Bergin 2021). Once a protostar has formed, the protostellar luminosity determines the physical condition of its natal dense core, which in turn affects the further chemical characteristics.
High-mass young stellar objects (HMYSOs; M ≥ 8 M ) undergo a dense cold phase, and rapidly evolve (Motte et al. 2018). Since the production of the ice mantle on dust grain is sensitive to the density and timescale, the chemical richness and complexity depend on the initial conditions of ice formation in the early stage of each star-forming region .
If an HMYSO is formed via a scaled-up version of its low mass counterpart, its luminosity is mostly an accretion luminosity (Fischer et al. 2022). The degree of accretion has a significant effect on the temperature of a hot core region and the UV radiation strength. The temperature is one of the key factors that control the chemistry. The temperature regulates the synthesis and destruction rates of species, determining the abundances of COMs, even affecting subsequent chemical networks in both ice and gas phases. Strong UV radiation generated by a luminous protostar promotes photochemistry, causing abundant COMs, especially in icy grain mantles (Öberg 2016).
In spite of the expectations, the chemical characteristics in hot cores of HMYSOs are largely unknown because of their rapid formation timescales (∼ 10 5 years) and large distances (d > 1kpc), making structural identification very challenging. Since the development of facilities which enable to resolve the hot cores such as Atacama Large Millimeter/submillimeter Array (ALMA), tens of COMs have been found in hot cores in different physical conditions (Herbst & van Dishoeck 2009;Jørgensen et al. 2020;McGuire 2022).
On the other hand, the maser lines from various molecular species, such as OH, H 2 O, CH 3 OH, and SiO, are useful probes for the high-mass star forming regions to trace energetic processes in very specific conditions associated with HMYSOs. In particular, the CH 3 OH maser at 6.6685192 GHz (5 1 -6 0 A + ) is exclusively associated with high-mass star-forming regions (HMSFRs; Minier et al. 2003;Breen et al. 2013) since it is pumped by strong far-infrared emission from warm dust heated by HMYSOs (e.g. Cragg et al. 2005). VLBI observations have revealed that the CH 3 OH masers could trace circumstellar disks, torus, or outflows (Bartkiewicz et al. 2009(Bartkiewicz et al. , 2016Fujisawa et al. 2014). Although the timescales and sequence of the evolutionary phases of HMYSOs associated with various masers are still a matter of debate, the 6.7 GHz CH 3 OH maser tends to be associated with the early phase of high-mass starformation with active accretion (Walsh et al. 1998;Cyganowski et al. 2009;Breen et al. 2010;de Villiers et al. 2015).
Since the majority of previous maser studies have focused on kinematics and measuring distances, the physical and chemical characteristics combined with maser properties have not been sufficiently studied. The spatial resolutions to study thermal emission are so low that the information given by maser observation (e.g., Green et al. 2009;Hu et al. 2016) is diluted when it is compared with thermal line observation, making the comprehensive interpretation difficult. In addition, while it has become known that the formation of COMs is a natural consequence of star formation, their forming mechanism in such environment is still unclear. To address these questions, we aim to study the HMYSOs in different physical conditions by taking advantage of their strong emission and larger samples than their low mass counterparts, and with constraints provided by the class II CH 3 OH maser association.
We investigate the twelve HMSFRs observed with ALMA in previously unexplored spatial resolution, to trace how physical properties affect their chemistry, in particular, COMs. The twelve regions are the sources of KVN (Korean VLBI Network) and VERA (VLBI Exploration of Radio Astrometry) array (KaVA) Large program (Kim et al. 2018). They are associated with class I or II CH 3 OH masers (6.7 GHz; e.g., Green et al. 2010;Breen et al. 2015;Hu et al. 2016, and 44 and 95 GHz;Kang et al. 2016). The class I CH 3 OH maser associations hint indirectly the active mass accretion process in our targets because the class I CH 3 OH maser emerges in the shocked regions by outflows and jets located along various offsets from the central source. The presence of protostars is also endorsed by the class II CH 3 OH maser emission since it is radiatively pumped by the strong infrared radiation field. Therefore, their association implies that the active accretion process is ongoing in the regions. In particular, the class II CH 3 OH masers are a strong evidence of heated dust material, which can develop hot cores enriched with COMs.
One of our samples, G25.82-0.17 was revealed as a clustered HMSFR with sources in different evolutionary stages, from starless cores to an ultracompact H II region (Kim et al. 2020). In G25.82-0.17 W1 core, a rotating signature has been found at 230.368 GHz CH 3 OH thermal transition, and an outflow shock feature was also detected in the SiO line. The 229.759 GHz CH 3 OH transition shows a mixture of thermal and maser emission. Hirota et al. (submitted to PASJ) investigated G24.33+0.14 region. From miltiple transitions of CH 3 OH including the 229.759 GHz class I CH 3 OH maser and the lines of 13 COMs, they detected the variability of line intensity responding to the accretion burst. In the observed spectral windows, a great number of unidentified lines of multiple COM candidates were also detected.
In this paper, we present the chemical diversity of COMs in twelve HMSFRs using ALMA observation. We also explore the relation between the class II CH 3 OH maser and COMs chemistry in hot cores. In Section 2, the observation is described, and the results are presented in Section 3 including core detection, line identification, class II CH 3 OH maser search, and rotation diagram of COMs. In Section 4 we discuss the roles of physical properties and class II CH 3 OH maser on the characteristics of COMs, the correlation between COMs, isotopic ratios and comparison with other regions. Finally, we summarize the contents in Section 5.  Table 1. The data was calibrated using a CASA package (version 4.7;McMullin et al. 2007). Imaging was performed on each spectral window using the CASA 5.4.0 tclean task with a Briggs weighting of 0.5. The spectral resolution is ∼1.3 km s −1 and the angular resolutions (beam size) range from 0. 27×0. 24 to 0. 38×0. 25 for continuum images and from 0. 27×0. 24 to 0. 40×0. 28 for line cubes.
Due to the crowded emission lines in some regions, which makes the continuum level estimation using linefree channels challenging, the STATCONT method (Sánchez-Monge et al. 2018) was used for continuum subtraction. It statistically estimates the continuum level in each pixel to provide the continuum-subtracted cube (which only contains spectrum data) and pure continuum image. Note-a Position offset from the phase center to the strongest continuum peak of the target. b RMS values of primary beam uncorrected continuum images are presented to estimate the sensitivity of the observation. c Distance is assumed as the same as G18.34+1.78. d G23.43-0.18 is observed toward the G23.43-0.18 C2 for the pointing center, southern from the strongest continuum peak at C1. * mark the regions observed from offset of the phase center.

Dust Continuum
We present 1.3 mm dust continuum images for the twelve HMSFRs. For all regions, multiple continuum peaks were detected, and many of them have been resolved as multiple cores for the first time . Although a few regions (G10.32-0.26, G10.34-0.14, G28.37+0.07) were not located at the phase center of observation, the cores were detected with the signal to noise ratios high enough to investigate their properties in primary beam uncorrected images ( Table  2). As a result, a total of 68 continuum peaks were found in the 12 regions. We define the small regions that present their own continuum peaks as continuum cores, which cover the cores at different evolutionary stages from dense starless core, prestellar core, embedded protostellar cores at extremely early stage, hot cores with the COMs emission and cores associated with ultracompact HII region (Kim et al. 2020). Since this study focuses on the chemistry of hot cores, we identify cores associated with the COMs emission in Section 3.3. Table 2 summarizes the information of the detected cores and derived properties.
Using the primary beam corrected continuum image, the peak position and intensity were found to estimate the H 2 column density (N(H 2 )); it is given by is the dust continuum peak flux density in Jy beam −1 , Ω is the solid angle of the beam (Ω=πθ maj θ min /4ln2), µ H2 is the mean molecular weight of 2.8 (Kauffmann et al. 2008), and m H is the mass of the atomic hydrogen. B ν (T dust ) is the blackbody intensity at T dust . For dust mass opacity, κ ν =0.1(250µm/λ) β cm 2 g −1 (Hildebrand 1983) with β=1.85 (Ossenkopf & Henning 1994) was adopted. A gas to dust mass ratio of 100 is assumed.
The uncertainty of the derived N(H 2 ) comes from the dust optical depth and temperature. To derive N(H 2 ), the optically thin dust emission is assumed. In the Rayleigh-Jeans limit (hν/kT << 1; h and k are planck and boltzmann constants, and ν is frequency) at 1.3 mm (ν ∼230 GHz), a low optical depth is a valid assumption if the T dust is much higher than ∼10 K (T >> hν/k). We adopted T dust from the excitation temperature (T ex ) of CH 3 OH (T ex (CH 3 OH)>100 K; Section 3.4) at CH 3 OH peak position, which are offset from the continuum peak position. The actual temperature at the protostellar position would be higher than the derived T ex , decreasing N(H 2 ). For some cores, the suppressed or missing COMs emission is observed at the continuum peak (Section 3.3 and Section 4.6). This can be caused by the high dust optical depths at dust continuum peak position. In this case, the N(H 2 ) may be underestimated. Multi-frequency continuum observation with comparable angular resolution will better estimate T dust and N(H 2 ).
For measuring the source size and integrated flux, the 2D gaussian fitting was performed for the observed continuum peaks using the imf it task in CASA. A threshold of 5 σ of RMS was applied. Masking was applied to constrain the fitting size in each core in the clustered regions. For cores surrounded by the diffuse cloud material, a zero-level offset fitting was applied to subtract diffuse emission.
With the estimated source size and integrated flux density at 1.3 mm at each continuum core, the total enclosed mass was derived (M total =S ν D 2 /κ ν B ν (T dust ), Hildebrand 1983). S ν is the 1.3 mm integrated flux, and D is the adopted distance of each region presented in Table 1. The optically thin dust emission and a gas to dust mass ratio of 100 are assumed. The most uncertain parameter for mass estimation is T dust . Since the continuum observation in multiple frequencies at comparable spatial resolution does not exist so far, we alternatively estimate the lower and upper limits of the mass in each continuum core. For all continuum cores, we assume 20 K as the lowest average temperature in the enclosed mass, giving the upper limit for M total . For lower-limits, toward the COMs-detected cores (see Table 3 and Table 4), the derived T ex (CH 3 OH) was adopted. Since the T ex was derived at the CH 3 OH position which is slightly off from the protostellar position, the actual average temperature of entire continuum core would be lower than T ex (CH 3 OH). Thus, the derived M total can be considered as the lower limit. The COMsundetected continuum cores may include the extremely embedded protostars, but they are probably cooler to have lower T ex (CH 3 OH) than COMs-detected cores in Table 5. Thus, 100 K was adopted to give the lower limit of M total for COMs-undetected continuum cores, which may be gravitationally bounded but chemically younger than COMs-detected cores. The enclosed masses of many COMs-undetected continuum cores are within the ranges of low (M 2M ) or intermediate (2M < M 8M ) protostar. This is a natural consequence in the sense of the formation of YSOs in a cloud; highmass stars are often born in a filamentary structure. In addition, inhomogeneous density distribution in a turbulent massive core may result in additional formation of low and intermediate mass protostars even during the monolithic collapse of the high-mass star formation.
In-depth studies of G24.33+0.14 and G25.82-0.17 regions were conducted using the ALMA Cycle 3 data (Kim et al. 2020;Hirota et al., submitted to PASJ). The masses of G24.33+0.14 C1 -C3 are estimated as 54.4,8.7,and 4.0 M (Hirota et al.,submitted to PASJ), assuming T dust of 100 K. We adopted T ex (CH 3 OH) (254.5 K) for G24.33+0.14 C1 as a low limit temperature, and 20-100 K for C2 and C3. The mass lower limits for C2 and C3 (Table 2) are larger than those of Hirota et al., (submitted to PASJ) by a factor of two in 100 K. However, this difference is caused by the different dust opacity κ ν assumption by a factor of two (1.0 cm 2 g −1 in Hirota et al. (submitted to PASJ) and 0.047 cm 2 g −1 in this study). Kim et al. (2020) derived the masses of three cores in G25.82-0.17 (referred as W1-W3 in Kim et al. (2020), and C1-C3 in this study). In the same κ ν , the derived mass for G25.82-0.17 C1 in Kim et al. (2020) is 20-84 M assuming T dust of 75-300 K. We derived the mass lower limit of ∼26 M adopting T ex (CH 3 OH) of 225.3 K, but it will be 19.1 M if T dust is 300 K. For C2 and C3, Kim et al. (2020) derived mass upper limits which are higher (366 and 99 M in 15 K) than ours (260 and 60 M in 20 K) due to the different assump-tion of T dust . Overall, the range of masses derived in this study is consistent with previous studies given the absence of the exact temperature information. To accurately constrain the physical properties of each cores, multi-band observation is essential.

Millimeter Class II CH 3 OH Maser Candidate
The twelve high-mass star-forming regions covered by this study are associated with 44 and 95 GHz class I or 6.7 GHz class II CH 3 OH masers (Green et al. 2010;Matsumoto et al. 2014;Breen et al. 2015;Hu et al. 2016;Kang et al. 2016). Since each region has been resolved into several continuum cores in our ALMA observation (Section 3.1), we specified the cores that are directly associated with previously detected 6.7 GHz class II CH 3 OH maser sources by the Methanol Multibeam Survey (MMB; Green et al. 2010;Breen et al. 2015) and VLA (Hu et al. 2016) observations ( Figure B1- Figure  B12). The association is marked in Table 3.
Two CH 3 OH transitions, 1 1 -0 0 E at 213.427 GHz and 5 1 -4 2 E at 216.946 GHz, were theoretically predicted to be class II CH 3 OH masers (Cragg et al. 2005). Since our observation covers these two millimeter maser candidate transitions, we checked if the maser characteristics are detected in our data sets. We used three criteria to determine a maser: 1) narrower line widths near the protostellar position, 2) higher brightness temperature (T b ), and 3) more spot-like features in the channel map compared to those of thermal lines.
In our observation, both the 213.427 and 216.946 GHz CH 3 OH emission shows a combination of compact and extended emission structures. According to the masing mechanism, the compact structures are suspected to maser emitting candidate regions. The 213.427 GHz (1 1 -0 0 E) transition shows a similar peak intensity to the 216.946 GHz (5 1 -4 2 E) transition at compact structure but a stronger intensity at the extended structure. This might be due to more abundant population in a lower energy state. Thus, here we mainly investigate 216.946 GHz transition. We use the channel maps of the 216.946 GHz CH 3 OH emission in each region  to search for the class II CH 3 OH maser signature. A more detailed analysis and discussion on the G24.33+0.14 region are conducted in a separate paper (Hirota et al., submitted to PASJ).
Toward the cores emitting CH 3 OH emission, the emission peak coordinates were found. At this position, the N(H 2 ) derived by a 1.3 mm continuum image, the peak intensity of CH 3 OH emission, velocity and brightness temperatures were estimated. The information is summarized in Table 3. The estimated brightness tempera-tures of the 216.946 GHz CH 3 OH lines range up to ∼140 K, which are not exceptional for thermal emission.
Given that the maser emission emerges in a very compact region and extremely narrow frequency range, we cannot rule out a possibility that the maser emission can be significantly diluted in our beam and frequency bin, considering our insufficient spatial and spectral resolutions. Since the emission shows both the extended and compact structures, this transition also has the characteristics of the thermal and maser emission mixture. As a result, the 216.946 GHz CH 3 OH maser is not firmly confirmed but remains as a maser candidate. Higher spatial and spectral resolution observations will resolve the features and confirm their existence.
The 216.946 GHz CH 3 OH maser candidates are listed in Table 4. In Table 4, the cores marked with TH (Thermal) are cases where the emission is detected at the position. We determine the COMs association of a core when the morphology of the COMs emission is concentrated toward the core. In other cases, the COMs emission is detected at the position of a core but it is at the periphery of the extended emission associated with a nearby core.
In addition to the two Class II CH 3 OH maser candidate transitions, the Class I CH 3 OH maser at 229.759 GHz was also detected in many of our sources (Figure 2,). The class I CH 3 OH maser is pumped by molecular collisions in shocked region and has been observed toward several high-mass protostars (Slysh et al. 2002;Cyganowski et al. 2011Cyganowski et al. , 2012Hunter et al. 2014). In our samples, the 229.759 GHz CH 3 OH emission shows strong maser features; spot-like features were detected off from the protostellar position but coincident with the outflow shocked region (Kim et al. 2020 for G25.82-0.17;Hirota et al., submitted to PASJ for G24.33+0.14 ). The detailed analysis of the class I maser emission will be presented in a separated paper (Baek et al., in prep).

Line Identification
In the five spectral windows in Band 6, many emission lines were detected in the vicinity of cores with chemical diversity. Toward hot cores, it has been known that the emission is dominated by COMs with relatively common and strong detections of CH 3 OH, CH 3 OCH 3 , and HCOOCH 3 , of which their intensity and relative abundances differ by cores (Herbst & van Dishoeck 2009). Since our observation detected many other COMs, we define a COM-rich hot core if the core has more than 5 COMs, which include 13 CH 3 OH and CH 2 DOH in addition to the three COMs.   While the emitting regions of COMs emission are compact toward the continuum sources, their peaks are slightly offset from the continuum peak positions. The reduced intensity of COMs at the continuum peak position could be caused by two reasons: 1) high dust opacity, obscuring the COMs emission at the continuum peaks (Harsono et al. 2018;Lee et al. 2019b) and 2) chemical effects, destroying COMs inside the water snow line (Lee et al. 2020). In addition to the reduced intensity, kinematic motions, such as rotation and infall, are detected at the continuum peak of some cores, producing blended line profiles or inverse P-Cygni profiles (Section 4.6). As a result, it is difficult to extract correct chemical properties at the continuum peak position. To identify the COMs in a consistent manner, we extract the spectrum of each COMs-detected core at the peak position of the 216.946 GHz CH 3 OH emission, assuming that all COMs emission originates from the same region where the simplest COM, CH 3 OH, exists.
From a total of 68 cores, 28 cores (∼41%) are associated with COMs emission. Assuming that the 68 continuum peaks are the site of high-mass star-formation, the ratio gives hot core emergence relative to the dense cores. If the total lifetime of HMYSOs is ∼10 6 years, the hot core lifetime is estimated ∼4.5×10 5 years, slightly higher than previously derived statistical high mass protostar lifetime (Table 2 of Motte et al. 2018). This can be because we targeted the actively accreting cores traced by the class I and II CH 3 OH masers. Compared with low-mass counterparts, we derived slightly lower fraction, which is approximately 58% (Yang et al. 2021).
To identify the detected emission lines in the extracted spectrum of each core, the eXtended CASA Line Analysis Software Suite (XCLASS; Möller et al. 2017) is used. XCLASS is a toolbox for modeling the spectra. It estimates physical properties of each species considering gas and dust attenuation in LTE conditions by searching for the best-fit model for observed spectra using MAGIX optimization algorithms and molecular databases from laboratory experiments (CDMS and JPL).
To synthesize the spectra, XCLASS mainly requires five parameters; source size, T ex , column density of the molecule (cm −2 ), line velocity width (FHWM; km s −1 ), and velocity offset relative to V LSR (km s −1 ). In addition, hydrogen column density (N (H)=2×N (H 2 ) in cm −2 ) and dust property (dust mass opacity κ 1.3mm in cm 2 g −1 and spectral index β) can be defined to consider the continuum optical depth. This may diminish the line intensity in the optically thick dusty disk/inner envelope. For the N H used in XCLASS, we adopted 1.3 mm continuum intensity at the the CH 3 OH peak position, where the spectrum is extracted, assuming that T dust is the same as T ex (CH 3 OH) under the LTE condition. Since we extract the spectrum in a pixel at the CH 3 OH peak, which is close to the continuum peak, the gas and dust should be well coupled due to the high density. Therefore, the remaining free parameter to fit the observed spectra is the column density of a given molecule. For line widths and offsets, since the spectra is extracted offset from the protostar position in a rotating hot core, ranges of line width (3 -10 km s −1 ) and velocity offset from V LSR (-8 -+8 km s −1 ) are given to match the observed spectra.
To securely identify the lines, we followed the criteria of Jørgensen et al. (2020) using XCLASS; (1) synthesized spectrum of more than two transitions for the given species is systemically accounted for the observation, and (2) the predicted spectrum does not appear in any transitions that no observed signal is seen. An exception is NH 2 CHO (Formamide), of which one transition is covered in our spectrum; however, the strong intensity is expected and its deuterated species NH 2 CDO is identified for this species. We identify the molecules in 3σ detection, except for the following molecules which are detected as upper limits; CH 3 NCO (Methyl isocyanate) in G18.34+1.78SW C2 Figure 34), (CH 3 ) 2 CO (Acetone) in G23.43-0.18 C1 (Figure 39), CH 2 CHCN (Vinyl Cyanide) in G28.37+0.07 C2 (Figure 43), and CH 3 OCH 3 (Dimethyl ether) in G28.37+0.07 C5 (Figure 45). Due to the heavily blended spectrum in some cores, the best fit parameters are found by fitting two or three molecules together. Figure 2 shows the observed and fitted spectra of G25.82-0.17 C1. Many COMs, including several oxygen and nitrogen bearing molecules, are identified. The rotational transitions in vibrationally and torsionally excited states are also detected, and isotopologues of carbon, oxygen and nitrogen molecules and their deuterated species are found. The column density derived by the XCLASS fitting is summarized in Table 6. Figure  3 shows the abundances of identified COMs relative to CH 3 OH. Up to nineteen COMs are identified in a wide range of energy states (up to ∼800 K).
Generally, O-bearing COMs are more abundant than N-bearing COMs in terms of greater detection rates and higher abundances. The most frequently detected COM is CH 3 OCH 3 , found in all regions except for G29.91-0.03 C1. The next abundant species are, in turn, HCOOCH 3 (Methyl formate), C 2 H 5 OH (Ethanol), and CH 3 CHO (Acetaldehyde). For N-bearing molecules, besides CH 3 CN (Methyl Cyanide), of which only its carbon isotope is covered in our observation, NH 2 CHO, CH 3 CH 2 CN (Ethyl Cyanide) and CH 2 CHCN are abundant.
We assume the filling factor is 1 (source size > beam size). While this assumption is reasonable for most of COMs, some of the COMs and CH 3 OH maser candidate transitions might have smaller emission distributions than the beam size. Thus, the derived beamaveraged values would be lower limit or underestimated compared with the actual values.

Rotation diagram
To investigate the excitation conditions of COMs, we estimate the rotation temperature and column density toward the 216.946 GHz CH 3 OH emission peak position. The integrated intensity was measured by gaussian fitting. The molecular data (A ij , E u , g u ) is adopted from the Splatalogue 1 database (Table A1). The best fit and uncertainty are calculated using the bootstrap method, which conducts linear fitting 10,000 times with randomly generated integrated intensity considering observation and gaussian fitting errors.
The results are listed in Table 5 and Figure C13. The T ex (CH 3 OH) range from 110 to 260 K; plausible for hot cores. The derived T ex of 13 CH 3 OH tends to be lower than that of the main isotope.
We assume that the COMs coincide with CH 3 OH. However, the T ex derived for other COMs exhibit some scatters from that of CH 3 OH, of which many are higher than that of CH 3 OH ( Figure 4 and Table 5). Together, given that the emitting region is smaller in the integrated intensity map, the COMs emission could be originated from the hotter and closer to the protostar compared to CH 3 OH. Larger COMs tend to show higher T ex than CH 3 OH, implying that they have higher sublimation temperatures and trace the inner and hotter region. However, we note that the T ex of COMs are loosely constrained in the limited number of detected lines in our spectral coverage.

Detection rate of COMs
We find significant diversity in the detection rates of COMs among the CH 3 OH detected cores. Up to nineteen O-and N-bearing COMs were identified in each core, including rotational transitions in excited vibrational and tortional states. The isotopologues, including 13 C, 15 N, 18 O, and deuterated species, were also detected. We note that the number of detected COMs is not dependent on the RMS and SNR of the 1.3 mm 1 https://splatalogue.online/ continuum image (Table 1 and Table 2), supporting that the observation sensitivity and continuum level do not affect the COMs detection rate.
As shown in Figure 3, O-bearing COMs have higher detection rates compared to N-bearing COMs. Among the 28 cores at which we found CH 3 OH (representative O-bearing COMs), N-bearing COMs are detected in 16 cores (∼57 %). The most frequently detected molecules following CH 3 OH (detected in 28 cores) are CH 3 OCH 3 (27 cores) and HCOOCH 3 (25 cores). The least frequently detected molecules are CH 3 NH 2 (Methylamine; 8 cores), CH 3 NCO (Methyl isocyanate; 10 cores), and HCOCH 2 OH (Glycolaldehyde; 11 cores). Except for HCOCH 2 OH, these higher detection rates of O-bearing COMs imply that the formation network for O-bearing COMs is more efficient than that of N-bearing COMs.
Cores with only O-bearing species detected with no detection of N-bearing ones might be in an earlier stage of star formation, since nitrogen species are known to be formed via the surface and gas-phase chemistry at a longer timescale (Lee et al. 2020). If a core experiences a longer timescale in the warm-up phase prior to the sublimation of COMs into the gas-phase, the formation of more complex COMs, including N-bearing species, can be promoted (Garrod 2013;Garrod et al. 2022).
We note that the higher detection rate of O-bearing COMs could be caused by our spectrum extraction at the CH 3 OH peak positions. We assume that O-and N-bearing COMs are present in the same gas parcel, but this might not always be the case. Csengeri et al. (2019) showed the spatial differentiation between the two groups toward the spatially resolved envelope of high-mass core within the G328.2551-0.5321 clump, interpreting that the O-bearing COMs trace the shocked region caused by the infall in addition to the protostellar heating, while N-bearing COMs are more concentrated toward the protostellar position. Higher angularresolution and more sensitive observations are required to investigate the spatial distributions of the detected COMs and possible chemical differentiation. Figure 5 illustrates the relation of the number of detected COMs with N (H 2 ) and T ex (CH 3 OH). The COMs detection rate steeply grows with N (H 2 ) above N(H 2 ) > 2×10 24 cm −2 , implying that the formation of COMs is advanced in denser cores. This is consistent with the observational results of the low mass counterpart, hot corinos (Yang et al. 2021).
Although the cores with higher detection rates of COMs tend to have higher T ex (CH 3 OH), the COMs detection rate is less sensitive to T ex (CH 3 OH) than N (H 2 ). For cores with the T ex (CH 3 OH) higher than 200 K, no clear trend is found between the number of detected COMs and the T ex (CH 3 OH). This could be interpreted as the derived T ex is more related to the current physical condition rather than the conditions for the formation/sublimation of COMs. It has been suggested that COMs mixed in the ice mantles tend to be sublimated when the main constituent of bulk ice, such as H 2 O, evaporates (Öberg & Bergin 2021). A detailed experiment of Garrod (2013) demonstrates that the ice sublimation temperature largely depends on their bonding structures (CO-like or water-like species). However, the mixed ice with various constituents, including secondary or tertiary species to the dominant constituent water ice, can change the sublimation temperature (Burke & Brown 2015), diluting the feasibility of the current temperature as the sublimation temperature tracer of individual species. Instead, a higher number of COMs appear to be detected because in a higher T ex the molecular desorption tends to be promoted at larger radii of hot cores, which allows the COMs exist in gas-phase, resulting in increasing the detection rate.
Among 28 cores with COMs, ten are associated with 6.7 GHz class II CH 3 OH maser. We also checked the COMs detection for those cores. The 6.7 GHz class II CH 3 OH maser flux density does not linearly respond to the T ex (CH 3 OH), whereas 216.946 GHz intensity which is thought to be a mixture of the thermal and maser emission tends to have a correlation with the T ex (CH 3 OH) (Figure 6). However, the emerging condition of the 6.7 GHz class II CH 3 OH maser emission itself is very sensitive to the current physical conditions in both density and temperature (e.g., Burns et al. 2020). The 6.7 GHz CH 3 OH maser is also an indication for the evolutionary stages of the central objects. Survey observations have revealed that the sources detected in the 6.7 GHz CH 3 OH maser are more evolved than quiescent IRDC cores but less evolved than ultracompact HII regions (Paulson & Pandian 2022).
In our targets, the cores associated with 6.7 GHz class II CH 3 OH maser show a high COMs detection rate (Figure 7). As shown in Figure 5 and also discussed in Section 4.3, the number of detected COMs mainly depends on the density in the early stages. However, the thermal history and variation of UV radiation may also influence the COMs production and destruction rates. The 6.7 GHz class II CH 3 OH maser emerges in a region with the active accretion process, which elevates the temperature and produce more UV radiations. Therefore, the  . The relaltion between the number of detected COMs and N (H2) column density. Filled and unfilled symbols present 6.7 GHz CH3OH detected and undetected cores, respectively. Colors for the filled symbols depict the 6.7 GHz maser peak flux in VLA observation (Hu et al. 2016).
high COMs detection rates in 6.7 GHz class II CH 3 OH maser detected cores suggest that more COMs can be detected in actively accreting cores. This result implies that 6.7 GHz class II CH 3 OH maser could be used as a COMs-rich hot core tracer.
To sum up, the chemical diversity among hot cores could be originated by the degree of the accretion of each source as well as the different physical conditions of cores.

The effect of physical properties
We test the effect of physical parameters on the chemical evolution. Figure 8 and Figure D14 present the column density of each molecule (N(COMs)), while Figure   9 and Figure D15 show the COMs abundances with respect to the CH 3 OH (X(COMs)) along with N(H 2 ) and T ex (CH 3 OH), respectively.
We divide the cores into two types by the number of detected COMs (i.e., > 5 and < 5). For all cores with < 5 COMs, CH 3 OH, 13 CH 3 OH, CH 2 DOH, CH 3 OCH 3 , and HCOOCH 3 were detected, and for cores with > 5 COMs (COM-rich cores), up to 19 COMs were detected. We focus on the COM-rich cores to consistently investigate the evolution of COMs complexity in the range of corresponding physical conditions. In Figure 8, we measure the slopes of the N(COMs) with respective to the N(H 2 ) in the log-log space for COM-rich cores. For most of main isotopologues, the uptrend becomes slow down when log 10 N(H 2 ) > 24.7 cm −2 , indicative of the high dust optical depths, preventing the COMs emission from escaping the system, as well as the high line optical depths of the species. The cores with the largest N(H 2 ) are G24.33+0.14 C1 and G27.36-0.16 C1 (Table 3, Figure 41 and Figure 42), of which the COMs emission is significantly suppressed at their continuum peak positions. The N(COMs) at the CH 3 OH peaks might be also affected by the strong continuum emission, leading to a slightly flattened slope at the densest ends in the N(COMs) vs N(H 2 ) plots (Figure 8). An extreme case is G49.49-0.39 C3, which shows systematically lower column densities in all species compared to other COMrich cores (Figure 52). G49.49-0.39 C3 is located close to G49.49-0.39 clustered region with highest continuum levels (Table 2). Therefore, this core might be highly embedded in the cluster as well as its own envelope although the continuum emission is clearly resolved and makes a unique peak toward G49.49-0.39 C3. As a result, we derived the slopes of N(COMs) with respect to N(H 2 ) only using the COM-rich cores and log 10 N(H 2 ) < 24.7 cm −2 .
In the observed range of physical condition, for COMrich cores, the characteristics of COMs are divided into two groups. The first group includes CH 3 OH, CH 3 OCH 3 , HCOOCH 3 , CH 3 CHO, and CH 3 CH 2 CN, which were detected in most hot cores. The N(COMs) of this first group COMs are directly proportional to the N(H 2 ), and thus, X(COMs) are almost constant with the N(H 2 ). Those X(COMs) also show relatively flat distribution of their abundance with T ex (CH 3 OH) except for CH 3 CHO.
The second group includes (CH 3 ) 2 CO, HCOCH 2 OH, C 2 H 5 OH, and aGg'-(CH 2 OH) 2 for O-bearing COMs, and CH 3 NCO, CH 3 NH 2 , CH 2 CHCN, and NH 2 CHO for N-bearing COMs. They show steeper gradients in N(COMs) versus N(H 2 ) and proportional relation in Figure 8. Relation between the column densities of COMs and N(H2), which is estimated using 1.3 mm continuum images. Grey empty and black filled circles represent the cores with the number of detected COMs smaller and greater than 5, respectively. Purple solid line is the linear fit in log-log space for the COM-rich cores, except for an outlier, G49.49-0.39 C3 (skyblue). Red solid line is the same as the purple line but fitted only at log10N(H2) < 24.7 cm −2 , excluding the three densest cores. X(COMs) versus both N(H 2 ) and T ex (CH 3 OH), unlike the first group.
This difference found in the two groups can be caused by the different formation timescales. At a low N(H 2 ) (e.g., ∼10 24 cm −2 ), the second group species have low column densities, if detected, and grow with N(H 2 ) proportionally. The less efficient formation of the second group species than the first group species in the low density probably cause the difference in slope. Moreover, the efficient destruction of the first group species at high densities or the high optical depth of the first group species lines may flatten the slope of the first species group (Figure 8).
Depending on the formation routes in both ice-and gas-phases, some of the second group species are formed from the destruction of the first group species, resulting in shallow gradients for the first group species at the high density regime. Our scenario is also supported by the comparison between this study and PEACHES. Yang et al. (2021) reported the elevated ratios of HCOOCH 3 and CH 3 OCH 3 to CH 3 OH in increasing continuum brightness temperature which traces the total gas column density (Figure 15 in Yang et al. 2021), while those in this study show almost constant levels. Compared with the PEACHES samples, the low mass hot corinos in Perseus molecular cloud, the cores in this study would be much denser because they are mostly high mass hot cores associated with class II CH 3 OH maser (Table 2 and Table 4). Therefore, the differences between PEACHES survey and ours indicates that chem- Figure 9. Relation for the abundance ratios of COMs relative to CH3OH versus N(H2), which is estimated using 1.3 mm continuum images. Grey empty and black filled circles represent cores with the number of detected COMs smaller and greater than 5, respectively. Red solid line is the linear fit in log-log space for the COM-rich cores.
ical timescale becomes shorter in denser cores, thus the reaction could proceed to form the other molecules resulting in the fast production of the second group species from HCOOCH 3 and CH 3 OCH 3 .

Comparison with other sources
We compare our abundance estimates of COMs with those found in various physical environments from recent interferometric observations. Figure 10 shows the abundances of detected COMs relative to CH 3 OH. To show their spreads, the grey boxes indicate the abundance ranges of our samples. The dark grey horizontal lines are the median and shades are the middle 25 -75 % values.
The COMs detections are collected in various physical conditions, which are illustrated by the colored shapes in Figure 10. AFGL 4179 is taken as the high-mass hot core associated with an O-type star (Bøgelund et al. 2019), Sgr B2(N2) as COMs-rich hot core at the Galactic center Jørgensen et al. 2020), WB 89-789 for a chemically rich hot core with low metallicity located in extremely outer region of our galaxy, , Orion KL hot core and com-pact ridge regions as recently shocked regions (Tercero et al. 2018), IRAS 16293 B as a low mass hot conino in the Class 0 stage (Jørgensen et al. , 2020, CA-LYPSO survey results toward 11 Class 0 and I sources , PEACHES survey results for 28 Class 0 and I protostars with COMs in Perseus molecular cloud (Yang et al. 2021), HH 212 observation to trace the COMs in low mass disk atmosphere (Lee et al. 2019a), and V883 Ori to trace the COMs freshely evaporated from the disk midplane by the accretion burst (Lee et al. 2019b). The recently updated astrochemical models are also included, considering non-diffusive mechanism (Garrod et al. 2022). Note that abundances of the hot core WB 89-789 are multiplied by 4. According to Shimonishi et al. (2021), when the COMs abundances are supplemented by a factor of 4, they show great similarity with those of Galactic hot core, which can be accounted for its low metalicity located in outer Galaxy.
From our samples, we sort the O-and N-bearing molecules in order of median of abundance. Except for the lower abundance molecules found in Orion KL, which cause different slopes from C 2 H 5 OH to HCOCH 2 OH, overall trends agree within slightly more than an order of magnitude for both high and low mass sources.
This similarity is very striking when we consider that the physical conditions of collected sources vary significantly, 0.16-47 L and 0.12-32 L for hot coninos from CALYPSO ) and PEACHES (Yang et al. 2021), 3 L for IRAS 16293B (Jacobsen et al. 2018), 9 L for low-mass disk harboring protostar in HH 212 (Lee et al. 2019a), 400 L for bursting protostar V883 Ori (Lee et al. 2019b), and ∼2×10 5 L for hot cores AFGL 4167 and Sgr B2(N2) (Bøgelund et al. 2019;Bonfand et al. 2019). The fact that over an order of magnitude range in abundance detected in huge luminosity differences over six orders of magnitude indicates the current luminosity is not a dominant factor affecting the characteristics of COMs.
This similarity seems more robust in O-bearing COMs than N-bearing COMs, but we cannot confirm this due to the lack of the N-bearing COMs detection toward some regions (Orion KL, HH212 and V883 Ori). This consistency with respect to CH 3 OH abundance has been found in different star-forming regions, in various COMs van Gelder et al. 2020;Belloche et al. 2020;Jørgensen et al. 2020) and the astrochemical model . The overall trend of COMs abundance suggests that the physical conditions and production mechanisms required for COMs are analogous, and this COMs chemistry is predominantly set by the ice chemistry in the cold prestellar core stage. One additional scenario is hinted by Shimonishi et al. (2021). They compared the COMs abundances of a hot core in extremely outer Galaxy, Galactic hot core, and low metallicity hot core in Large Magellanic Cloud. The abundances between hot cores in outer Galaxy and inner Galactic region show a great similarity when the values of hot core in outer Galaxy are multiplied by a factor of 4, which corresponds to its poor metallicity level compared with near the solar neighborhood. In contrast, the abundances for the extreme hot core and that in Large Magellanic Cloud do not resemble, despite the similarity of the poor metallicity. Based on this, the found global COMs abundance agreement found in our collected samples may come from the coherent initial condition of our Galaxy which decides the overall COMs ratios.
On the other hand, in each species, scatters can reflect the roles of local conditions and subsequent physics, which control the desorption mechanisms and gas-phase chemistry following to the evaporation.
Our sources shows consistency with the CALYPSO group 3 samples which are classified as CN-COMs-rich hot corinos, PEACHES survey, and V883Ori for the detected COMs. One common characteristic for them is the burst accretion. Interestingly, the COMs emitting regions of CALYPSO group 3 sources are larger than their hot corino sizes with the temperature range of 100 K < T < 150 K estimated from their current luminosities (Figure 13 in Belloche et al. 2020). This is a strong indication of the prior burst accretion events that evaporates the COMs beyond the water snow line set by the current luminosity. For PEACHES, hot corino survey in Perseus molecular cloud known as one of the most active SFR, half (14 out of 28) of their sources have the evidence of prior or current burst accretion; the CO or H 2 O snowlines were observed in larger radii than those can be explained by the current bolometric luminosity (Hsieh et al. 2019;Yang et al. 2021). V883 Ori is a low mass protostar currently in the FUor ourburst phase with a bolometric luminosity of ∼400 L with an enhanced accretion rate by a few orders of magnitude than typical low mass protostars (Cieza et al. 2016). In addition, our samples include ten cores in active accretion, traced by the class II CH 3 OH maser association.
The systematically higher abundance compared with the hot core in AFGL 4176 is found in hot corino surveys (CALYPSO and PEACHES). It may indicate no distinctive correlation between the occurrence of COMs and the current luminosity, which is consistent with the previous results van Gelder et al. 2020;Yang et al. 2021). Instead, the strength of accretion events, which episodically varies during the star formation, may more affect the chemical characteristics.
The abundances of O-bearing COMs in Orion KL regions are exceptionally low. One can expect the density of Orion KL regions to be lower than other regions. It has been suggested that the COMs enhancement in Orion KL region is not caused by the internal heating source, but externally heated by a recent explosive shock within 550 years and/or outflow shocks driven by an HMYSO Orion source I (Zapata et al. 2011). In particular, the compact ridge region can be excited by the older outflows from Orion source I (Tercero et al. 2018). Unfortunately, the behavior of N-bearing COMs cannot be directly compared with our observation because of the absence of N-bearing COMs observations with high resolutions comparable to ours, toward Orion KL to date.
The high (CH 3 ) 2 CO abundance derived in our study stands out compared with other regions. This may in-dicate that the chemical process of (CH 3 ) 2 CO could be distinct from other O-bearing COMs. One suggestion is that (CH 3 ) 2 CO is more sensitive to radiation or shocks compared to other O-bearing COMs. Peng et al. (2013) observed the (CH 3 ) 2 CO toward Orion BN/KL, revealing that the distribution of (CH 3 ) 2 CO is closer to N-bearing COMs than other O-bearing COMs. They suggested that the formation or destruction routes of (CH 3 ) 2 CO can involve N-bearing species. This may be the case for our sources since the SiO emission, which traces the outflow shocks responding to recent burst accretions, is readily detected in many of our samples. The (CH 3 ) 2 CO abundance range of our sources covers the values of Sgr B2(N2) and V883 Ori, while the abundances are systematically higher than all other samples. Sgr B2(N2) is a well known hot core in active star-forming region with a strong radiation field in the Galactic center ). V883 Ori is in the actively accreting phase with FUor type outburst.
Meanwhile, for N-bearing COMs, the abundance of our sources are systematically lower than those found in hot core Sgr B2(N2) but higher than low mass hot corino IRAS 16293 B. One exception is CH 3 NCO. The abundance of our samples is comparable with Sgr B2(N2), while it is higher than that of IRAS 16293 B by more than an order of magnitude. Since N-bearing COMs are known to be sensitive to the strength of UV radiation field and/or shock, the abundance lower than that of Sgr B2(N2) and higher than that of IRAS 16293 B can be interpreted as intermediate strengths of UV radiation field and shock, induced by accretion and outflow, respectively. Therefore, the accretion rates of our sources might be in between the two regions. The other scenario is that our cores are chemically younger than Sgr B2(N2) but older than IRAS 16293 B. Given that most COMs formation occurs in the warm-up phase (T ∼ 20-40 K, Garrod 2013) and N-bearing molecules require a longer timescale than O-bearing molecules to be synthesized in the cold cloud (Lee et al. 2020), the ratio of O-bearing and N-bearing COMs can be a tracer of the timescale for the warm-up phase. Our samples might have experienced the warm-up phase timescales shorter than Sgr B2(N2) but longer than IRAS 16293 B.
The model in Garrod (2013), where the chemical complexity operates efficiently via a diffusive mechanism in the warm-up phase, shows discrepancy with our results for the O-bearing COMs. However, their updated model (Garrod et al. 2022), which considered non-diffusive mechanisms, shows better agreement with our observations. Figure 11. Heat maps of the correlation using the column densities (left) and abundances to N(H2) (right) between two species. The colors show the level of correlation which correspond to the Pearson r coefficient shown inside in each box. The scatter plots with measured slopes are presented in Figure E16 and E17.

Correlation between COMs
. We test the correlation of COMs to investigate how the molecules are chemically linked. The deduced molecular column densities span over 2 orders of magnitudes, which allows us to investigate their relative abundances in wide range of physical condition. The relative abundances of some species can provide the information on the branching ratio if they are indeed chemically related. Figure 11 displays the maps for correlation coefficients for all pairs of main isotopologues using the Pearson coefficient r. The correlation coefficient between two species in the column densities (left panel) and abundances relative to N(H 2 ) (right panel) ranges from 0.45 to 0.95 and from 0.49 to 0.95, respectively. For all pairs, the scatter plots with slopes measured in logarithmic scale are presented to infer the power-law index between two species ( Figure E16 and Figure E17). CH 3 OCH 3 and HCOOCH 3 , the most abundant species following to CH 3 OH in our observation, show a tight correlation. This is consistent with the previous observations in different physical conditions Jaber et al. 2014;Belloche et al. 2020;Yang et al. 2021), implying their strong relevance in various chemical pathways. Under the grain-surface chemistry, their close relation is explained by a common precursor. CH 3 OCH 3 and HCOOCH 3 are dominantly formed by the radical-radical reaction of CH 3 O with CH 3 and HCO radicals, respectively (Garrod 2013;Oberg 2016). The gas-phase process has been also proposed, for CH 3 OCH 3 being a precursor of HCOOCH 3 (Balucani et al. 2015). In our observation, both species also show great correlations with CH 3 OH ( Figure E16), supporting that the origin of CH 3 O radical is likely the photodissociation of CH 3 OH both in grain and gas (Garrod 2013; Balucani et al. 2015).
A strong correlation between C 2 H 5 OH and aGg'-(CH 2 OH) 2 is also explained by a shared precursor. On the grain-surface, C 2 H 5 OH, (CH 2 OH) 2 , and HCOCH 2 OH can be efficiently formed via the radicalradical reactions involving CH 2 OH radical, another product of the photodissociation of CH 3 OH (Vasyunin & Herbst 2013;Garrod 2013).
Moreover, in our analysis, C 2 H 5 OH (CH 3 CH 2 OH) tends to show high Pearson r coefficients (r > 0.8) with other O-bearing COMs containing CH 3 , indicative of a common radical involving their formation routes. Meanwhile, the coefficient in HCOCH 2 OH and aGg'-(CH 2 OH) 2 with others are relatively low. This can be due to the small samples with detection, or another nontrivial pathway which is less involved with the other molecules. A recent laboratory experiment showed that HCOCH 2 OH and (CH 2 OH) 2 can be formed via a sequential atom addition reaction on the ice mantle started from CO deposition . In their experiment, both species form with little CH 2 OH radical, and HCOCH 2 OH is produced with HCO radicals followed by hydrogenation and further forming (CH 2 OH) 2 . In addition to this formation network, gasphase formation route is also opened (Section 4.4.1).
The correlation coefficient between CH 2 CHCN and CH 3 CH 2 CN is not as high as others, due to the plateau observed toward the CH 3 CH 2 CN column density above 10 16 cm −2 , probably caused by the high optical depth of this molecule. According to Garrod et al. (2017), the appreciable amount of CH 2 CHCN is not achieved on the grain-surfaces because the successive atom addition reaction for CH 2 CHCN quickly proceeds to CH 3 CH 2 CN. Instead, after the ejection of CH 3 CH 2 CN from the dust grain, CH 2 CHCN primarily survives in gas due to the CH 3 CH 2 CN destruction. In our analysis (Section 4.2), we grouped the detected COMs into two by the response to the physical condition, probably due to their different formation timescales. CH 3 CH 2 CN belongs to the first group, of which the abundance relative to CH 3 OH maintains a constant level with respect to N(H 2 ) (Figure 9). This indicates the co-existence of CH 3 CH 2 CN and CH 3 OH, with similar formation timescale, primarily via grain-surface networks. In contrast, as belonging to the second group, CH 2 CHCN abundance shows positive gradient, indicating a relatively longer formation timescale of CH 2 CHCN than CH 3 CH 2 CN, which is consistent with the model prediction.
CH 3 NH 2 , CH 3 NCO, and NH 2 CHO contain peptide (-NH-(C=O)-) like bonds, which are regarded as potential compounds of amino acids. Thus, they are of interest to prebiotic chemistry because they are thought as precursors for more complex biomolecules. Their formation routes are still under debate, however, solid state reactions are suggested for their efficient formation starting with NH and NH 2 (Garrod et al. 2008;Garrod 2013;Barone et al. 2015;Kaňuchová et al. 2016;Belloche et al. 2017;Ligterink et al. 2018;Quénard et al. 2018;Gorai et al. 2021).
From our study, one of the observed characteristics of the three species is their abundances, which are in between those of IRAS 16293B and Sgr B2(N2) ( Figure  10); the difference may be attributed to the degree of UV radiation. In correlation plots (Figure 11), among three, CH 3 NCO and NH 2 CHO have the closest correlation, but general correlation is not robust compared to other pairs. Perhaps this is due to the relatively small number of samples with detection (∼10).
We note that the column densities of some abundant species may be underestimated by the high optical depths of their lines. A similar response to the physical parameters such as temperature, density, or IR radiation field can give a high correlation by chance . For the former issue, in our samples, the high optical depths of COMs seemingly affect only 3 out of 28 cores (Figure 8), thus they do not significantly affect the result. For the latter one, we distinguish a possible correlation using Figure D14. Last, the single linear fitting may not be appropriate to account for the relation between two species, as reported by El-Abd et al. (2019), which is not considered in this study. 4.4.1. CH3CHO and HCOCH2OH Figure 12. The correlation between the CH3CHO and HCOCH2OH of eleven cores where the two COMs are detected. The sold and dashed lines are the abundance ratio of 3.5 and 1.6, respectively, which imply the branching ratio suggested gas-phase formation route including their uncertainty (Vazart et al. 2020).
We compare our observations with theoretical computations. As one of the examples for key molecular species, we here compare correlation between HCOCH 2 OH and CH 3 CHO. Figure 12 shows that the correlation of column densities between CH 3 CHO and HCOCH 2 OH for eleven cores where both species were detected. The found proportional relation suggests that the two molecules may be chemically related. If they form from a sharing parent molecule, the slope of their proportional relation may provide their branching ratio. Skouteris et al. (2018) and Vazart et al. (2020) showed that the destruction of C 2 H 5 OH could efficiently form HCOCH 2 OH and CH 3 CHO with branching ratios of 3:1 in the gas-phase (solid and dashed lines in Figure 12 with the modeling uncertainty), which explains our data reasonably well. Skouteris et al. (2018) used the reaction coefficient measured at a sufficiently high temperature for the destruction of C 2 H 5 OH (100-300 K), and they reproduced the observation by using the gas-phase network followed by the sublimation of C 2 H 5 OH formed in ice. The uncertainty of their pathway is caused by the abundance of OH radicals for the destruction of C 2 H 5 OH at the first step of the gas-phase chemistry.
On the other hand, the recently updated model by Garrod et al. (2022) adopts non-diffusive grain-surface and ice-chemistry. Among three models with different warm-up timescales, the slow warm-up model better matches with our observation in COMs abundances (Figure 10) and ratio between CH 3 CHO and HCOCH 2 OH. This model more concerns the relavance of HCOCH 2 OH with HCOOCH 3 and CH 3 COOH (than with CH 3 CHO), as the three species are structural isomers. Nevertheless, the ratio between CH 3 CHO and HCOCH 2 OH in slow warm-up model is consistent with our data (∼4 ; Table 18 of Garrod et al. 2022).
Therefore, both ice-and gas-phases chemistry reproduce the observed ratio between HCOCH 2 OH and CH 3 CHO (e.g., Skouteris et al. 2018;Vazart et al. 2020;Garrod et al. 2022) in our samples. It suggests that the surface chemistry generally determine the COMs abundances, but the gas-phase production of aldehyde from the C 2 H 5 OH destruction also plays a role in hot cores.

Isotope ratio
The level of fractionations of COMs is sensitive to the density, temperature, and timescale so that they provide insight of when and where they are formed. There are a few mechanisms that can alter the isotopic ratios, including isotopic exchange reactions both in the gas-and ice-phases, selective photodissociations, and timescale spending at cold temperature (Furuya 2018;Jørgensen et al. 2020;Öberg & Bergin 2021). From our data sets, we examine the hydrogen, carbon, oxygen, nitrogen isotopic ratios.
The isotopic ratios of carbon, oxygen and nitrogen are presented as functions of N(H 2 ) and T ex (CH 3 OH). As discussed in Section 4.2, column densities of main isotope species toward three highest column density cores (log 10 N(H 2 ) > 24.7 cm −2 ; red dotts) might be underestimated due to the high optical depth, resulting in isotope ratios overestimated (upper panels of Figure 13). Except for the measurements in the three densest cores, no clear trend is found in isotopic ratios with respective to N(H 2 ). Instead, the isotopic ratios are tentatively proportional to T ex (CH 3 OH) (lower panels of Figure 13). Although we speculate the high optical depth for main isotopologues in the three densest cores, the detected linear trends are not significantly affected by the three cores.
The selective photodissociation in COMs unlikely cause the proportional trends since the abundances of COMs are very low compared to CO and N 2 , for which the self-shielding effect is important (Visser et al. 2009;Heays et al. 2014). More favored mechanism is the isotope change reactions. We found that the chemical processes become efficient in cores with high N(H 2 ) and T ex (CH 3 OH) (Figure 9 and Figure D15). Based on more notable correlation between the isotopic ratio and T ex (CH 3 OH), the efficient gas and surface chemistry, which could be induced by the high gas temperature and strong radiation field, may promote isotopic fractionation, as discussed toward IRAS 16293 B ). Nevertheless, the absolute isotope ratios of CH 3 OH in our samples (> 0.05) are much higher than those found in Sgr B2(N2) (0.04) and IRAS 16293B (0.015) Jørgensen et al. 2016).
We note that the levels of fractionation are much higher than the Galactic (Wilson & Rood 1994) and local interstellar medium (ISM) values (0.015 for 13 C/ 12 C; 0.0036 for 15 N/ 14 N; 0.0018 for 18 O/ 16 O; Nomura et al. 2022). Since we observed interior region, close to the massive protostar, these unusually high isotopic ratios may come from significantly different physical and chemical conditions deviated from the Galactic scale environment. In addition, the isotope fractionations differ by species, but the values adopted for the Galactic and local ISM ratios are derived from atoms, and thus, they could be very different from the ratios in COMs (e.g., Nomura et al. 2022;Drozdovskaya et al. 2022).
Deuterium fractionation has been used to constrain how long the cores stayed at cold temperature. Since D atom is twice heavier than H atom, thus stronger in bond interactions, the deuterated ions or molecules are favorably formed from HD (Öberg & Bergin 2021). The deuteration of molecules begins with two reactions: (1) H + 3 + HD ↔ H 2 D + + H 2 + 232 K and (2) CH + 3 + HD ↔ CH 2 D + + H 2 + 654 K. In turn, the produced ions become the main reactants to form simple molecules and COMs. Since the forward reactions are exothermic, the reverse reaction is inactive at cold temperature (< 30 K for H + 3 and < 300 K for CH + 3 ), which results in increasing D/H ratio if a core stays in the cold-phase long.
In Figure 14, we could not find a clear trend of the D/H ratios with respective to N(H 2 ) within the range of observed properties: the D/H ratio found in CH 3 OH tentatively proportional to T ex (CH 3 OH) and inversely proportional to N(H 2 ). We need to remind that T ex does not represent the temperature when COMs were formed. Instead, the current temperature could be a measure of the efficiency of atomic exchange reaction after desorption.
Compared with other regions, in most of our cores, the D/H ratios (0.001-0.13 for CH 3 OH) are higher than that of Sgr B2(N2) (<0.004 for O-and N-bearing COMs; Belloche et al. 2016). The low deuteration level of Sgr B2(N2) has been interpreted by the high temperature in prestellar stage in the extreme condition of Galactic Center region ). Among our samples, G49.49-0.39 C2 and G10.34-0.14 C3 show the lowest D/H ratios of 0.001 and 0.009, respectively ( Figure  51 and Figure 32). Except for G49.49-0.39 C4, cores in G49.49-0.39 tend to show low D/H ratios (≤0.02). Since deuterium fractionation increases if a core has stayed at the cold environment longer, the observed low ratios of  G10.34-0.14 C3 and cores in G49.49-0.39 region indicate that the cores mostly stayed in the high temperature environment through their formation history. The low deuteration level can be also caused by the shorter presellar core stage owing to more rapid collapse than other sources; in short timescale, the isotopic exchange reaction could not sufficiently occur (Faure et al. 2015).
This D/H ratio can be different by species. The fractionation of COMs tend to be higher than that of H 2 O (Persson et al. 2014;Taquet et al. 2019) because COMs would be formed on grain-surface after the CO freezeout on grain surfaces: H 2 O ice can be formed on grain surfaces even before the dense core is formed, while CO freeze-out occurs after the temperature drops below 18 K in cold and dense cores (Bisschop et al. 2006). Thus, ice structure can have two types of layers: inner H 2 Odominated ice layer with a lower D/H ratio and outer CO-based ice layer with a higher D/H ratio (Jensen et al. 2021). In addition, among COMs, the fractionation level also depends on the CO depletion timescale during core formation as well as core density (Bergin 2014;Öberg & Bergin 2021) since the frozen CO on the grain surfaces can participate in the ice chemistry producing COMs.
The D/H ratios of CH 3 OH and NH 2 CHO measured in IRAS 16293B are ∼0.02, lower than those of most our sources. The ratios of different COMs range over 0.005-0.2 in IRAS 16293B (Drozdovskaya et al. 2022). Within a source, the varied deuteration of individual species imply the different formation timescales of those molecules (Jørgensen et al. 2020). The D/H ratios in our sources are generally in agreement with the COMs ratios found in IRAS 16293B. However, some cores (G18.34+1.78SW C1, G18.34+1.78SW C2, and G18.34+1.78 C8 for CH 3 OH and G49.49-0.39 C3 for NH 2 CHO) show notably high D/H ratios. The D/H ratio can be enhanced when core stay longer at the dense and cold environment shielded from the interstellar UV radiation field before protostars form at the core center. In addition, the level of deuteration can also change with time via gas-phase reactions induced by the photodesorption.

COMs intensity at the continuum peak position
In our observation, the CH 3 OH emission lines are weakened at the continuum peak position (Figure 15). The suppressed CH 3 OH emission at the continuum peak position might be caused by the dust obscuration due to the high dust opacity near the protostar. This has often been seen in the ALMA view toward early stages of both low-and high-mass YSOs (e.g., Lee et al. 2019b;Ko et al. 2020;Kim et al. 2020). The observation toward low mass protostar TMC1A resolved the inner disk (Harsono et al. 2018) with 1.3 mm continuum and isotopologues of CO. The emission of CO isotopologues is suppressed toward the central position of the protostar. They are even completely devoid at the dust continuum peak. This observation is best reproduced by the exis-tence of CO emission at the central part, but obscured by the large (1 mm) size dust grains near the protostar.
This effect of dust obscuration is also investigated in detail in the NGC 1333 IRAS 4A binary system. The 4A1, which showed a stronger continuum emission than 4A2, is free of COMs emission whereas 4A2 is associated with the rich spectrum at the submillimeter wavelengths (De Simone et al. 2020). The polarization observation with ALMA (0.85-0.89 mm and 1.3 mm) and JVLA (6.3-16.7mm) showed a 90 • angle offset between the two observations only toward 4A1. This is explained by the dust grains with different optical depths seen in two wavelengths; the high opacity foreground material of 4A1 causes high extinction. This leads to observed orthogonal polarization angles between two frequencies.
In contrast to the ALMA observation, a similar brightness of CH 3 OH emission is observed toward the two sources in the centimeter regime, at which the dust opacity is almost negligible (De Simone et al. 2020).
This can be the case for some of our samples since they are in the early stage of star formation and embedded in thick envelope material. Different detection rates have been seen within a clustered region; a core with strong continuum emission shows weaker or no COMs emission, compared to other cores with weaker continuum emission. One of such examples is the G29.91-0.03 region (Figure 47 and Figure 48); C1 shows a higher N (H 2 ) at both the CH 3 OH peak and continuum peak positions than those of C2 by a factor of two. However, only the CH 3 OH emission was detected in C1, while three COMs were detected in C2. The other example is G18.34+1.78SW; 15 COMs were detected in C2, while only three COMs were detected in C1, where N (H 2 ) is higher than that of C2 by a factor of 1.3 (Figure 37, Figure 38 and Table 4).
Extended toward COMs non-detected cores with higher continuum peaks, the G28.37+0.07 region is particularly interesting. At G28.37+0.07 C1 the strongest continuum emission is observed, but no COM is detected (Table 3 and Table 4). In contrast, G28.37+0.07 C2 showed a rich spectra with 15 COMs (Figure 43), associated with the 6.7 GHz class II CH 3 OH maser.
However, the presence of any embedded protostar or gravitational boundness is not checked toward G28.37+0.07 C1.
If the G29.91-0.03 C1, G18.34+1.78SW C1, and G28.37+0.07 C1 cores are in the same situation with the aforementioned case toward NGC 1333 IRAS 4A system, they are still highly embedded and the density near the central protostar might be too high for the molecular emission to escape. As a result, to avoid the effect of thick dust material on the COMs lines, we extract the spectra from a position slightly off the continuum peak position.
To investigate the COMs richness and complexity without hinder of the high dust opacity, the more transparent windows at longer wavelengths can be useful to study the COMs emission close to the central protostar unless the non-thermal effect is high (McGuire et al. 2020).
Another scenario for the offset of COMs peak from the central position is additional desorption mechanisms for COMs. We speculate the thermal evaporation of COMs due to their compact spatial structure around the protostar. There are several observations showing that COMs can be also evaporated by shocks by outflowing winds originated very near the protostar (Lee et al. 2019a) or infalling material (Csengeri et al. 2019). However, the portion of detected COMs desorbed by the non-thermal mechanism cannot be separated because it requires a detailed physical and chemical modeling toward individual cores, which is beyond the scope of this work.

SUMMARY
We aim to study the chemical characteristics of COMs in hot cores and the roles of physical conditions on their chemistry. The results of our study is summarized as follows.
1. We investigate twelve high-mass star forming regions associated with class I and class II CH 3 OH masers, using the ALMA band 6 observations. In twelve regions, total 68 continuum peaks are resolved. Among them, COMs are detected in 28 hot cores, in which 10 cores are associated with 6.7 GHz class II CH 3 OH maser.
2. Up to 19 COMs are identified in a core, including rotational transitions in torsionally excited state. 13 C, 18 O, and 15 N isotopologues and deuterated species are also detected.
3. COMs detection rate increases in cores with high N(H 2 ). We divide COMs into two groups of which the first one shows flat ratio of their abundance relative to CH 3 OH with N(H 2 ), in contrast to positive slopes found in second group. Similar tendency is found with T ex but less sensitive than N(H 2 ). It indicates that the COMs complexity is efficiently elevated in denser and hotter core.
4. The cores associated with the class II CH 3 OH maser show higher COMs detection rates, suggesting that the 6.7 GHz class II CH 3 OH maser can be a COMs rich hot core tracer. 5. By the comparison of our sample with various types of COM-rich sources, we found a general similarity in the relative COMs abundances to CH 3 OH over six orders of magnitudes in source luminosity. It suggests that the COMs chemistry is not dominantly affected by the current luminosity, but is predominantly set by the ice chemistry in the cold prestellar core stage. On the other hand, different COMs abundances within an order of magnitude may be influenced by the local physical conditions. 6. The correlations among COMs detected in our observation are mostly explained by the current chemical models, however, the distinct roles of grain-surface and gas-phase chemistry are needed to be examined by models to understand what causes the relative abundances of COMs including CH 3 CHO and HCOCH 2 OH.

ACKNOWLEDGMENTS
We thank the anonymous referee for detailed and helpful comments. We thank J.-H. Kang for comments on the data reduction and manuscript. We Thank J. Kim and KaVA team for sharing the distance information toward a part of the regions. We also thank Y.-L. Yang for sharing the results of PEACHES survey and related discussions. This work was supported by the National Research Foundation of Korea (NRF) Table 2 continued Note-a SNR is measured in primary beam uncorrected images. b For N(H2), we adopt the Tex(CH3OH) which is derived from the rotation diagram of CH3OH for COMs-detected cores. Note that Tex(CH3OH) is derived at the CH3OH emission peak located off from the continuum peak where the high dust optical depths affects the line intensity. The actual Tex(CH3OH) would be higher than the adopted temperature. The T dust and Tgas are assumed to be the same. For COMs-undetected continuum peaks, a range of T dust temperature of 20 to 100 K are adopted for prestellar to protostellar cores, assuming that the cores without COMs emission are cooler than those with COMs emission. c Convolved source size. d The lower and upper limits are derived. For all components, the lowest average temperature is assumed as 20 K to provide the upper limit of the total mass. For COMs detected cores (Table 4), the Tex(CH3OH) is adopted ( Table 5). The actual average temperature of each core would be lower than the derived Tex(CH3OH), which is measured at the CH3OH peak position, thus the derived mass is the lower limit. For COMs undetected continuum components, the highest temperature of 100 K is assumed for mass lower limit.  Table 3 continued Note-a N(H2) is obtained by 1.3mm dust continuum image at which the spectra is extracted.   b Green et al. (2010); Breen et al. (2015) and Hu et al. (2016).

1.09±0.48E+14
· · · · · · · · · · · · · · · · · · G18.34+1.78SW 2 90.3±51.8 9.18±3.79E+13 · · · · · · · · · · · · · · · · · · G23.43-0.18 1 87.6±44.9 0.86±0.27E+14 · · · · · · · · · · · · · · · · · · G23.43-0.18 2 112.2±128.0 9.37±11.93E+13 · · · · · · · · · · · · · · · · · · G24.33+0.14 1 117.4±77.6 9.54±6.77E+14 · · · · · · 205.2±182.4 3.39±4.25E+17 · · · · · · Table 5 continued                   Figure B1 for the CH3OH emission peak position from which the spectra were extracted. Figure 30. The same as Figure 28 except for G10.34-0.14 C1. Check Figure B2 for the CH3OH emission peak position from which the spectra were extracted. Figure 31. The same as Figure 28 except for G10.34-0.14 C2. Check Figure B2 for the CH3OH emission peak position from which the spectra were extracted. Figure 32. The same as Figure 28 except for G10.34-0.14 C3. Check Figure B2 for the CH3OH emission peak position from which the spectra were extracted.  Figure B3 for the CH3OH emission peak position from which the spectra were extracted.  Figure B3 for the CH3OH emission peak position from which the spectra were extracted.  Figure B3 for the CH3OH emission peak position from which the spectra were extracted.  Figure B3 for the CH3OH emission peak position from which the spectra were extracted. Figure 37. The same as Figure 28 except for G18.34+1.78SW C1. Check Figure B4 for the CH3OH emission peak position from which the spectra were extracted.  Figure B4 for the CH3OH emission peak position from which the spectra were extracted.  Figure B5 for the CH3OH emission peak position from which the spectra were extracted.  Figure B5 for the CH3OH emission peak position from which the spectra were extracted. Figure 41. The same as Figure 28 except for G24.33+0.14 C1. Check Figure B6 for the CH3OH emission peak position from which the spectra were extracted.  Figure B8 for the CH3OH emission peak position from which the spectra were extracted. Figure 43. The same as Figure 28 except for G28.37+0.07 C2. Check Figure B9 for the CH3OH emission peak position from which the spectra were extracted.  Figure B9 for the CH3OH emission peak position from which the spectra were extracted. Figure 45. The same as Figure 28 except for G28.37+0.07 C5. Check Figure B9 for the CH3OH emission peak position from which the spectra were extracted. Figure 46. The same as Figure 28 except for G28.37+0.07 C6. Check Figure B9 for the CH3OH emission peak position from which the spectra were extracted.   Figure B10 for the CH3OH emission peak position from which the spectra were extracted.  Figure B11 for the CH3OH emission peak position from which the spectra were extracted.  Figure B12 for the CH3OH emission peak position from which the spectra were extracted.  Figure B12 for the CH3OH emission peak position from which the spectra were extracted.  Figure B12 for the CH3OH emission peak position from which the spectra were extracted.  Figure B12 for the CH3OH emission peak position from which the spectra were extracted.    Table A1 continued     Table 3. The complete figure set (12 images) is available in the online journal.           Figure E17 depict the pair plots of column density of COMs and their abundances relative to N(H 2 ) for all main species detected in this work. Pearson r coefficients are estimated in all pairs, and slopes are measured in log-log space to give power-law indices. The red solid lines are the fitted line with linear function with the corresponding slope. Figure E16. Correlation between the column density of COMs. Pearson r coefficient of each pair is presented. The slope measured in log-log space are presented as s to give the power-law indicies of two given species. Figure E17. Correlation between the abundances of COMs to N(H2). Pearson r coefficient of each pair is presented. The slope measured in log-log space are presented as s to give the power-law indicies of two given species.