The complex chemistry of hot cores in Sagittarius B2(N): Influence of cosmic-ray ionization and thermal history

As the number of complex organic molecules (COMs) detected in the interstellar medium increases, it becomes important to place meaningful constraints on the formation pathways of these species. The molecular cloud SgrB2(N) is host to several hot molecular cores in the early stage of star formation, where a great variety of COMs are detected in the gas phase. Because of its exposure to the extreme conditions of the the Galactic center region, SgrB2(N) is one of the best targets to study the impact of environmental conditions on the production of COMs. Our main goal is to characterize the physico-chemical evolution of SgrB2(N)'s sources in order to explain their chemical differences and constrain their environmental conditions. The chemical composition of SgrB2(N)'s hot cores, N2, N3, N4, and N5 is derived by modeling their 3mm emission spectra extracted from the EMoCA imaging spectral line survey performed with ALMA. We derive the density distribution in the envelope of the sources based on the masses computed from the ALMA dust continuum emission maps. We use the radiative transfer code RADMC-3D to compute temperature profiles based on the COM rotational temperatures derived from population diagrams. We use published results of 3D RMHD simulations of high-mass star formation to estimate the time evolution of the sources properties. We employ the chemical code MAGICKAL to compute time-dependent chemical abundances in the sources and investigate how physical properties and environmental conditions influence the production of COMs. We find that chemical models with a cosmic-ray ionization rate of 7e-16s-1 best reproduce the abundances with respect to methanol of ten COMs observed toward SgrB2(N2-N5). We also show that COMs still form efficiently on dust grains with minimum dust temperatures in the prestellar phase as high as 15K, but that minimum temperatures higher than 25K are excluded.


Introduction
To date more than 200 molecules 1 have been discovered in the interstellar medium (ISM) or in circumstellar envelopes of evolved stars (see, e.g., McGuire 2018). These molecules are mostly organic (that is they contain at least one atom of carbon) and among them about one third are composed of six or more atoms and are considered as complex in terms of astrochemistry (Herbst & van Dishoeck 2009). As the inventory of complex organic molecules (COMs) expands, they have attracted a lot of attention, in particular recently due to their possible link to amino acids, the chemical building blocks of proteins which contribute to life as we know it on Earth. Numerous extraterrestrial amino acids have been discovered in meteorites (Botta & Bada 2002) and one in comets (Altwegg et al. 2016). However, the simplest amino acid, glycine (NH 2 CH 2 COOH), has not 1 See https://cdms.astro.uni-koeln.de/classic/molecules yet been detected in the ISM. Nonetheless, many COMs have been discovered toward the warm and dense regions associated with high-mass star formation, known as hot cores, by investigating their rotational spectrum at (sub)millimeter wavelengths. Despite these numerous detections, the precise origins of these complex species, as well as the mechanisms leading to their formation, are still subject to debate. The rise of astrochemical models solving coupled kinetic rate equations built from networks of thousands of chemical reactions has allowed significant progress in our understanding of the complex chemistry in the ISM. For instance, it has been shown that dust-grain chemical processes play an important role in the formation of COMs (Garrod & Herbst 2006;Garrod et al. 2007;Garrod 2008;Herbst & van Dishoeck 2009), the abundances of some of which cannot currently be explained by gas-phase chemistry. Therefore, it is crucial to better understand the chemical processes that take place in the interstellar ices at early time and low temperature to A27, page 1 of 47 explain the later high gas-phase chemical abundances observed toward regions forming high-mass stars. The numerous COMs detected toward hot molecular cores, right after desorption from the dust grains into the gas phase, may be used to trace the physical conditions of their parent molecular cloud since they hold information on the collapse history.
The Galactic center (GC) cloud Sagittarius B2 (Sgr B2), located at ∼100 pc in projection from the central supermassive black hole Sgr A * , exhibits one of the richest molecular inventories observed to date and thus appears to be an excellent target to search for COMs and probe interstellar chemistry. Many of the first detections of interstellar molecules at radio and (sub)mm wavelengths were made toward Sgr B2 (Menten 2004), in particular some of the most complex species detected in the ISM, such as acetic acid CH 3 COOH (Mehringer et al. 1997), glycolaldehyde CH 2 (OH)CHO (Hollis et al. 2000), acetamide CH 3 C(O)NH 2 (Hollis et al. 2006), aminoacetonitrile NH 2 CH 2 CN (Belloche et al. 2008), and N-methylformamide CH 3 NHCHO . With a mass of ∼10 7 M in a diameter of ∼40 pc (Lis & Goldsmith 1990), Sgr B2 is also one of the most prominent regions forming high-mass stars in our Galaxy. Its two main centers of star-formation activity, Sgr B2(N)orth and Sgr B2(M)ain, both host a cluster of (ultra)compact HII regions (Mehringer et al. 1993;Gaume et al. 1995;De Pree et al. 1998, 2015 and several 6.7 GHz class II methanol masers (Caswell 1996), which uniquely mark the locations of massive young protostellar objects.
Given its exceptional characteristics and its strong, ongoing star-formation activity (e.g., Bonfand et al. 2017;Sánchez-Monge et al. 2017;Ginsburg et al. 2018), Sgr B2 stands out from the other clouds located in the central molecular zone (CMZ; inner 200-500 pc) of the Galaxy. Indeed, recent studies find that despite its large reservoir of dense gas, the CMZ appears to be deficient in star formation overall in comparison to other high-mass star-forming regions in the galactic disk with similar density and age (e.g., Rathborne et al. 2014;Kruijssen et al. 2014;Kauffmann et al. 2017;Walker et al. 2018). This apparent inactivity suggests that the critical density for star formation in the CMZ is ∼10 7 cm −3 (see, e.g., Rathborne et al. 2014;Kruijssen et al. 2014;Kauffmann et al. 2017;Ginsburg et al. 2018; Barnes et al. 2019), much higher than in the rest of the galactic disk (∼10 4 cm −3 , see, e.g., Lada et al. 2010Lada et al. , 2012. The environmental conditions in the CMZ are also known to be extreme compared to the rest of the galactic disk, in particular with a stronger interstellar radiation field (ISRF, Lis et al. 2001;Clark et al. 2013), and higher cosmic-ray fluxes (Oka et al. 2005;van der Tak et al. 2006;Yusef-Zadeh et al. 2007;Clark et al. 2013). The exceptional physical conditions of the CMZ provide an excellent test of current models of high-mass star formation and can advance our understanding of star-formation mechanisms within extreme environments. Due to its proximity to the GC, Sgr B2 provides us with an interesting case study to investigate how extreme physical and environmental factors influence the high-mass star formation process, whether they tend to enhance or suppress star formation activity.
In a previous analysis of the 3 mm imaging line survey Exploring Molecular Complexity with ALMA (EMoCA) we reported the detection of three new hot cores, Sgr B2(N3), N4, and N5 (Bonfand et al. 2017). These new hot cores are all associated with 6.7 GHz class II methanol masers and Sgr B2(N5) also with an ultra-compact HII (UCHII) region. From the analysis of the 3 mm emission spectra observed toward the new hot cores, we derived their chemical composition, revealing the presence of up to 23 different species of which about half are complex (Bonfand et al. 2017). Because the spectrum observed toward the main hot core Sgr B2(N1) is close to the confusion limit and severely affected by line blending, here we focus our analysis on the four other hot cores Sgr B2(N2-N5). That all of these hot cores most probably originate from the same cloud material, with similar chemical and physical properties (e.g., initial chemical abundances, gas-to-dust mass ratio), makes comparing them with each other particularly meaningful. By comparing the predictions of astrochemical models with the observations, we want to identify the most efficient pathways leading to the formation of COMs (see for instance Belloche et al. (2017) for the case of N-methylformamide) and explore the impact of the environment on the production of these complex species.
The extreme environmental conditions of the CMZ most likely also affect the chemistry of Sgr B2(N)'s hot cores. In particular, in dense regions shielded from direct UV irradiation, cosmic rays play an important role in the gas-phase chemistry as they represent the main source of ionization. The direct ionization of atomic and molecular hydrogen results in the formation of secondary electrons (see Eq. (1)), which can cause additional ionization (Goldsmith & Langer 1978). In addition, collisions of H and H 2 with secondary electrons inside dense clouds lead to the production of ultra-violet (UV) radiation (Prasad & Tarafdar 1983). This cosmic ray-induced UV field can dissociate molecules, both in the gas phase and on dust-grain surfaces inside the clouds, where external UV photons cannot penetrate.
The so-called cosmic-ray ionization rate (CRIR), ζ H 2 , describes the total rate of ionization per H 2 molecule per second, including secondary ionizations. The first theoretical determination of the CRIR was made by Hayakawa et al. (1961) and revised later by Spitzer & Tomasko (1968). The latter derived a lower limit of ζ H = 6.8 × 10 −18 s −1 , which corresponds to ζ H 2 ∼ 1.3 × 10 −17 s −1 according to the simple approximation ζ H 2 ∼ 2ζ H (Glassgold & Langer 1974) commonly used for dense gas in which most hydrogen is in its molecular form. This value of ζ H 2 is often referred to as the standard CRIR. An accurate measurement of the CRIR is difficult to obtain through direct detection because low-energy particles (<1 GeV) are prevented from entering the heliosphere by the solar wind, making them not directly measurable from Earth (Parker 1958). However, the CRIR can be estimated with astrochemical models, by analyzing specific molecules whose formation and destruction processes are driven by cosmic rays. Given its relatively simple chemistry, the H 3 + ion is commonly used as an indirect probe for constraining the CRIR: attempts to constrain the CRIR in diverse regions of the ISM, huge uncertainties still remain. Therefore, the CRIR is usually treated as a free parameter in astrochemical models, ranging from ζ H 2 ∼ 10 −17 s −1 to 10 −14 s −1 . The impact of variations in the CRIR value on calculated chemical abundances can then be used to identify functional probes for the CRIR in dense clouds and determine the ionization rate that best reproduces the observations (see, e.g., Albertsson et al. 2018;Allen et al. 2018).
In this paper we use the chemical kinetics code, MAG-ICKAL (Garrod 2013), coupled with the radiative transfer code RADMC-3D (Dullemond et al. 2012) and the published results of radiation-magnetohydrodynamical (RMHD) simulations of high-mass star formation (Peters et al. 2011) to model the physico-chemical evolution of the four dense molecular cores, Sgr B2(N2-N5). We investigate the impact of physical properties and environmental conditions on the formation of specific COMs observed toward Sgr B2(N2-N5) in order to constrain the CRIR and the minimum dust temperature reached during the prestellar phase that best characterize Sgr B2(N).
The article is structured as follows: the observations and methods of analysis are presented in Sect. 2. In Sect. 3 we determine the physical properties of Sgr B2(N2-N5) and build up the physical models used as inputs for the chemical modeling discussed in Sect. 4. The results of our standard models are presented in Sect. 5. In Sect. 6 we discuss the impact of varying the minimum dust temperature and the CRIR on the production of COMs. We compare the model results to the observations to constrain the environmental conditions in Sgr B2(N2-N5). Finally our main results are summarized in Sect. 7.

The EMoCA survey
The EMoCA imaging spectral line survey was conducted with ALMA in its cycles 0 and 1 at high angular resolution (∼1.6 ) with a sensitivity of ∼3 mJy beam −1 per 488 kHz (1.7-1.3 km s −1 ) wide channel. The survey is divided into five spectral setups covering the frequency range between 84. 1 and 114.4 GHz (Table E.1). At these frequencies, the size (half power beam width -HPBW) of the primary beam of the 12 m antennas varies between 69 at 84 GHz and 51 at 114 GHz (Remijan 2015). The phase center is located half way between the two main hot molecular cores Sgr B2(N1) and N2 (α J2000 = 17 h 47 m 19. s 87, δ J2000 = −28 o 22 16. 0). A detailed description of the observations, the calibration, and the imaging procedures is presented in Belloche et al. (2016). Belloche et al. (2016) used the spectra observed toward the main hot cores Sgr B2(N1-N2) to split the line and continuum emission into separated datacubes. In comparison, the spectra observed toward the hot cores Sgr B2(N3-N5) contain significantly less emission lines (that is more emission-free channels, Bonfand et al. 2017). We used the spectra extracted toward their peak positions to split the spectral lines from the continuum emission more accurately, using the CLASS 2 software. In each spectral window of each setup, a first-order baseline is subtracted to the channels that are free of strong line emission to separate the spectral lines from the continuum. The noise level measured in the new continuum-subtracted spectra is listed in  Belloche et al. (2016) in the channel maps of the full field of view for comparison.

Radiative transfer modeling of the line survey
In order to derive the chemical composition of Sgr B2(N3-N5) we performed a radiative transfer modeling of the EMoCA line survey. Given the high densities observed in the hot cores (∼10 7 cm −3 , Bonfand et al. 2017), collisions in the gas phase are frequent enough for the local thermodynamic equilibrium (LTE) approximation to be valid. We used Weeds (Maret et al. 2011), which is part of the CLASS software, to produce LTE synthetic spectra that we compared to the observed continuum-subtracted spectra, after correction for the primary beam attenuation. Weeds solves the radiative transfer equation, taking into account the finite angular resolution of the interferometer, the line opacity, line blending, and the continuum background. The spectroscopic predictions used to model the spectra mainly originate from the CDMS (Cologne Database for Molecular Spectroscopy, Endres et al. 2016) and JPL (Jet Propulsion Laboratory, Pearson et al. 2010) catalog. They are the same as in Belloche et al. (2016Belloche et al. ( , 2017 and Müller et al. (2016), except for ethanol, methanol, and methyl cyanide, for which we used new CDMS predictions. Each molecule is modeled separately, adjusting the following parameters: column density, rotational temperature, angular size of the emitting region (assumed to be Gaussian), velocity offset with respect to the systemic velocity of the source, and linewidth (full width at half maximum -FWHM). For each molecule, a population diagram is plotted to derive the rotational temperature and 2D Gaussians are fit to integrated intensity maps of unblended transitions to measure the size of the emitting region. 1D-Gaussian fits to emission lines are used to derive the linewidth and velocity offset. The column density is adjusted manually until a good fit to the data is obtained. Finally, the contribution from all species is added linearly to build up the complete synthetic spectrum. More details about the methods of analysis can be found in Bonfand et al. (2017).

Physical properties of Sgr B2(N)'s hot cores
In order to model the time-dependent chemical composition of Sgr B2(N)'s sources, we need to determine their physical properties and characterize their time evolution, from the cold ambient cloud phase (that is prior to the free-fall collapse of the cloud) to the dense and warm hot core phase. In this section we derive the physical properties of the sources (H 2 column density, mass, density, size, temperature, and luminosity) based on the observations. In particular, temperature and density will be used to derive the physical profiles used as inputs for the chemical modeling (Sect. 4).

H 2 column densities and masses
We derived H 2 column densities from the dust thermal emission measured in the ALMA continuum emission maps toward each source, using the following equation (Bonfand et al. 2017): with µ H 2 = 2.8 the mean molecular weight per hydrogen molecule (Kauffmann et al. 2008), m H the mass of atomic hydrogen, Ω beam = π 4 ln 2 × HPBW max × HPBW min the solid angle of the synthesized ALMA beam, B ν (T d ) the Planck function at the dust A27, page 3 of 47 A&A 628, A27 (2019) Table 1. Physical properties of four hot cores embedded in Sgr B2(N).
Hot Notes. The uncertainties are given in parentheses. (a) Radius of the mean synthesized beam of the ALMA data (Eq. (6)). (b) Average peak H 2 column density for a mean synthesized beam radius R (Eq. (3)). The uncertainty corresponds to the standard deviation weighted by the error on the ALMA peak flux density, S beam ν , and on the correction factor for the free-free emission (Bonfand et al. 2017). (c) Gas mass of the envelope contained in R (Eq. (5)). (d) Radius of the COM emission region (Sect. 3.2). (e) Excitation temperature derived from population diagrams (Sect. 3.2). ( f ) Dust mass density at the radius r 0 assuming a standard gas-to-dust mass ratio of 100 (Eq. (10)). (g) Density of total hydrogen (n H = n(H)+2n(H 2 )) at r = r 0 (Eq. (15)). (h) Estimated current luminosity of the source (Sect. 3.3). (i) Age of the source estimated from its current luminosity (Sect. 3.4, Fig. 2c). ( j) Mass of the protostar derived from its current luminosity (Sect. 3.4, Fig. F.3). (k) Initial radius of the envelope of the source (Sect. 3.5.2). ( ) In the case of Sgr B2(N3), which is not detected in the EMoCA continuum data at 3 mm, we used the peak flux density reported by Sánchez-Monge et al. (2017), measured at 1.3 mm in an ALMA synthesized beam of 0.4 (that is R = 2003 au), to derive the H 2 column density and mass. temperature T d = 150 K (assuming T d ∼ T g , the gas kinetic temperature, see Sect. 3.2), κ ν the dust absorption coefficient per unit of mass density of gas (in cm 2 g −1 ), and S beam ν the peak flux density measured in Jy beam −1 toward each hot core in the ALMA continuum emission maps, after correction for the primary beam attenuation and the free-free contribution (Bonfand et al. 2017). The dust mass opacity is given by the power law: with κ 0 the reference dust mass absorption coefficient at the frequency ν 0 . In our previous analysis (Bonfand et al. 2017) we derived a dust emissivity exponent β of 1.2 ± 0.1 from the combined analysis of our ALMA data and the dust continuum emission map obtained with the submillimeter array (SMA) at 343 GHz (Qin et al. 2011). According to the simulations of dust coagulation in cold dense cores conducted by Ossenkopf & Henning (1994), a dust emissivity exponent of β = 1.2 suggests an intermediate dust opacity spectrum between the models of dust grains without ice mantles and those with thin ice mantles (see Fig. F.1). Here we adopt the model with grains without ice mantle and for gas densities of 10 6 cm −3 , which corresponds to a dust emissivity exponent of β = 1.3, with κ 0 = 0.0199 cm 2 g −1 (of gas) at ν 0 = 230 GHz (that is λ 0 = 1.3 mm, see Fig. F.1). Table 1 shows the H 2 column densities calculated at T d = 150 K for Sgr B2(N2-N5). In the case of Sgr B2(N3), which is not detected in the ALMA continuum maps at 3 mm (Bonfand et al. 2017), we used the peak flux density measured by Sánchez-Monge et al. (2017) at 242 GHz in ALMA data obtained at 0.4 resolution.
For each hot core, we derived the gas mass of the envelope contained in the 1.6 ALMA mean synthesized beam (or 0.4 for Sgr B2(N3)), from the H 2 column density as follows (Hildebrand 1983): with D = 8.34 kpc the distance to Sgr B2(N) from the Sun (Reid et al. 2014) and R the radius of the mean ALMA synthesized beam: with θ b = 1.6 the mean synthesized beam of the ALMA observations (or 0.4 for Sgr B2(N3)). The resulting masses are given in Table 1.

Source sizes and temperatures
In our previous analysis of Sgr B2(N) we derived the size of the COM emission region, θ s , for the hot cores Sgr B2(N3-N5) (Bonfand et al. 2017). For each source we selected in the observed spectra the transitions that have a high signal-to-noise ratio (S/N), are well reproduced by the LTE model, and are not severely contaminated by other species. We fit 2D Gaussians to their integrated intensity maps to derive the emission size. We derived deconvolved sizes of 1.0 ± 0.3 and 1.0 ± 0.4 for Sgr B2(N4) and N5 respectively, which are slightly smaller than the size of Sgr B2(N2) (∼1.2 , with values ranging between 0.8 and 1.5 ; Belloche et al. 2016Belloche et al. , 2017. Sgr B2(N3) is more compact, with an emission size of 0.4 . For each hot core the radius of the emitting region, r 0 , is computed using Eq. (6) and is given in Table 1. From the COM emission Bonfand et al. (2017) also derived excitation temperatures by plotting population diagrams based on the transitions that are well detected and not severely contaminated by lines from other species. At high densities, where the LTE approximation is valid, the level populations can be described by a single excitation temperature, T rot . We found that T rot varies from ∼145 to 190 K for Sgr B2(N3-N5). The temperatures for Sgr B2(N2) were determined by Belloche et al. (2016Belloche et al. ( , 2017 and Müller et al. (2016). We assume that the excitation temperature is equal to the kinetic temperature and we adopt a gas temperature T 0 of 150 K at the radius r 0 (see Table 1).

Luminosity and temperature profile
We used the radiative transfer code RADMC-3D (version 0.41), to compute dust temperature profiles in the envelopes of Sgr B2(N2-N5), based on their current masses (Sect. 3.1) and temperatures (Sect. 3.2). RADMC-3D performs radiative transfer calculations using the Monte-Carlo method for a given dust distribution and computes the associated dust temperatures. Each source is modeled as a single protostar surrounded by a 1D spherically symmetric envelope with a radius of 10 6 au (that is 4.8 pc), divided into 5 × 10 4 cells. The protostar is defined as a point-source (that is its radius is not taken into account) and we A27, page 4 of 47 M. Bonfand et al.: Complex chemistry in Sgr B2(N)'s hot cores assumed simple black body radiation by setting its effective temperature. For simplicity we adopted a single type of interstellar grains and do not include scattering. We used the dust opacities from Ossenkopf & Henning (1994) for dust without ice mantle and coagulated with gas densities of 10 6 cm −3 as found to best fit Sgr B2(N)'s dust properties (Sect. 3.1).
The dust mass density in the envelope of each source is assumed to follow a power-law distribution: with ρ 0 the dust mass density at the radius r 0 . We assume α = −1.5, holding for a free-falling envelope (Shu 1977). Assuming the following gas-to-dust mass ratio, M d = M g /χ d with χ d = 100, we can write: with M g (R) the gas mass inside of R derived in Sect. 3.1. By replacing ρ d (r) by its expression (Eq. (7)) we obtain: Finally, by integrating Eq. (9) we obtain: The resulting dust mass densities, ρ 0 , are given in Table 1. The dust mass density profiles derived in the envelope of Sgr B2(N2-N5) using Eq. (7) is shown in Fig. F.2. Given the high densities reached in the envelopes, photon packages may get trapped inside optically thick regions, considerably slowing down the simulations. To prevent this effect, we used the Modified Random Walk Method (MRW, Fleck & Canfield 1983) that allows RADMC-3D to predict where the photon package will go next and save computation time.
Finally, for each source the luminosity of its central protostar is adjusted until the calculated dust temperature in the envelope matches the observational constraint, T d (r 0 ) = T 0 . This gives us an estimate of the current luminosity of the sources. Figure 1 shows the dust temperature profiles computed by RADMC-3D in the envelope of the sources along with the total luminosities of their central protostars, which are also listed in Table 1.

Protostellar evolution
In Sect. 3.3 we estimated the current luminosity of Sgr B2(N2-N5). In order to model the time-dependent chemical evolution of these sources, we need to determine the evolution of their protostellar properties, from the earlier stage of star formation, when a young protostar has just formed and started to warm up the surrounding gas. To this aim, we use published results from theoretical simulations of high-mass star formation. Peters et al. (2011) report the results from the first 3D radiation-magnetohydrodynamical (RMHD) simulations of high-mass star formation including ionization feedback. They modeled the free-fall collapse of a magnetized rotating molecular cloud containing 1000 M . The central core, with a gas density of ρ = 1.27 × 10 −20 g cm −3 within a radius of 0.5 pc, is surrounded by an envelope characterized by a density distribution ∝r −1.5 . They assumed an initial temperature of 30 K. After the first 20 kyr of simulation, many sink particles are formed, simulating a group of young stars all contributing to the radiative feedback (see their Fig. 1a, run E). In the initial accretion phase the first sink particle reaches a high accretion rate >10 −3 M yr −1 , which reduces to ∼10 −4 M yr −1 after 0.7 Myr (see their Fig. 1b, run E). At the end of the simulation, the first sink particle is by far the most massive one. Although its accretion rate drops significantly when secondary sink particles form, it continues accreting material until the end of the simulation. Based on the accretion rate evolution of the most massive sink particle, without taking into account the influence of the formation of secondary particles, we approximated the relationṀ (t), characterizing the time-dependent evolution of the accretion rate for a young high-mass protostar. Figure 2a (dashed line) shows the accretion rate as a function of time for a protostar that starts accreting material from its surrounding envelope at t 0 = 0. By integrating the relationṀ (t) we obtain the protostellar mass as a function of time, M (t). Figure 2b (dashed line) shows that after 10 6 yr the protostar reaches a final mass of about 30 M .
In order to estimate the time-dependent evolution of Sgr B2(N2-N5)'s luminosities, L (t), from the relation M (t) (Fig. 2b), we used the mass-luminosity relation derived by Hosokawa & Omukai (2009) for spherically accreting protostars with a constant accretion rate ofṀ = 10 −4 M yr −1 (see Fig. F.3). The resulting total luminosity is plotted as a function of time in Fig. 2c (dashed line). It shows that after 10 6 yr the protostar reaches a final luminosity of ∼10 5 L , which is not high enough to reproduce the luminosities estimated for Sgr B2(N2-N5) (up to 3.9 × 10 5 L , Sect. 3.3). In order to form more massive and thus more luminous protostars, we simply assumed that the accretion rate does not drop below 10 −4 M yr −1 (Fig. 2a, solid line). Using this simple assumption, protostars with masses >10 2 M and luminosities up to ∼2 × 10 6 L are formed after 10 6 yr (Figs. 2b and c, solid lines). We used the relation L (t) (Fig. 2c) to estimate the age, t source , of Sgr B2(N2-N5) based on their current luminosity (see Table 1).
The approach described in this section provides us with a simple model of the physical evolution of young high-mass protostars. Applied to Sgr B2(N2-N5) it allows us to investigate the evolution of their chemical composition using astrochemical models. Our physical model is based on the simple assumption that in all sources, objects with similar final mass are formed, which might not be the case in reality (see discussion in Sect. 6.6). Our treatment is not intended to reflect the complexity of cloud-collapse dynamics but rather to provide a simple framework in which to model the behavior of the cloud chemistry under such conditions. Only a deeper analysis of the whole Sgr B2 molecular cloud would provide a more accurate model for individual sources (see for instance the 3D radiative transfer model of Schmiedeke et al. 2016, which addresses the physical structure of Sgr B2 as a whole, from 45 pc down to 100 au).

The physical model
Modeling the time-dependent chemical evolution of Sgr B2(N2-N5) requires to know their physical structure and evolution over the whole star formation formation process, including the early cold phase preceding the hot core phase. However, the early initial stage of high-mass star formation, prior to the protostar's formation, and its timescale are still poorly known. Based on the fraction of star-forming versus quiescent clumps detected in the APEX Telescope Large Area Survey of the Galaxy (ATLASGAL), Csengeri et al. (2014) estimated an upper limit of ∼7.5±2.5 × 10 4 yr for the quiescent deeply embedded phase (n ∼ 4 × 10 5 cm −3 ) prior to the onset of free-fall collapse. Wilcock et al. (2012) derived statistical lifetimes for starless infrared dark clouds (IRDCs) of ∼2.3 × 10 5 yr from the Herschel Infrared Galactic Plane (HI-GAL) survey. This is a factor three longer than the prestellar phase described by Csengeri et al. (2014) for objects with densities ranging from ∼10 4 to 10 5 cm −3 . More recently, Jeffreson et al. (2018) performed theoretical calculations to estimate the lifetimes of giant molecular clouds based on the time-scales derived for the large-scale dynamical processes that affect the ISM. For clouds located in the CMZ, at galactocentric radii from ∼45 to 120 pc, they estimated median cloud lifetimes of 1.4-3.9 × 10 6 yr.
We divide the early evolutionary phase of high-mass star formation into two distinct stages based on their physical conditions. The first stage traces the evolution of the chemistry for 1 Myr, which corresponds to an intermediate time between the IRDC starless phase derived by Wilcock et al. (2012) and the molecular cloud lifetime determined by Jeffreson et al. (2018). This stage is characterized by low densities and low temperatures in the envelope of the source (see Sect. 3.5.1), which undergoes quasi-static contraction leading to the formation a centrally peaked prestellar core. The second stage starts with the ignition of a young accreting protostar warming up its surrounding envelope. In this stage the envelope undergoes a free-fall collapse (see Sect. 3.5.2), starting from the conditions reached at the end of the first stage. The two stages of our physical model along with the associated timescales are schematized in Fig. 3. The physical properties characterizing each stage are computed in Sects. 3.5.1 and 3.5.2.

Quasi-static contraction
The quasi-static contraction phase prior to the free-fall collapse comprises the evolution of density, visual extinction, and temperature in the envelope of Sgr B2(N2-N5) from t start = −1 × 10 6 yr to t 0 = 0, corresponding to the formation of the stellar embryo and the onset of the free-fall collapse. For simplicity we assumed that the radius does not change and remains constant throughout this phase (that is we look at the evolution of the density, visual extinction, and temperature of a given parcel of gas at a fixed radius, r start , see also Sect. 3.5.2). The quasi-static contraction phase starts at low density, with n H−start = 3 × 10 3 cm −3 where n H is the total hydrogen density (that is n H = n(H) + 2n(H 2 )), and low visual extinction, with A v−start = 2 mag. These are the same initial condition as in Garrod (2013) and subsequent papers. We made the simple assumption that density and extinction at a given radius increase linearly as a function of time in log-log space (see Fig. F.4 and Table E.4).
In order to characterize the behavior of the dust temperature during the cold, low density phase preceding the warming up of the envelope, we adopted the visual-extinction-dependent treatment presented in Garrod & Pauly (2011). Assuming that the rate of cooling of the dust is equal to the rate of radiative heating due to its exposure to a standard interstellar radiation field, ISRF = 1 G 0 in units of the Draine field (Draine 1978), they found that for A v < 10 mag, the dust temperature can be described with a third-order polynomial: (11) Following Garrod & Pauly (2011), at higher visual extinctions (A v > 10 mag) we use the dust temperature profile given by Zucconi et al. (2001 Table 2). The parameters that are common to all sources are highlighted in boldface. The evolution of the physical parameters during each stage is indicated in green with arrows. The evolution of T d during the first stage depends on the adopted minimum dust temperature (see Sect. 3.5.1). For all sources, the quasi-static contraction phase starts with T d−start = 16 K because A v−start = 2 mag, and the dust temperature then decreases as the visual extinction increases. However, dust temperature measurements carried out at infrared wavelengths with the Herschel Space Observatory report somewhat higher temperatures toward the GC region (20-28 K, see Sect. 6.2). Therefore, in order to keep the model physically meaningful and account for the higher dust temperatures expected toward the GC region, we defined an arbitrary minimum temperature, T min , which represents the lowest temperature that is allowed in the chemical simulations (see also Sect. 4.2). In Fig. 4 we plotted the dust temperature as a function of time for Sgr B2(N2-N5) with different T min . The gas kinetic temperature, T g , is held constant throughout the whole quasi-static contraction phase, with T g = T min .
The physical conditions (density, temperature, and visual extinction) achieved at the end of the quasi-static contraction phase are given for each source in Table 2. They represent the initial conditions for the free-fall collapse.

Free-fall collapse
In our physical models, the free-fall collapse phase starts with the ignition of a young high-mass protostar. As the newly formed protostar evolves, it starts heating up the surrounding gas and accretes material from its envelope. In order to characterize the evolution of the physical conditions in the collapsing envelope, we trace the trajectory of a parcel of gas gradually infalling toward the central protostar with the free-fall velocity with M g (r) the mass of gas inside of r. The starting point, r start , of the parcel of gas is chosen such that it reaches r 0 at t = t source (see Sect. 3.4), partaking in the free-fall governed by the enclosed mass. The trajectory, r(t), of a parcel of gas infalling through the envelope of Sgr B2(N2-N5) with the free-fall speed v ff is shown in Fig. F 2.14 (4) 2.19 (6) 1.00 (4) 1.86 (4) 8.27 (7) 1.48 (8) 3.81 (7) 5.91 (7) A v−final (mag) 500 500 500 500 2.63 (5) 4.71 (4) 3.92 (5) 2.82 (5) Notes. Physical properties derived or assumed for each source at three different stages of the high-mass star formation process (see also Fig. 3). X (Y) means X × 10 Y . (a) T d−start = 16 K is valid for the models with T min = 10 and 15 K. For T min = 20 K, T d−start = 20 K.
We derived the gas density along the trajectory of the freefalling parcel of gas from the power-law profile: with n H = 2 n H 2 , the density of total hydrogen, and n 0 the reference density at r 0 . For simplicity we neglected the time variation of n H (r) during the free-fall collapse phase. From Eq. (8), we can write: which, in combination with Eq. (13) yields: The visual extinction along the trajectory of the free-falling parcel of gas is computed from the gas density as follows (Bohlin et al. 1978): with r init the initial radius of the source envelope computed such that n H (r init ) = n start ; thus r init = r 0 n start n 0 − 1 1.5 (see Table 1). During the free-fall collapse phase the visual extinction increases with density, until A v−max = 500 mag. For simplicity we neglected the time variation of A v (r) during the free-fall collapse.
As the protostar's total luminosity rises, the temperatures in the envelope increase, progressively dominating the dust heating over the external radiation field. Based on the evolution of the luminosity over time, L (t) (Sect. 3.1), we computed for each source dust temperature profiles for different luminosities, that is at different stages of star formation (see Fig. F.6). Based on the relation r (t) (Fig. F.5) we derived T d (r) the temperature evolution along the trajectory of a parcel of gas free-falling through the envelopes of Sgr B2(N2-N5) during the collapse phase. The free-fall collapse phase stops when the temperature reaches T max = 400 K. Given the high densities and high visual extinctions reached in the envelope of Sgr B2(N2-N5), gas and dust temperatures are assumed to be well coupled during the free-fall collapse phase such that T g (r) = T d (r).
The gas density, temperature, and visual extinction computed along the trajectory of the parcel of gas infalling through the envelopes of Sgr B2(N2-N5) are plotted in Figs. F.7 and F.8. The physical conditions (radius, density, temperature, and visual extinction) achieved at the end of the free-fall collapse phase are given for each source in Table 2.

A three-phase astrochemical code
In order to interpret the COM abundances observed toward Sgr B2(N)'s hot cores (see Sect. 5.1), we used the astrochemical code MAGICKAL (Model for Astrophysical Gas and Ice Chemical Kinetics and Layering, Garrod 2013) to calculate chemical abundances in the envelopes of Sgr B2(N2-N5). MAGICKAL is a single-point model that allows us to simulate in one stage the time-dependent physico-chemical evolution of a given source, from the early cold phase of the star-formation process to the warm up of the dense cores. The calculated chemical abundances characterize a range of distances from the cold extended region to the warmer innermost part of the protostellar collapsing envelope.
The chemical network is based on that of Belloche et al. (2017). It contains 1333 distinct chemical species of which neutrals can exist either as gas-phase, grain-surface, or icemantle species, while charged species are considered in the gas phase only. The grain surface is defined as the outermost layer of ice that resides on top of the other ice layers or directly onto the grain surface itself. The rest of the ice mantle, composed of the deeper ice layers below the ice surface, is treated as a separate phase. The network comprises 13 374 chemical reactions and processes, which include gas-phase ion-molecule reactions and the related electronic recombination of the resulting protonated ions. Photodissociation and photoionization processes are also included, affecting both the gas-phase and grain-surface chemistry with UV photons from either the ambient radiation field or the CR-induced UV field. Gas-phase and ice-surface chemistry are coupled through the accretion of gas-phase material onto the grains, which may be released back into the gas phase via thermal or non-thermal desorption processes (including chemical reactive desorption, Garrod et al. 2007). Diffusion allows reactions within the mantle as well as exchange of chemical species between surface and mantle via swapping mechanisms. Details on the chemical processes and reactions rates can be found in Garrod et al. (2008), Pauly (2011), andGarrod (2013).
The initial gas-phase material is assumed to be mostly composed of atoms and atomic ions, except for H 2 . Initial abundances are taken from Garrod (2013). As Sgr B2(N)'s hot cores presumably originate from material with the same elemental abundances, the same initial abundances are used for all hot cores.

The chemical models
In this section we present the grid of chemical models (see Table 3) run to investigate the time-dependent chemical evolution of Sgr B2(N2-N5) and explore the influence of environmental conditions on the calculated chemical abundances. The first parameter we vary is the minimum dust temperature, T min , that can be reached in the simulations, in particular during the cold quasi-static contraction phase prior to the onset of the free-fall collapse (see also Sect. 3.5.1). The minimum dust temperature adopted in the chemical simulations is critical because the production of the main ice constituents (H 2 O, CO, CO 2 , and CH 3 OH), essential to the formation of COMs, is highly sensitive to the dust temperature. At low temperatures (∼10 K), mainly hydrogenation reactions are important on the surface of dust grains. As the temperature increases, the surface diffusion rates of heavier species become more important leading to a massive increase in the surface production of simple species by addition of radicals (OH, HCO, NH, NH 2 , CH 3 , CH 2 OH, CH 3 O). However, hydrogen atoms evaporate quickly and hydrogenation reactions become less efficient because the rate of desorption of H is much faster than the rate of reaction between H atoms. In order to explore the impact of the minimum dust temperature on the production of COMs we run "low temperature" models with T min = 10 K, "intermediate temperature" models with T min = 15 K, and finally, to better represent the higher temperatures expected toward the GC, "high temperature" models with T min = 20 K (see Sect. 6.2).
The second aspect explored in the chemical simulations is the impact of the CRIR on the production of COMs. We ran chemical simulations using CRIRs ranging from the standard value ζ H 2 0 = 1.3 × 10 −17 s −1 to 1000 × ζ H 2 0 = 1.3 × 10 −14 s −1 . Given the uncertainties on the column-density dependence of the CRIR in dense regions (see Sect. 6.3), we assumed for simplicity that no attenuation occurs in the envelope of the sources. Due to its location in the CMZ, Sgr B2 is also known to be affected by X-ray emission (see, e.g., Terrier et al. 2018). As the effects of X-rays are often considered to be similar to cosmic ray-related processes (see, e.g., Viti et al. 2003), here the total ionization rate may be seen as a combination of both effects.
For each hot core we ran a grid of 18 models, one for each set of physical parameters (CRIR and T min , see Table 3). The interstellar radiation field is held constant throughout the model grid, assuming a standard value of 1 G 0 . The influence of the ISRF strength on the chemical simulations is discussed in Appendix D.

Observed chemical composition of the hot cores
We used Weeds as described in Sect. 2.3 to derive the chemical composition of Sgr B2(N3-N5) by modeling the continuum subtracted emission spectra observed toward each source. We present the results obtained for 11 COMs that include N-, O-, and S-bearing species, to provide a broad census of the chemical composition of the sources. The parameters of the best-fit LTE models are listed in Table E.2 along with the column densities derived for each molecule. The column densities derived for Sgr B2(N2) are also given for comparison Müller et al. 2016). In the case of methyl cyanide (CH 3 CN), we investigated its isotopologs 13 CH 3 CN and CH 13 3 CN because the vibrational ground state transitions of CH 3 CN are optically thick. We assumed the isotopic ratio 3 CN] = 21 derived by Belloche et al. (2016) to obtain the column density of CH 3 CN.
The column densities derived for the 11 COMs in the four hot cores are plotted in Fig. 5a along with their uncertainties, estimated as described in Appendix A. In order to compute abundances relative to H 2 (Fig. 5b), we used the H 2 column densities derived for the compact region with warm dense gas from which the COM emission arises (see Table E.3). The abundances relative to methanol (CH 3 OH) are plotted in Fig. 5c. The overall chemical composition of Sgr B2(N3) is relatively similar to that of Sgr B2(N5), with at most a factor two of difference between the abundances (relative to CH 3 OH) of all ten COMs. Sgr B2(N4) shows higher abundances of O-and S-bearing species relative to CH 3 OH than Sgr B2(N3) and N5, except for NH 2 CHO, but significantly lower abundances relatives to H 2 . The chemical composition of Sgr B2(N2) differs significantly from that of the other hot cores. For instance, C 2 H 5 CN is about four times more abundant (with respect to CH 3 OH) in Sgr B2(N2) than toward the other hot cores; in addition, the ratio [C 2 H 5 CN]/[C 2 H 3 CN] is three times higher. The abundances relative to CH 3 OH of the O-bearing species CH 3 OCHO and CH 3 OCH 3 are significantly lower toward Sgr B2(N2) compared to the other sources, but similar to Sgr B2(N3) in terms of abundances relative to H 2 . Finally, NH 2 CHO shows an astonishing ratio relative to CH 3 OH about 15 times larger than toward Sgr B2(N3-N5).

Model results
In order to interpret the results derived from the observations (Sect. 5.1) we used the astrochemical code MAGICKAL as described in Sect. 4 to compute time-dependent chemical abundances in the envelopes of Sgr B2(N2-N5). In this section we present the results obtained for the standard models N(2-5)-T15-CR1. All abundances refer to fractional abundances with respect to the total hydrogen density (n H = 2n(H 2 ) + n(H)) unless specified otherwise. The effects on the chemical abundances of variations in the minimum dust temperature reached in the A27, page 9 of 47 A&A 628, A27 (2019) 10 -8 10 -7 10 -6 cold prestellar phase (that is the quasi-static contraction phase) and the CRIR are discussed in Sects. 6.2 and 6.3, respectively. The final gas-phase fractional abundances calculated for the 11 investigated COMs in all models are given in Table B.1.

Building up ice mantles
The 1 Myr cold phase prior to the free-fall collapse allows the gas-phase chemistry to quickly convert free atoms into stable simple molecules. Meanwhile, gas-phase material is accreted onto the interstellar dust grains to form ice mantles. The rate of growth of the ice mantles is not linear as a function of time. It is determined by the net rate of deposition of gas-phase material onto the grains and by the degree of surface coverage of bare grains. For the four standard models, N(2-5)-T15-CR1, an ice thickness of one layer is achieved at a time of 6.4 × 10 3 yr into the simulation (that is t = −0.9936 × 10 6 yr, see Fig. B.1). The rate of mantle deposition is relatively similar for models N2-T15-CR1, N4-T15-CR1, and N5-T15-CR1, over the whole quasi-static contraction phase. After 1 Myr (that is at t 0 = 0), about 110 ice layers have been formed on dust grains for the three models compared to the 161 layers formed in model N3-T15-CR1. This is due to the large densities computed at r start in the envelope of Sgr B2(N3) during the contraction phase (up to 2 × 10 6 cm −3 , see Fig. F.4a). The accretion of material onto the dust grains is thus more efficient in model N3-T15-CR1, leading to rapid deposition of new ice layers. For all models, the net deposition rate continues to increase until the end of the quasi-static contraction phase, indicating that the depletion of gas-phase material has not reached completion yet. Figure 6a shows the chemical composition of the grain mantles at each ice layer based on the surface ice that gets incorporated into the mantle. A timescale is indicated such that a time associated to a given ice layer represents the time at which the ice layer was deposited onto the dust grains. The chemical composition within the first 50 ice layers of the grain mantles is similar for the four standard models. Water is clearly the dominant ice constituent for all models, while CO and CO 2 are also present in significant quantities (see also Fig. 6b). Given the relatively high dust temperatures (T min = 15 K), the slow H-to-H 2 conversion on the grains results in high fractional abundances of atomic hydrogen in the gas phase over the whole quasi-static contraction phase (a few 10 −3 , see Fig. 6c). We find that the surface reaction OH + CH 4 → H 2 O + CH 3 is the main mechanism responsible for the formation of water ice, rather than the reaction of OH with atomic or molecular hydrogen expected at lower dust temperatures (T = 10 K, see, e.g., Garrod & Pauly 2011, see also Sect. 6.2). At the beginning of the simulations, the rapid conversion of atomic carbon to CO in the gas phase causes a sharp decline in the abundance of CH 4 ices, which in turn affects the abundance of water ice (Figs. 6a and c). Due to the high dust temperature (T min = 15 K) the gas-phase CO accreted onto the grains is mobile enough to react with OH radicals to form CO 2 + H, limiting the fraction of OH going toward H 2 O. The abundance of solid-phase CO 2 reaches a peak within the first 10 5 yr of the simulation, representing 60-65% of the water abundance (Fig. 6b). Then, due to the high accretion rate of gas-phase CO onto the dust grains, a rapid switch-over occurs and CO dominates over CO 2 . H 2 CO and CH 3 OH form by the successive hydrogenation of surface CO. Table 4 gives the chemical composition of the ice mantles at the end of the quasi-static contraction phase for the four standard models. The composition of ices observed in the line of sight of the field star Elias 16 (Whittet et al. 1998) and toward the high-mass protostar W33A (Gibb et al. 2000) are also given for comparison. Models N2-T15-CR1, N4-T15-CR1, and N5-T15-CR1 build up ices with similar chemical composition, also comparable (within a factor ∼2) to that of Elias 16. The ice composition in model N3-T15-CR1 differs from that of the other models, in particular with higher abundances of H 2 CO and CH 3 OH (see also Fig. 6a). These differences are due to the higher densities computed during the contraction phase at r start in the more compact inner region of Sgr B2(N3) compared to the other sources. High densities have an impact on the accretion rate of gas-phase species onto the dust grains, such that after 1 Myr simple species start to deplete onto dust grains in N3-T15-CR1 (Fig. 6c).

Production of COMs
In this section we investigate the set of chemical reactions that drive the production of the 11 COMs listed in Table E.2. All activation-energy barriers (when there is one) and binding energies may be found in Garrod (2013) 6. Evolution of the ice-mantle and gas-phase chemistry during the cold quasi-static contraction phase. The results from the standard models N(2-5)-T15-CR1 are displayed from top to bottom. Panel a: calculated fractional ice-mantle composition by ice layer. Panel b: cumulative fractional composition of the ices with respect to water, summed over the grain mantle up to an ice thickness given on the abscissa. Panel c: gas-phase fractional abundances (with respect to total hydrogen) of atomic and simple diatomic species as a function of time with the same time axis as panels a and b to facilitate the comparison.
shows the N-bearing species C 2 H 5 CN, C 2 H 3 CN, CH 3 CN, NH 2 CHO, and CH 3 NCO, as well as CH 3 SH. Panel b shows the O-bearing species CH 3 OH, C 2 H 5 OH, CH 3 OCHO, CH 3 OCH 3 , and CH 3 CHO. In each panel the solid lines indicate the gasphase species while the dotted lines of the same color indicate the solid-phase species, combining both ice-surface and ice-mantle abundances. The calculated fractional abundances are plotted as a function of temperature in the envelope of the source (in log scale) for easy comparison between all models. A timescale is also indicated at the top of each plot. The figure shows that in most cases, the solid-phase abundances of the investigated COMs directly determine the final A27, page 11 of 47 A&A 628, A27 (2019) Notes. Abundances of the main ice constituents, summed over all the ice layers formed at the end of the quasi-static contraction phase. These results are for the standard models and are given in percentage of the water ice value.
References. (a) Observational values from Gibb et al. (2000) and references therein.
gas-phase abundances. This suggests that the gas-phase abundances observed in the warm envelopes of Sgr B2(N2-N5) are dominated by the thermal desorption of the dust-grain ice mantles. In all four models, most COMs are already present on the grains with significant abundances before the warm-up phase (T = 15 K). It is the case for the cyanides C 2 H 5 CN and CH 3 CN, as well as C 2 H 5 OH, CH 3 OH, and CH 3 SH, indicating that their production mostly relies on the early cold chemistry. Other species still form efficiently on the grains up to ∼50 K, in particular the O-bearing species CH 3 CHO, CH 3 OCHO, CH 3 OCH 3 , NH 2 CHO, and CH 3 NCO, suggesting that their production relies on the warm-up stage chemistry. The peak gas-phase fractional abundances of the ten investigated COMs are given in Table 5 along with the temperatures at which they are achieved for each of the standard models. For comparison, we list also the fractional abundances reached at the end of the simulations (that is T = 400 K). Most species reach their peak abundance right after desorption into the gas phase (except for C 2 H 3 CN and CH 3 NCO), then they quickly reach steady fractional abundances. The O-bearing species CH 3 OH, C 2 H 5 OH, NH 2 CHO, CH 3 OCHO, as well as CH 3 CN, C 2 H 5 CN, and CH 3 SH desorb from the grains at high temperatures (120-147 K). The species which desorb at lower temperatures (via thermal or chemical desorption), or which are formed in the gas phase via ion-molecule reactions, are particularly exposed to gas-phase destruction through photodissociation by CR-induced UV photons or reactions with gas-phase ions.
There is little gas-phase formation of methanol; it is mostly formed on the grains via successive hydrogenation of the CO accreted from the gas phase: The high gas-phase fractional abundance of CH 3 OH at low temperature is mostly due to chemical desorption from the dust grains. Model N3-T15-CR1 shows a modest increase in the gas-phase abundance of CH 3 OH around T ∼ 40 K, caused by the electronic recombination of HC(OH)OCH 3 + , product of the reaction between protonated methanol (CH 3 OH 2 + ) and H 2 CO when the latter is released abundantly into the gas phase.
Ethanol, C 2 H 5 OH, is predominantly produced on the grains, mostly via the sequence: (18) followed by later desorption into the gas phase.
Formamide, NH 2 CHO, is mainly produced on the grains via the successive hydrogenation of OCN accreted from the gas phase: The abundance of methyl cyanide, CH 3 CN, mostly derives from the grain-surface hydrogenation reaction: where CH 2 CN is formed via the grain-surface atomic-addition At low temperatures (T < 50 K), CH 3 CN is also produced in the gas phase via the electronic recombination of protonated methyl cyanide (CH 3 CNH + ), product of the ion-molecule reaction between CH 3 + and HCN. Methyl formate, CH 3 OCHO, begins the free-fall collapse phase with modest abundances on the grains. It is formed more efficiently at later times (T ∼ 20-30 K), mainly via the grain-surface radical-radical addition reaction: followed by later desorption into the gas phase. Model N3-T15-CR1 shows a modest increase in the gas-phase abundance of CH 3 OCHO around 40 K, when (H 2 CO) is abundantly released into the gas phase and reacts with CH 3 OH 2 + to form HC(OH)OCH 3 + . Gas-phase CH 3 OCHO is then produced via the electronic recombination of HC(OH)OCH 3 + . Methyl mercaptan, CH 3 SH, is predominantly formed on the grains via the successive hydrogenation reactions: During the cold quasi-static contraction phase, CH 3 SH forms via the successive hydrogenation of surface CS accreted from the gas phase. Above ∼20 K, H 2 CS is accreted onto the grains directly from the gas phase, where it is formed via the electronic recombination of H 3 CS + , deriving from the ion-molecule reaction S + + CH 4 .
Vinyl cyanide, C 2 H 3 CN, is formed on dust grains through the successive hydrogenation of HC 3 N accreted from the gas phase: Then, as hydrogenation continues on the grain surfaces, C 2 H 3 CN is quickly converted to C 2 H 5 CN as follows: Ethyl cyanide, C 2 H 5 CN, can thus maintain a steady abundance on the grains, dominating C 2 H 3 CN at all times. Only a small fraction of solid-phase C 2 H 3 CN contributes to the gasphase abundance as it is quickly destroyed via a ion-molecule reaction right after desorption. At late times C 2 H 3 CN may form in the gas phase via the reaction: Due to their low binding energies to water ices, CH 3 NCO (3575 K, Belloche et al. 2017), CH 3 OCH 3 (3675 K, Garrod 2013), and CH 3 CHO (2275 K, Garrod 2013) desorb at low temperatures (that is around 60-80 K). Methyl isocyanate, CH 3 NCO, starts the free-fall collapse phase with modest A27, page 12 of 47 M. Bonfand et al.: Complex chemistry in Sgr B2(N)'s hot cores Fig. 7. Calculated fractional abundances for 11 COMs detected toward Sgr B2(N2-N5). Panels a: N-bearing species C 2 H 5 CN, C 2 H 3 CN, CH 3 CN, NH 2 CHO, and CH 3 NCO, as well as CH 3 SH. Panels b: O-bearing species CH 3 OH, C 2 H 5 OH, CH 3 OCHO, CH 3 OCH 3 , and CH 3 CHO. The results for the standard models, N(2-5)-T15-CR1, are shown from top to bottom, as a function of the temperature in the envelope of the sources during the free-fall collapse phase. In each panel the timescale shown at the top is derived from the time-dependent evolution of the temperature in the envelopes of the sources (see Sect. 3.5.2). In each panel, the solid lines show the fractional abundances (with respect to total hydrogen) in the gas phase while the dotted lines show the abundances of the same species on the grains (ice surface+mantle).
abundances on the grains. It is more efficiently produced around 20 K through the radical-addition reaction: At low temperatures (T < 50 K), CH 3 NCO may be destroyed on the grains via hydrogenation reactions to form CH 3 NHCO. In addition, CH 3 NCO is quickly destroyed via gas-phase ion-molecule reactions right after desorption, leading to very low gas-phase fractional abundances at the end of the simulations (<10 −14 with respect to total hydrogen, for all standard models but N3-T15-CR1). In model N3-T15-CR1, CH 3 NCO is formed much more abundantly on the grains, mainly because of the high abundance of surface CH 3 accreted from the gas phase. Therefore it cannot be destroyed completely once it is released into the gas phase, explaining the higher fractional abundances of gas-phase CH 3 NCO observed in this model.    Notes. X(Y) means X × 10 Y . (a) Gas-phase fractional abundances (with respect to total hydrogen) of the standard models at the end of the simulations. (b) Gas-phase peak fractional abundances (with respect to total hydrogen) reached during the free-fall collapse phase. (c) Temperature at which the peak fractional abundance is achieved.
Dimethyl ether, CH 3 OCH 3 , is mostly formed on the surface of dust grains up to T ∼ 20-30 K, through the radical-addition reaction: In model N3-T15-CR1, the gas-phase fractional abundance of CH 3 OCH 3 increases when CH 3 OH desorbs from the dust grains around 130 K. This is due to the electronic recombination reaction CH 3 OCH 4 + + e − → CH 3 OCH 3 + H, where CH 3 OCH 4 + is a product of the reaction between CH 3 OH and CH 3 OH 2 + . Finally, acetaldehyde, CH 3 CHO, is formed on the grains up to T ∼ 30 K, via the surface reactions: Surface CH 3 CHO may be destroyed either by reacting with CH 3 to form CH 4 + CH 3 CO, or with NH 2 to form NH 3 + CH 3 CO. In particular in models N3-T15-CR1 and N4-T15-CR1, only a small fraction of solid-phase CH 3 CHO contributes to the gas-phase fractional abundance. At lower temperature (T ∼ 24 K), surface CH 3 CHO is released into the gas phase via chemical desorption from the reaction CH 2 CHO + H. Once in the gas phase, CH 3 CHO becomes the dominant reaction partner for gas-phase ions, damping ionic abundances and thus limiting the destruction of other gas-phase species. In models N3-T15-CR1 and N4-T15-CR1, CH 3 CHO is formed directly in the gas phase via the reaction C 2 H 5 + O → CH 3 CHO + H, when C 2 H 5 desorbs from dust grains around 90 K.
The gas-phase fractional abundances of the 11 investigated COMs obtained at the end of the chemical simulations (that is at T = 400 K) are given for the standard models in Table 5 and plotted in Fig. 8. For each of the 11 COMs the abundances calculated in models N2-T15-CR1 and N5-T15-CR1 are relatively similar. Model N4-T15-CR1 shows significantly lower abundances of CH 3 CHO, CH 3 OCHO, and NH 2 CHO than the other models. Model N3-T15-CR1 significantly differs from the other models. It shows higher abundances of cyanides, and has an astonishing [C 2 H 5 CN]/[C 2 H 3 CN] ratio >18 times lower than the other models. The final fractional abundance of CH 3 OH is 2-5 times higher than in the other three models. Finally, N3-T15-CR1 shows an astonishing high abundance of CH 3 NCO, more than five orders of magnitude higher than for the other models.

Physical properties of Sgr B2(N)'s hot cores
In Sect. 3.1 we computed H 2 column densities from the dust continuum emission measured toward Sgr B2(N2-N5) in the ALMA continuum emission maps. The results reported in Table 1 show that at a given radius, Sgr B2(N2) has the highest density, lying a factor ∼3 above N3 and N4 which have similar densities. The density computed for Sgr B2(N5) lies in between that of the other sources. Given the uncertainties on the calculated densities (Table 1), the differences between the hot cores may be marginal. The results discussed above slightly differ from our previous analysis of the same ALMA data (Bonfand et al. 2017), where for the H 2 column density calculations we assumed a dust mass opacity coefficient (κ 100 GHz = 5.9 × 10 −3 cm 2 g −1 ) ∼1.1 times smaller than the value we use at the same frequency (100 GHz their Table 3) from the integrated flux measured over a source size of 0.59 (see Table 1 of Sánchez-Monge et al. 2017). We note that they used a mean molecular weight of µ = 2.33, which corresponds to that of the mean free particle, so that their calculated value has to be divided by ∼1.2 to obtain the H 2 column density. This gives N H 2 = 3.7 × 10 24 cm −2 , a factor ∼4 higher than the value we derived in this paper. The difference is due to the different assumptions made to compute the H 2 column density. They assumed a lower dust temperature of 100 K and used a dust mass opacity coefficient (κ 230 GHz = 0.011 cm 2 g −1 ) ∼1.8 times smaller than our value at the same frequency.
In Sects. 3.1 and 3.3 we used H 2 column densities to compute masses and densities, which are important inputs for the chemical modeling. One must keep in mind that the uncertainties in H 2 column densities may have an impact on the chemical model results.

Influence of the minimum dust temperature on the chemistry
As mentioned in Sect. 3.5.1, dust temperature measurements carried out at infrared wavelengths with the Herschel Space Observatory report somewhat higher dust temperatures toward the GC region than other typical regions forming high-mass stars. For instance, Longmore et al. (2012) derived dust temperatures ranging from 19 to 27 K from the center to the edge of the galactic center cloud G0.253+0.016 (also known as "the Brick"). Guzmán et al. (2015) reported measurements of T d = 20-28 K toward Sgr B2, somewhat higher than the 17 K measured for dense cores embedded in the NGC 6334 molecular cloud lying at a distance of 1.3 kpc from the Sun (Chibueze et al. 2014).
Recently Belloche et al. (2016) suggested that the low deuteration fractions of COMs observed toward Sgr B2(N) could be explained by either an elemental abundance of deuterium toward the GC region lower than in the solar neighborhood, or higher temperatures during the prestellar phase that would have decreased the efficiency of the deuteration process. Deuteration levels of methanol similar to that of Sgr B2(N) measured toward the aforementioned high-mass star-forming region NGC 6334I located at a distance of ∼7 kpc from the GC favor the second hypothesis (Bøgelund et al. 2018). Gas kinetic temperatures ranging from 50 K to more than 100 K were found on the basis of imaging of H 2 CO transitions with the APEX 12 m telescope (Ao et al. 2013;Ginsburg et al. 2016). Ott et al. (2014) showed that the kinetic temperature distribution of molecular clumps within the region between the supermassive black hole Sgr A * and Sgr B2 peaks at ∼38 K based on a survey of NH 3 conducted with the Australia Telescope Compact Array (ATCA). In comparison, lower temperatures have been derived from NH 3 lines toward several massive clumps outside the GC region (T g = 10-30 K, Giannetti et al. 2013).
In this section we explore the impact of varying arbitrarily the minimum dust temperature, T min (see also Sect. 3.5.1), on the formation of COMs, to account for the high dust temperatures observed toward the GC. Figure 9 shows the chemical composition of the ice mantles built up during the 1-Myr cold contraction phase for the "low temperature" model N2-T10-CR1 (panel a), and the "high temperature" model N2-T20-CR1 (panel b). Significant differences can be seen in the overall chemical composition of the ice mantles compared to the standard model (that is N2-T15-CR1, Fig. 6a, top panel), indicating that the minimum dust temperature adopted in the chemical simulations is of critical importance. The ice mantle formed in the low-temperature model is enriched in CO and depleted in CO 2 compared to the standard model. This is due to the less efficient CO-to-CO 2 conversion when the dust temperature drops below ∼14 K in the low-temperature model. At low temperature, H is converted more efficiently into H 2 due to its lower thermal desorption rate, thus decreasing the abundance of gas phase atomic hydrogen compared to the standard model. Surface OH mostly reacts with H 2 to form water ice rather than with CO to form CO 2 . More CO is thus available to form H 2 CO and CH 3 OH, enhancing their abundances in the outer ice layers compared to the standard model. , assuming a standard CRIR and different T min : 10 K in blue, 15 K in purple, and 20 K in red (see Table 3).
The high-temperature model produces thinner ice mantles than N2-T10/15-CR1, with fewer than 70 ice layers formed on dust grains after 1 Myr. The ice mantles are enriched in carbon dioxide as most surface CO and OH are destroyed to form CO 2 , decreasing the amount of H 2 O, H 2 CO, and CH 3 OH formed on the grains. The atomic carbon accreted to the grains mostly reacts with surface O 2 instead of reacting with atomic hydrogen to form CH, which affects the abundance of CH 4 . Furthermore, at high temperature reactions involving atomic hydrogen are less efficient because the H coverage of the grains is reduced due to rapid desorption, cutting down the formation of complex species.
The minimum dust temperature adopted in the chemical simulations also has an impact on the abundances of complex species at later times during the free-fall collapse phase, as the investigated COMs have been shown to be mainly produced on the grain surfaces (see Sect. 5.2.2). Figure 10 compares the final gas-phase fractional abundances of the low-, intermediate-, and high-temperature models (see also Figs. B.2-B.5). In most cases, a higher T min results in lower fractional abundances. In particular, the O-bearing species CH 3 OH, CH 3 NCO, CH 3 CHO, CH 3 OCH 3 , and CH 3 OCHO are very sensitive to T min with, for a given source, fractional abundances lower by more than two orders of magnitude in the high-temperature models compared to the low-temperature model. The fractional abundance of CH 3 NCO varies by more than four orders of magnitude between the low-and high-temperature models for N2, N3, and N5. The fractional abundance of CH 3 SH is barely sensitive to changes in T min , with at most a factor 3.2 of difference between the low-, intermediate-, and high-temperature models for all sources. This is due to its production mechanism, which relies on gasphase H 2 CS which sticks to the dust grains at T ∼ 20 K (see Sect. 5.2.2).
A more accurate representation of the dust temperatures in the envelope of Sgr B2(N2-N5), in particular in the outer parts exposed to external UV photons, requires to take into account the dust heating via the interstellar radiation field. The influence of the ISRF strength on the chemical model results is explored in Appendix D. It shows that assuming a radiation field stronger than the standard ISRF value does not improve the agreement between calculated and observed abundances with respect to CH 3 OH for Sgr B2(N2-N5) (see discussion in Sect. 6.4).  Finally, a more realistic treatment of the grain-surface chemistry would require considering a grain-size distribution rather than single dust-grain radius (0.1 µm) adopted in our simulations. Pauly & Garrod (2016) investigated the effects of a broad grain-size distribution on the ice-surface chemistry by implementing a distribution of five distinct initial grain sizes to their dark cloud models, where the radii of each grain population increase as ice mantles form on the grains. These authors found that when using a uniform dust-grain temperature across all grain sizes the surface chemistry does not significantly differ from that of models assuming a single dust-grain radius. However, for chemical models with a more explicit treatment of the grain temperatures, in which T d varies also with the grain-mantle growth of each individual grain-size population, Pauly & Garrod (2016) found that the temperature difference plays a significant role in determining the total amount of ice formed and its composition for each grain-size population due to the temperature-sensitive nature of grain-surface reactions.

Influence of the cosmic-ray ionization rate on the chemistry
In this section we explore the impact of varying ζ H 2 on the formation of COMs. Figure 11 for instance shows the composition of the ice mantles built up during the cold phase prior to the free-fall collapse for Sgr B2(N2), with ζ H 2 = 10 × ζ H 2 0 (N2-T15-CR10), 100 × ζ H 2 0 (N2-T15-CR100) and 1000 × ζ H 2 0 (N2-T15-CR1000). Panel a shows that the chemical composition of the ices exposed to a CRIR enhanced by a factor of 10 are relatively similar to that obtained for the standard model N2-T15-CR1 (Fig. 7a, top panel). Bigger differences can be seen in the composition of the ice mantles for higher ionization rates. For instance in the models N2-T15-CR100 and N2-T15-CR1000, thicker ice mantles are formed on the grains (up to 153 ice layers for N2-T15-CR100 after 1 Myr, Fig. 11b), than in the low-CRIR models (ζ H 2 ≤ 10 × ζ H 2 0 ). At higher ionization rates, water is still the main ice constituent, but the grain mantles are strongly depleted in CO 2 compared to the low-CRIR models, because most surface CO 2 is photodissociated via CR-induced UV photons. The abundance of CO in the ice is also diminished because gas-phase CO reacts with He + before accretion onto dust grains. Gas-phase molecular hydrogen is rapidly dissociated, leading to higher abundances of atomic hydrogen in the gas phase compared to the standard models. Surface reactions involving atomic hydrogen are thus more efficient, as seen in model N2-T15-CR100 which produces ice mantles enriched in CH 4 and NH 3 compared to the other models. NH 3 maintains a fairly stable abundance through the ice layers since the surface NH 3 photodissociated to NH 2 is quickly hydrogenated to form NH 3 again. Finally, the ice mantles in the high CRIR models (ζ H 2 ≥ 100 × ζ H 2 0 ) are significantly depleted in H 2 CO and CH 3 OH, which are efficiently destroyed by reacting with surface atomic hydrogen.
Besides their impact on the production of ice-mantle constitutents, cosmic rays also play an important role in gas-phase chemical processes, producing the ions that react with gas-phase species leading to their destruction. The results of models N2-T15, N3-T15, N4-T15, and N5-T15 with different CRIRs are plotted in Figs. B.6, B.7, B.8, and B.9, respectively. If we focus for instance on the results obtained for model N2-T15 with different CRIR, Fig. B.6 shows that higher gas-phase fractional abundances of C 2 H 5 OH and CH 3 NCO are obtained at the end of the simulation for a CRIR enhanced by a factor ten compared to the standard value. The final gas-phase fractional abundance of CH 3 CHO obtained for N2-T15-CR10 is one order of magnitude lower than that of N2-T15-CR1. This is because of reactions with atomic hydrogen destroying CH 3 CHO on the surface of dust grains. Only a small fraction of the CH 3 CHO formed on the grains thus desorbs around 60 K and contributes to its final gasphase abundance. The CH 3 CHO present in the gas phase at lower temperatures (that is before thermal desorption) mainly derives from chemical desorption, but once in the gas phase it is quickly destroyed as it reacts with ions.
At higher CRIR, in the model N2-T15-CR100, C 2 H 5 OH reaches a final gas-phase fractional abundance of about two orders of magnitude higher than in the standard model. CH 3 NCO is formed approximately three orders of magnitude more abundantly than in the standard model although it is strongly destroyed by reacting with atomic hydrogen. Its final gas-phase abundance derives from the CR-induced photodissociation of CH 3 NHCO which desorbs from the grains around T ∼ 90 K. Finally in the model with the highest ionization rate, N2-T15-CR1000, all COMs are affected by cosmic rays (Fig. B.6, bottom panel). In particular, the abundances of COMs with little to no gas-phase formation routes decrease right after desorption from the grain surfaces, leading to low final gas-phase fractional abundances (Fig. 12). This is the case for C 2 H 5 OH, CH 3 OH, CH 3 OCHO, CH 3 OCH 3 , NH 2 CHO, C 2 H 5 CN, and CH 3 SH, which are all destroyed in the gas phase as they react with one of the most abundant ions, H 3 O + , produced from water after it desorbs into the gas phase:  fractional abundance a N2-T15    Table 3. Species with fractional abundances lower than 10 −15 are not visible in the histograms.
The gas-phase abundance of CH 3 OCH 3 decreases by more than three orders of magnitude right after it desorbs from the dust grains. At later times, CH 3 OCH 3 is formed in the gas phase via dissociative recombination of CH 3 OCH 4 + , a product of the ion-molecule reaction CH 3 OH 2 + + CH 3 OH, when CH 3 OH desorbs from dust grains, around 130 K. The main contribution to the final gas-phase abundance of C 2 H 3 CN and CH 3 CN also comes from gas-phase reactions. C 2 H 3 CN is mostly formed via the gas-phase reaction CN + C 2 H 4 → C 2 H 3 CN + H, when C 2 H 4 desorbs from dust grains around 70 K, and the dissociative recombination of C 2 H 6 CN + . The later reaction is particularly efficient since C 2 H 6 CN + is produced from the ionmolecule reaction H 3 O + + C 2 H 5 CN. The gas-phase abundance of C 2 H 3 CN thus increases as C 2 H 5 CN is destroyed. CH 3 CN reaches its peak gas-phase abundance at early times, around 57 K, after which it is strongly destroyed by ion-molecule reactions. The main contribution to CH 3 CN's final gas-phase abundance comes from gas-phase recombination of CH 3 CNH + , product of the ion-molecule reaction CH 3 + + HCN. Figure 12 shows that the final gas-phase fractional abundances of the 11 investigated COMs are very sensitive to the CRIR for Sgr B2(N2-N5). For instance, the abundances of cyanides rise as the CRIR increases up to 100 × ζ H 2 0 . For higher ionization rates (ζ H 2 ≥ 500 × ζ H 2 0 ), the final gasphase abundance of C 2 H 5 CN decreases significantly, leading to [C 2 H 5 CN]/[C 2 H 3 CN] < 1 in all models (see also Fig. 16). The final fractional abundances calculated by models N3-T15 and N4-T15 are strongly diminished at high CRIR (1000 × ζ H 2 0 ) compared to the low ionization-rate models, except for CH 3 CHO. Model N2-T15-CR1000 shows significantly higher fractional abundances of C 2 H 5 OH, CH 3 CHO, and CH 3 NCO compared A27, page 18 of 47 to the standard model N2-T15-CR1. The impact of varying the CRIR on the calculated fractional abundances depends on the source, suggesting that the chemical differences discussed here may be exacerbated by the different physical properties of the cores.

Main trends and uncertainties
The gas-phase fractional abundances obtained at T = 150 K (which corresponds to T 0 , see Sect. 3.2) for the 11 investigated COMs in all models are given in Table B.2. The results of all chemical models are also plotted in Fig. 13 compared to the observed abundances with respect to CH 3 OH. It shows that in most cases the chemical models tend to underestimate the observed abundances. In particular, this is the case for CH 3 OCHO, for which the calculated abundances are lower than the observed ones toward all four sources. The abundance of CH 3 NCO with respect to CH 3 OH is also underestimated by all models, by at least two orders of magnitude compared to the observed abundances. In Sect. 5.2.2 we have seen that solidphase CH 3 NCO desorbs from the grains around 70 K, however, the measured rotational temperatures toward Sgr B2(N2-N5) are higher than 140 K (Sect. E.2). This suggests that some important chemical reactions involved in the production of CH 3 NCO may be missing in the chemical network (see, e.g., Quénard et al. 2018), or that the binding energy adopted for CH 3 NCO in our chemical models (3575 K, Belloche et al. 2017) is underestimated (see Belloche et al. 2019, for models computed with a higher binding energy). For this reason we ignored CH 3 NCO in the rest of the discussion.
In Fig. 13 no error bars are plotted on the model results although it has been shown that calculated chemical abundances may be affected by rate coefficient uncertainties (see, e.g., Wakelam et al. 2005). Given the large number of parameters in the chemical code, a full analysis taking into account both observational and theoretical uncertainties (see, e.g., Wakelam et al. 2006) would be too difficult and is not attempted here. Besides, additional uncertainties affecting the model results derive from the lack of experimental and theoretical data on grain-surface processes. For instance, Iqbal et al. (2018) found that the abundances of both gas-phase and surface species calculated under typical dark cloud conditions are affected by large errors due to the uncertainties in the diffusion energy of surface species, in particular for the molecules predominantly formed at the grain surface like CH 3 OH. Enrique-Romero et al. (2016) presented quantum chemical computations showing that grainsurface radical-radical addition reactions do not necessarily lead to the formation of complex species because these radicals are trapped by the water-ice molecules with an orientation that favors other two-product reactions. This suggests that gas-phase reactions may play a more important role in the formation of COMs than currently assumed in chemical models. Recently Skouteris et al. (2018) computed rate coefficients for two gas-phase ion-molecule reactions proposed as main formation route for ceCH 3 OCH 3 : CH 3 OH + CH 3 OH 2 + → (CH 3 ) 2 OH + + H 2 O followed by the reaction of (CH 3 ) 2 OH + with NH 3 to form CH 3 OCH 3 + NH 4 + . The latter reaction, not included in our chemical network, could help in improving the agreement between the calculated and observed [CH 3 OCH 3 ]/[CH 3 OH] ratios, which differ by at least one order of magnitude in our models for N3-N5 (see Fig. 13).

Constraining the physical parameters T min and ζ H 2
We used the method presented by Garrod et al. (2007) to evaluate the level of confidence of each model with respect to the observations in order to determine which set of physical parameters (T min and CRIR, that is ζ H 2 ) best characterizes Sgr B2(N)'s hot cores (see Appendix C for details about the method). We draw our attention first to the overall abundances of COMs with respect to H 2 , and then examine in more detail the influence of T min and CRIR on the chemical composition itself, that is the COM abundances relative to CH 3 OH or CH 3 CN. Figure 14 shows the confidence levels for the ten investigated COMs (excluding CH 3 NCO) relative to H 2 , for Sgr B2(N2-N5) taken individually and for all sources taken together. Panels a, b, and c show the confidence levels for all molecules, only the O-bearing molecules, and only the cyanides, respectively. The confidence level matrices differ from one source to the other, maybe reflecting the limits of our strong assumption that all sources share the same accretion history. Given that all sources should have been exposed to the same cosmic-ray flux and probably shared a similar thermal history during the prestellar phase, we focus our attention on the matrices corresponding to all sources taken together. The cyanides appear to be much more sensitive to T min and CRIR than the O-bearing species. The cyanides constrain the CRIR value to be about 50 times the standard one, and the minimum temperature to be below 20 K. The best confidence levels for the O-bearing species are also obtained for these values, but the constraints are less sharp.
Given that a model with T min = 20 K still gives acceptable results for the O-bearing species, we ran additional simulations with T min = 25 and 28 K for N2, with ζ H 2 = 50 × ζ H 2 0 (see Fig. F.9a). Ignoring the cyanides, we obtained confidence levels of about 10 and 5% for the models N2-T25-CR50 and N2-T28-CR50, respectively. These values are much lower than the 23% obtained for N2-T20-CR50. This results from the high dust temperatures (T min ≥ 25 K) preventing gas-phase atoms and simple molecular species such as CO to stick to the grains, damping the formation of ice mantles and thus the production of complex species. We expect the models of N3-N5 to behave in the same way and conclude that the COM abundances relative to H 2 in N2-N5 certainly exclude T min > 20 K.
We now focus on the COM chemical composition itself and show in Fig. 15 the matrices of confidence levels computed for nine COMs relative to CH 3 OH (panel a), the O-bearing molecules relative to CH 3 OH (panel b), and the cyanides relative to CH 3 CN (panel c). This refinement reveals that the O-bearing COM chemical composition relative to CH 3 OH is sensitive to the CRIR value, with the best-fit value being 50 times the standard value, in agreement with the constraints obtained from the cyanide abundances relative to H 2 (see Fig. 14). A minimum temperature below 20 K seems to be favored as well. The cyanide chemical composition with respect to CH 3 CN is less sensitive to the two investigated parameters (bottom panel of Fig. 15c). However, a more detailed investigation reveals that the abundance ratio of C 2 H 5 CN to C 2 H 3 CN, two species that are chemically linked (see Sect. 5.2.2), is sensitive to the CRIR. Figure 16 shows that this ratio decreases as the CRIR increases in the four sources Sgr B2(N2-N5). The observed ratios are best reproduced by the models with 50 × ζ H 2 0 ≤ ζ H 2 ≤100 × ζ H 2 0 , due to the increase in the gas-phase abundance of C 2 H 3 CN at late times (that is high temperatures) in these models. This behavior is consistent with what we observe, in particular toward Sgr B2(N2) where the rotational temperature derived for C 2 H 3 CN is higher than that  of C 2 H 5 CN. The flat gas-phase abundance profile of C 2 H 3 CN at T > 40 K in the standard models with lower CRIR (Fig. 7) would imply that it traces more extended regions than C 2 H 5 CN, which would be inconsistent with what is observed in Sgr B2(N2). Our analysis above shows that the chemical models with a CRIR enhanced by a factor 50 compared to the standard value best reproduce the observations for the four hot cores taken together. However, one must keep in mind that for simplicity our models assume that cosmic rays are not attenuated in the envelope of the sources, although the large visual extinctions of Sgr B2(N2-N5) (see Figs. F.4b and F.7b) may actually lead to attenuation of the cosmic-ray flux. Neufeld & Wolfire (2017) showed, based on H 2 and H 3 + column densities measured toward diffuse clouds in the galactic disk, that the CRIR decreases with increasing H 2 column density with a best-fit dependence scaling as N(H 2 ) −a with a = 0.92 ± 0.32 for N(H 2 ) ∼ 10 20 -10 22 cm −2 . Panel a: matrice of confidence levels of the models with respect to the observed abundances of ten COMs relative to H 2 . The first four panels in this column show, from top to bottom, the results for Sgr B2(N2), N3, N4, and N5, respectively. Bottom panel: average matrice for all four hot cores taken together. Panel b: same as (a) but only for the O-bearing species CH 3 OH, CH 3 OCHO, CH 3 OCH 3 , CH 3 CHO, C 2 H 5 OH, and NH 2 CHO. Panel c: same as (a) but only for the cyanides CH 3 CN, C 2 H 5 CN, and C 2 H 3 CN.
This implies that the cosmic-ray flux is attenuated by a factor ∼70-600 from the diffuse regions with N(H 2 ) = 10 21 cm −2 to the dense gas with N(H 2 ) = 10 23 -10 24 cm −2 . The attenuation of the cosmic-ray flux might have a significant impact on the chemical model results. For instance, Rimmer et al. (2012) showed that a column-density-dependent CRIR improves the agreement between chemical model predictions and observations for the Horsehead nebula, compared to a standard model with a constant CRIR. For a more sophisticated treatment of the cosmic-ray driven chemistry in hot core models including cosmic-ray attenuation see, Willis et al. (in prep.). Moreover, our models do not take into account direct cosmicray bombardment of dust-grain ice mantles although laboratory experiments showed that it can trigger a rich chemistry that may lead to the formation of more complex species even at very low temperatures (Hudson & Moore 2001;Rothard et al. 2017). Recently, Shingledecker et al. (2018) used chemical models to investigate the effect of direct cosmic-ray collisions with dust grains on the solid-phase chemistry, including the formation of suprathermal species in the ices from collision with energetic particles. For instance, for a standard CRIR (10 −17 s −1 ), they found that the abundance of CH 3 OCHO, which is systematically underproduced in all our models (Fig. 13), is significantly enhanced in the gas phase, as well as on the grains compared to the standard chemical model (that is without cosmic-ray driven reactions).
As reported above, our models firmly exclude minimum dust temperatures of 25 K or higher during the prestellar phase of N2-N5. Such high temperatures would prevent the production of COMs at the level observed in these sources. Our models show that a value of 15 K, higher than the low values used in previous studies (≤10 K, see, e.g., Müller et al. 2016;Belloche et al. 2017;Garrod et al. 2017), still leads to an efficient production of COMs roughly consistent with the observations, while for 20 K signs of disagreement emerge between model and observations, especially for the cyanide abundances relative to H 2 and the O-bearing COM chemical composition relative to CH 3 OH. A more accurate treatment of the dust temperatures in the envelope of Sgr B2(N2-N5) for different values of the ISRF does not modify the results much compared to our simple parametrization with an artificially-set minimum temperature (see Table D.1). A value of 15 K is close to the lowest dust temperatures measured in the GC region (19-20 K, see Sect. 6.2). It could even be that models with T min =17-18 K produce COMs in amounts that are Confidence level (%) Fig. 15. Panel a: matrice of confidence levels of the models with respect to the observed abundances of nine COMs relative to CH 3 OH. The first four panels in this column show, from top to bottom, the results for Sgr B2(N2), N3, N4, and N5, respectively. The bottom panel shows the average matrice for all four hot cores taken together. Panel b: same as (a) but for the abundances relative to CH 3 OH of the O-bearing species CH 3 OCHO, CH 3 OCH 3 , CH 3 CHO, C 2 H 5 OH, and NH 2 CHO only. Panel c: same as (a) but for the abundances relative to CH 3 CN of the cyanides C 2 H 5 CN and C 2 H 3 CN only.
still consistent with the abundances measured in Sgr B2(N)'s hot cores. In addition, the coldest dust grains in the GC region may be masked by warmer outer layers and thus have been missed by Herschel because of its limited angular resolution. Therefore, it is likely that the constraint we derived on T min from the COM abundances is fully consistent with the thermal properties of the GC region.

Chemical differentiation between N-and O-bearing complex molecules
Chemical differentiation between O-and N-bearing species has been observed toward several high-mass star-forming cores (Blake et al. 1987;Wright et al. 1992;Wyrowski et al. 1999;Beuther et al. 2005;Allen et al. 2017). Recently, Csengeri et al. (2018) investigated the young high-mass star-forming region G328.2551-0.5321 where CH 3 OH and other O-bearing species show two distinct bright emission peaks spatially offset from the peak of the cyanides (see also Csengeri et al. 2019). The authors suggested that in this source, CH 3 OH traces the accretion shocks resulting from the interaction between the collapsing envelope and the accretion disk, while the cyanide emission may trace the more compact hot region radiatively heated by the protostar. In this scenario, the O-bearing COM emission detected in the envelope of G328.2551-0.5321 arises from accretion shocks rather than from the region radiatively heated by the protostar itself. If such a mechanism is at play in Sgr B2(N) at scales smaller than the beam of the EMoCA survey, then infering the luminosity of the protostars from the rotational temperature measured at the radius of the COM emission as we did in Sect. 3.3 would be wrong. The thermal history derived for Sgr B2(N2-N5) (see Sect. 3.5.2) would then be incorrect, and this could have a significant impact on the results of our chemical modeling.
Observations at higher angular resolution that are attained in the ReMoCA survey ) are needed to investigate whether a chemical differentiation between N-and O-bearing species occurs in Sgr B2(N)'s hot cores. so, we used the relation L (t) obtained from RMHD simulations of high-mass star formation (Peters et al. 2011, see also Sect. 3.4), combined with the mass-luminosity relation given by Hosokawa & Omukai (2009). The estimated t source values are unlikely to be accurate because our analysis assumes that the evolution of the accretion rate,Ṁ (t), is the same for all sources, although the density measured at a given radius toward Sgr B2(N2) is a factor 3.2, 3.5, and 2.0 higher than that of N3, N4, and N5, respectively (Fig. F.7b). In addition, the simulations of Peters et al. (2011) were set up for a molecular cloud with a central core gas mass density of ρ g = 1.27 × 10 −20 g cm −3 within a radius of 0.5 pc, which is about 10-40 times lower than the values we derive for Sgr B2(N2-N5) at the same radius ( Fig. F.7b).
Using the mass-luminosity relation provided by Hosokawa & Omukai (2009) for a constant accretion rate of 10 −4 M yr −1 , the final mass of ∼30 M obtained using the model from Peters et al. (2011) leads to a final luminosity of ∼10 5 L , significantly lower than the values we derive for Sgr B2(N2-N5) (Sect. 3.3). Our simplified treatment, in which we assumed that the accretion rate does not drop below 10 −4 M yr −1 to form more massive and more luminous objects (Sect. 3.4), results in significant differences in the time evolution of Sgr B2(N3)'s physical parameters (n H , A v , and T ) compared to the other three sources. Given its lower luminosity, N3 is apparently younger than the other sources (see Table 1). In the framework of our model, this results in a faster increase of the density and temperature along the trajectory of the modeled parcel of gas during the free-fall collapse phase of N3. Figure F.9 shows that our best models (that is with ζ H 2 = 50 × ζ H 2 0 ) better reproduce the chemical abundances (with respect to H 2 ) of C 2 H 5 OH and CH 3 OH in Sgr B2(N3) compared to the other sources. However, the modeled abundances of CH 3 OCHO and CH 3 OCH 3 are 2-3 orders of magnitude lower than the abundances measured toward N3. This might suggest that the relationsṀ (t), M (t), and L (t) we derived in Sect. 3.4 do not represent accurately the physical evolution of Sgr B2(N2-N5). Furthermore, in our previous analysis of Sgr B2(N)'s hot cores (Bonfand et al. 2017), we proposed a tentative evolutionary sequence in which N4 appears younger than N3, based on the molecular outflow signature detected toward N3, while N4 does not show any evidence of outflow yet. This is in contradiction with the present analysis in which we find that N3 is the youngest of the four investigated hot cores. This is another indication that the short evolution timescale we obtain here for N3 may not represent well the physical evolution of the source. However, the origin of some wing emission seen in the spectrum of N4 is not understood (see discussion in Sect. 3.7 in Bonfand et al. 2017) and a search for a bipolar outflow in N4 should be carried out with more robust outflow tracers, such as CO, whose J = 1-0 transition is not covered by the EMoCA survey.
A27, page 23 of 47 A&A 628, A27 (2019) We also point out that the ages we obtain for Sgr B2(N2), N4, and N5 in Sect. 3.4 are about four times longer than the hot core lifetime we estimated in our previous analysis of Sgr B2(N)'s hot core population (∼6 × 10 4 yr, Bonfand et al. 2017) on the basis of a very simple statistical argument involving the number of hot cores and UCHII regions detected in the same area and assuming a lifetime of ∼10 5 yr for the UCHII phase (Peters et al. 2010). This might suggest that the accretion rate,Ṁ (t), of these sources is/was actually higher than what we have assumed here.
In order to better represent the extreme physical conditions (densities/masses, luminosities) of Sgr B2(N)'s hot cores and give an accurate estimate of their age, RMHD simulations treating specifically their different physical properties will be needed. In addition a more self-consistent treatment of the accretion rate, mass, and luminosity evolution (that is, taking into account the time-dependent evolution of the accretion rate to derive the mass-luminosity relation) is also required.

Conclusions
We analyzed the results of the 3 mm imaging line survey EMoCA conducted with ALMA in its cycles 0 and 1 to characterize the hot core population embedded in Sgr B2(N). We derived the chemical composition of the three hot cores Sgr B2(N3-N5) focusing on 11 COMs. We investigated the timedependent evolution of the protostellar properties and we derived the evolution of density, temperature, and visual extinction in the envelopes of Sgr B2(N2-N5). We used the astrochemical code MAGICKAL to simulate the time-dependent chemical evolution of the sources based on their physical evolution. We investigated the impact of the cosmic rays and the minimum dust temperature reached during the cold prestellar phase on the chemical abundances calculated by the models. We compared the results to the observations to constrain the physical parameters which best represent the environmental conditions characterizing Sgr B2(N) and the GC region. Our main results are the following: 1. Sgr B2(N3) and N5 share a similar chemical composition relative to methanol. The chemical composition of Sgr B2(N2) differs significantly from the other hot cores Sgr B2(N3-N5). 2. We derived gas densities of 1.4 × 10 7 cm −3 toward Sgr B2(N2) at a radius of 6010 au (∼0.03pc), while the values derived for N3, N4, and N5 at the same radius are a factor 3.2, 3.5, and 2.0 lower, respectively. 3. On the basis of the rotational temperature derived for the COMs and using the results of RMHD simulations of highmass star formation, the current luminosities of the hot cores are estimated to be 2.6 × 10 5 L , 4.5 × 10 4 L , 3.9 × 10 5 L , 2.8 × 10 5 L for Sgr B2(N2), N3, N4, and N5, respectively. These results hold only if our assumption that the COM emission arises from the region radiatively heated by the protostar is valid. 4. The production of the cyanides C 2 H 5 CN and CH 3 CN as well as CH 3 OH and C 2 H 5 OH mostly relies on the early cold grain-surface chemistry (T d 15 K), while the O-bearing species CH 3 CHO, CH 3 OCHO, CH 3 OCH 3 are predominantly formed during the warm-up phase of the protostellar evolution. 5. COMs form efficiently on grains during the prestellar phase of Sgr B2(N2-N5) with minimum dust temperatures as high as 15 K. Minimum dust temperatures 25 K are firmly excluded by our chemical models. 6. The cosmic-ray ionization rate adopted in the chemical simulations is critical, in particular during the warm-up stage where cosmic rays may rapidly destroy COMs once they are released into the gas phase. The best match between chemical models and observations is obtained for ζ H 2 = 7 × 10 −16 s −1 , that is a CRIR enhanced by a factor 50 compared to the solar neighborhood value. This is somewhat lower than extreme values expected toward the diffuse medium in the GC region (1 -11 × 10 −14 s −1 ). This difference may reflect the attenuation of cosmic rays in denser gas. The chemical composition calculated by our models for Sgr B2(N)'s hot cores strongly depends on their history since the initial stages in the prestellar phase. Our strong assumption about the accretion history being the same for all sources implies a much shorter timescale for Sgr B2(N3) which has a significant impact on its calculated chemical composition and results in a worse agreement with the observations. Furthermore, our model in which the three other hot cores, N2, N4, and N5 are by far older than N3 is inconsistent with the tentative evolutionary sequence suggested in our previous analysis of Sgr B2(N)'s hot cores based on their association with class II methanol masers, outflows, and UCHII regions (Bonfand et al. 2017). This emphasizes the need for RMHD simulations tailored to each hot core in order to have a more realistic description of their history and improve the reliability of the chemical model results.
The column densities estimated in Sect. 5.1 for Sgr B2(N3-N5) (see also Table. E.2) are based on simple assumptions, such as a single excitation temperature to characterize all transitions from a given molecule under the LTE approximation. In addition, because the hot core emission toward Sgr B2(N3-N5) is resolved in our data for only a few species, especially toward Sgr B2(N3) for which transitions from only two species could be used to derive the source size (Bonfand et al. 2017), we assume for each of these three sources that all molecules are emitted from the same region characterized by a single size. This may not be true as it is shown in Sect. 5.2.2 that all investigated COMs do not desorb from the dust grains at the same temperature in the chemical models, suggesting that they may trace different regions.
In order to estimate the uncertainties on the column densities we varied for each species detected toward Sgr B2(N3-N5) the emission size and rotational temperature within the uncertainties given by the 2D Gaussian fits to the integrated intensity maps and population diagrams, respectively. Comparing the column densities derived with our best-fit synthetic spectra (Table E.2) with those obtained varying temperature and source size, we find that column densities vary at most by a factor four. In the case of Sgr B2(N2) we simply assume uncertainties of 30% on the molecular column densities of all investigated species.      In order to determine which set of physical parameters (T min and CRIR) best characterizes Sgr B2(N)'s hot cores (see Sect. 6.4.2), we evaluate the success of each model in reproducing the observations using the method presented by Garrod et al. (2007). For each model, we compute a level of confidence, κ i , in the agreement between calculated and observed abundances (relative to H 2 , CH 3 OH, or CH 3 CN). For each calculated abundance ratio, the level of confidence is given by:

Appendix B: Results of the chemical simulations
which computes the logarithmic distance of disagreement between the calculated abundance ratio (R model,i ) and the observed one (R obs,i ). We define the standard deviation σ = 1 such that one standard deviation corresponds to one order of magnitude lower/higher than the observed abundance ratio. erfc is the complementary error function (erfc = 1−erf) such that κ i ranges between 0 and 1. For instance, a calculated abundance ratio lying one order of magnitude lower/higher than the observed one has a confidence level κ i = 31.7% (and κ i = 4.6% for two orders of magnitude etc.). For each model we take the mean of the confidence levels obtained for the molecules in the sample to define the overall confidence level of the model with respect to the observations. The results are given in Tables C.1-C.6.   Notes. Same as C.1 but for the abundances of nine COMs relative to CH 3 OH. (a) Because NH 2 CHO is not detected toward N4 and we estimated only an upper limit to its molecular column density, a level of confidence of unity is attributed if the calculated ratio [NH 2 CHO]/[CH 3 OH] is lower than the upper limit, otherwise it is treated normally.

Appendix D: Influence of the external radiation field on the chemistry
Investigations carried out on the extreme environmental conditions characterizing the GC region suggest that the ambient radiation field may be enhanced by a factor ∼500-1000 in the inner Galaxy (Lis et al. 2001;Clark et al. 2013), compared to the standard ISRF value (1 G 0 ) used in our simulations. A stronger ISRF may not have a big impact on the warm chemistry during the free-fall collapse phase, because external UV photons are rapidly absorbed as the visual extinction increases. However, it may be more critical during the preceding low-density phase, where dust grains are not efficiently shielded against external radiation. Recently, Hocuk (2017) where χ UV represents the strength of the UV field in units of the standard Draine field. In order to explore the impact of a stronger radiation field on the production of COMs, we run three additional chemical models, N2-UV100, N2-UV500, and N2-UV1000, with χ UV = 100 G 0 , 500 G 0 , and 1000 G 0 , respectively. No minimum threshold is fixed for the dust temperatures, which go as low as ∼14, 17, and 19 K in models N2-UV100, N2-UV500, and N2-UV1000, respectively. The gas temperature, T g , is held constant at 15 K during the whole quasi-static contraction phase. We adopt ζ H 2 = 6.5 × 10 −16 s −1 as found to be the CRIR that best reproduces the observations toward Sgr B2(N2-N5) (see Sect. 6.4.2). Figure D.1 compares the calculated gas-phase fractional abundances of ten COMs relative to CH 3 OH for the models with strong UV fields to the observed abundances. The results of model N2-CR50-UV1, with a standard radiation field and a minimum dust temperature T min = 15 K, are also shown for comparison. It shows that stronger UV fields increase the discrepancy between the calculated and observed abundances of cyanides with respect to CH 3 OH, while the abundances of NH 2 CHO, CH 3 CHO, and C 2 H 5 OH relative to CH 3 OH are only slightly affected by the UV field strength. Table D.1 gives the mean confidence level of the models with the observations for different χ UV . It shows that the highest confidence level for the abundances of nine COMs (excluding CH 3 NCO) with respect to CH 3 OH is obtained with the model assuming a standard ISRF (1 G 0 ). The confidence levels drop significantly for models with higher χ UV . The opposite behavior is obtained when computing confidence levels with only the cyanides C 2 H 3 CN and C 2 H 5 CN with respect to CH 3 CN (that is the agreement with observations increases with the UV field strength). Finally, all models have similar confidence levels when computed based only on the O-bearing species CH 3 OCHO, CH 3 OCH 3 , CH 3 CHO, C 2 H 5 OH, and NH 2 CHO. This does not allow us to conclude strictly on which interstellar radiation field better characterizes Sgr B2(N). This is consistent with the results we obtain by simply parametrizing the dust temperature profile with an artificially-set minimum temperature (see Sect. 6.4.2).
As pointed out by Hocuk (2017), the dust temperature expression of Garrod & Pauly (2011) results in higher temperatures at low visual extinctions than the values computed using the      Notes. (a) Number of lines detected above 3σ (Table E. Fig. 5b to plot chemical abundances relative to H 2 . The uncertainties are calculated based on the uncertainties on the average H 2 column densities given in Table 1 and on the source sizes. ( * ) In the case of N3 we use the H 2 column density given in Table 1, measured in the 1 mm ALMA continuum map at 0.4 resolution ( Table 1).

Index of power-law function Source n H (t)
A Bonfand et al. (2017) this paper no ice -10 6 cm −3 no ice -10 8 cm −3 thin ice -10 6 cm −3 thin ice -10 8 cm −3 3 ×10 5 3 ×10 4 3 ×10 3 3 ×10 2 ν (GHz) Fig. F.1. Dust mass opacity spectra for a standard gas-to-dust mass ratio of 100. The dashed lines show the results of Ossenkopf & Henning (1994)'s simulations after 10 5 yr of dust coagulation with gas densities of 10 6 cm −3 and 10 8 cm −3 and for dust grains without ice mantle and with thin ice mantle. The blue solid line shows the κ ν (of gas) used to compute H 2 column densities in Bonfand et al. (2017). The black solid line shows the κ ν (of gas) used in this paper, computed with Eq. (4), and based on Ossenkopf & Henning (1994)'s model for grains without ice mantle and for gas densities of 10 6 cm −3 (magenta dashed line).  (Shu 1977). For each source the observational constraint, ρ(r 0 ) = ρ 0 is plotted with errorbars (1σ). For each source the parcel of gas is held at a constant radius, r start (see Table 2), during the whole quasi-static contraction phase. During the free-fall collapse phase, the parcel of gas gradually falls toward the central protostar with the free-fall speed v ff (Eq. (12)).