Ice Age: Chemodynamical Modeling of Cha-MMS1 to Predict New Solid-phase Species for Detection with JWST

Chemical models and experiments indicate that interstellar dust grains and their ice mantles play an important role in the production of complex organic molecules ( COMs ) . To date, the most complex solid-phase molecule detected with certainty in the interstellar medium is methanol, but the James Webb Space Telescope ( JWST ) may be able to identify still larger organic species. In this study, we use a coupled chemodynamical model to predict new candidate species for JWST detection toward the young star-forming core Cha-MMS1, combining the gas – grain chemical kinetic code MAGICKAL with a 1D radiative hydrodynamics simulation using Athena ++ . With this model, the relative abundances of the main ice constituents with respect to water toward the core center match well with typical observational values, providing a ﬁ rm basis to explore the ice chemistry. Six oxygen-bearing COMs ( ethanol, dimethyl ether, acetaldehyde, methyl formate, methoxy methanol, and acetic acid ) , as well as formic acid, show abundances as high as, or exceeding, 0.01% with respect to water ice. Based on the modeled ice composition, the infrared spectrum is synthesized to diagnose the detectability of the new ice species. The contribution of COMs to IR absorption bands is minor compared to the main ice constituents, and the identi ﬁ cation of COM ice toward the core center of Cha-MMS1 with the JWST NIRCAM / Wide Field Slitless Spectroscopy ( 2.4 – 5.0 μ m ) may be unlikely. However, MIRI observations ( 5 – 28 μ m ) toward COM-rich environments where solid-phase COM abundances exceed 1% with respect to the column density of water ice might reveal the distinctive ice features of COMs.


INTRODUCTION
Over the past decades, the rich chemistry of the interstellar medium (ISM) has been unveiled largely through dedicated ground-based microwave and sub-mm observations of molecular rotational transitions.Complex organic molecules (COMs), typically defined as carbon-bearing molecules containing 6 or more atoms, are found to be ubiquitous in the gas throughout various evolutionary stages of star formation, from extremely cold prestellar cores to chemically-rich sources such as hot corinos (e.g.Blake et al. 1987;Arce et al. 2008;Bottinelli et al. 2010;Öberg et al. 2010;Bacmann et al. 2012;Jørgensen et al. 2012;Fayolle et al. 2015;Jørgensen et al. 2016).
Chemistry taking place both in the gas and on the surfaces of interstellar dust grains has been proposed to explain such a high degree of chemical complexity in star-forming regions.Earlier theories relied more or less exclusively on gas-phase mechanisms such as ion-molecule reactions to drive the chemical complexity (e.g.Hamberg et al. 2010;Vasyunin & Herbst 2013), following the sublimation of simple ice species from the grains.However, the timescales available for gas-phase COM production have been argued to be too short (e.g.Aikawa et al. 2008), encouraging the deeper exploration of an alternative explanation -grain-surface/ice chemistry.In this scenario, COMs would already have been formed on the grains or within their ice mantles prior to ice sublimation and subsequent detection in the gasphase.Although there has recently been renewed interest in gas-phase processes, including radiative association and neutral-neutral reactions (Balucani et al. 2015;Shannon et al. 2013), grain-surface chemistry is nevertheless considered to be -at minimum -a substantial contributor to COM production in star-forming regions.The recent modeling study of hot core chemistry by Garrod et al. (2022) indicates that in those hot sources, at least, COM production on grains outweighs that occurring in the gas phase.
Until recently, the grain-surface/ice chemistry was believed to be related to the warming of the dust, with intermediate temperatures (T dust 30 K) being sufficient to allow heavy radicals on the grains to become thermally mobile and thus reactive (Garrod & Herbst 2006).However, recent detection of COMs toward prestellar cores shows that large species can be synthesized during the very cold, early phase of star formation in which thermal diffusion of radicals is highly inefficient.Laboratory studies, and now gas-grain chemical models, show that the formation of COMs appears to be possible in such cold and quiescent environments through nondiffusive reactions, and/or radiolysis within ice mantles (Fedoseev et al. 2015;Chuang et al. 2016;Shingledecker et al. 2018;Garrod 2019;Jin & Garrod 2020).This new picture of grain-surface/ice chemistry during star formation suggests that much of the gas-phase COM material observed under hot conditions may have been formed during the very earliest stages of ice-mantle growth (Garrod et al. 2022), although some processing likely occurs at later times and higher temperatures.A recent laboratory study shows that species as complex as glycine (NH 2 CH 2 COOH) may also be formed under cold, prestellar conditions, given an appropriate ice composition (Ioppolo et al. 2020).
Although not directly detected, the presence of highly complex molecules on grain surfaces has frequently been inferred, based on a variety of evidence.For example, many laboratory studies indicated that UV and energetic radiation of ices results in the formation of COMs (e.g.Gerakines et al. 1996;Bernstein et al. 2002).The chemical richness revealed in hot cores by gas-phase observations may indirectly indicate the formation of COMs in interstellar ices; in this scenario, an embedded protostar heats up the surrounding envelope, the ensuing thermal release of simple solid-phase molecules would be accompanied by the release of complex species that are also present in the ice (see the review of Jørgensen, Belloche & Garrod 2020).
However, while the inventory of gas-phase COMs has continued to expand, our understanding on ice-mantle composition in interstellar clouds and star-forming sources has been limited to the most abundant, simple ice species (H 2 O, CO, CO 2 , NH 3 , CH 4 , H 2 CO, and CH 3 OH) in absorption from 3 to 15 µm.Several absorption features have been attributed to frozen COMs, although the identification with particular molecular species remains a hotly debated topic, in contrast to gas phase studies (reviewed by Boogert et al. 2015).The observational identification of larger ice species is challenging; in addition to the lack of infrared (IR) spectroscopic data from the laboratory, the spectral resolution and sensitivity of past astronomical instruments were insufficient to identify the weak and complicated IR bands of COMs, and to distinguish them from the deep, broad features presented by the abundant simple molecules.Furthermore, the detection of various ice species (see Figure 7 of Boogert et al. 2015) requires a line of sight of sufficiently high visual extinction, thus narrowing the number of sources in which COMs could plausibly be detected.
With the launch of the James Webb Space Telescope (JWST), our limited observational knowledge of solid-phase COMs is expected to expand rapidly, particularly through the approved JWST Early Release Science (ERS) program Ice Age.This program will trace the evolution of pristine and complex ice chemistry during different evolutionary stages of low-mass star formation, from prestellar cores to protoplanetary disks, through observations of the Ced 110 region of the Chamaeleon molecular cloud I (Cha I) complex, located at a distance of 192 pc (Dzib et al. 2018).Many background stars exist within this field, providing more than 140 lines of sight, which would allow the determination of spatial profiles for ice abundances.In particular, four of these lines of sight probe the ice within 10,000 AU of the Class 0 protostar "Chamaeleon-MMS1" (Cha-MMS1) which is a main target of this study.Due to the benefits of observing ice formation in this field, one of the main science goals of the Ice Age project is to identify new COM species within ice mantles.If successful, such identifications would provide a direct observational link between the well-known gas-phase presence of COMs in warm/hot sources and their expected production on dust grains at earlier times.
In preparation for the Ice Age project, we have modeled the time-dependent chemical evolution and structure of one of its target sources, Cha-MMS1, to guide the observational study of solid state COMs towards this object.This very young star-forming core is an interesting source from an astrophysical perspective.Since Reipurth et al. (1996) first identified this source, it has drawn attention as a promising candidate for the first hydrostatic core (FHSC) stage of evolution due to its weak mid-IR emission at 24µm and 70µm (Belloche et al. 2006).The FHSC is defined as an adiabatic core in a quasi-hydrostatic equilibrium state where the gravity is balanced by the thermal pressure gradient (Larson 1969).Because of the short lifetime of this stage, only one FHSC candidate, Barnard 1b-N (Gerin et al. 2015), has so far been successfully identified as such.High-resolution observations toward Cha-MMS1 have recently become available with ALMA, but the evolutionary stage of this source remains highly debated.Busch et al. (2020) detected a CO outflow towards Cha-MMS1 with high angular resolution observations.The deprojected outflow velocities exceed the range of the expected values from a FHSC outflow, implying that Cha-MMS1 has already gone through the FHSC phase; those authors suggest a dynamical age of 200-3000 yr post-FHSC.However, another recent paper, Maureira et al. (2020), reported weak or missing emission of HCO + , HCN, and CH 3 OH toward the Cha-MMS1 envelope.ALMA observations of CO, CS, H 2 CO, and CH 3 OH towards an outflow associated with Cha-MMS1 also show that the outflow is weaker than the ones produced by Class 0 objects (Allen et al. 2020).However, regardless of whether Cha-MMS1 is in the FHSC stage or is a young Class 0 source, we would expect the state of the ice mantles, especially in the more extended regions, to be largely unaffected by this change.
In this work, the chemical structure of Cha-MMS1 is modeled by using an astrochemical model combined with a radiative hydrodynamics simulation (RHD) of FHSC with a simple 1-D approach.This paper is structured as follows: The chemical model and its implementation within a dynamical simulation are described in § 2. The results of the models are discussed in § 3. Conclusions are summarized in § 4.

CHEMO-DYNAMICAL MODEL
The main modeling focus of this study involves the coupling of a three-phase chemical kinetics model to the outputs of a radiation hydrodynamics simulation (RHD) of the formation of the first hydrostatic core, in order to investigate the chemistry of Cha-MMS1.This chemo-dynamical modeling method enables us to explore the chemical responses to the dynamical change, achieving a more realistic chemical picture of the source.Thus far, several models have employed this approach to understand the dynamical / chemical evolution of FHSC.For example, Furuya et al. (2012) performed three-dimensional RHD simulation coupled with astrochemical models and showed that large organic molecules are associated with the FHSC.Hincelin et al. (2013Hincelin et al. ( , 2016) ) developed the same kind of simulation with the consideration of the effect of magnetic field.Although these models are powerful, the chemical treatments of these studies were focused on the gas and grain-surface chemistry only, using a so-called two-phase approach.However, full consideration of the chemistry including deeper icy mantle layers is needed particularly for the investigation of astronomical ices, because of the preservation of the chemical composition of earlier surface layers during the evolution of the cloud (Garrod & Pauly 2011).Here, we use a fully-coupled gas-phase / grain-surface / icy grain-mantle chemical model that also includes non-diffusive surface/ice-mantle chemistry (Jin & Garrod 2020).This allows us to build up simulated profiles of gasand solid-phase species through the core that may ultimately be compared with new data from the JWST Ice Age project.The approach used here is broadly similar to that of Barger et al. (2021), who applied a 1-D RHD method to the evolution of high-mass hot cores (note also that the gas-grain chemical model used by those authors included only diffusive grain-surface/ice-mantle chemistry, whereas the present model also employs non-diffusive chemistry).
Firstly, a one-dimensional collapse from a flat density profile (n H ≃ 6.3 × 10 4 cm −3 ) to a self-supported FHSC is simulated using the hydrodynamic code Athena++ (Stone et al. 2020).The model setup is very similar to that presented by Masunaga et al. (1998).Lagrangian particles are placed throughout the initial density profile, which then evolve over time to track the motions of mass-conserved fluid parcels, or shells, as the core collapses.The time-dependent physical conditions and radial positions of each particle (which we refer to as trajectories) determined by the dynamical code are then fed into the chemical model, MAGICKAL, to build up a corresponding chemical picture of the core.The chemical model is run independently for each trajectory, after the RHD simulation is complete; thus, there is no direct chemical interaction between the different trajectories, and the chemistry does not affect the dynamics in any way.
A total of 47 trajectories were selected to describe the time-dependent 1-D physical structure, for each of which the chemical evolution was calculated.Each trajectory has its own set of time-dependent density, temperature, visual extinction and radius values.Due to the 1-D nature of the simulations, the radial ordering of the trajectories is maintained.The visual extinction for each trajectory throughout its physical evolution (which is needed for the chemical rates and to determine the dust temperature) is calculated in post-processing, by evaluating the one-dimensional integral of the total hydrogen density profile to determine the total H column density, then applying the relation of Bohlin et al. (1978), i.e. (1) Here, N H is the column density from the instantaneous position of a particular trajectory to the outer edge of the core.To this value for A V , an assumed background visual extinction provided by material surrounding the core, A V,bac , is added, to account for external material that is not explicitly considered in the RHD calculations.A value of A V,bac = 2.5 mag is adopted, which corresponds to the observed visual extinction threshold of 5 mag for starless cores in the Cha I cloud (Belloche et al. 2011), within which Cha MMS-1 resides.
The initial conditions of the RHD simulation are intended to represent a gravitationally-bound dense condensation within the dark cloud (slightly more massive than the Jeans mass).The chemistry of this initial state should be expected already to have evolved somewhat; thus, a prior stage of collapse is also modeled, allowing the chemistry to evolve from that appropriate to a translucent cloud to that of the dense core from which the RHD model begins.This pre-RHD chemical model (referred to hereafter as the "pre-model") is run individually for each RHD-model trajectory, and consists of a simple (0-D) freefall collapse treatment, in which the gas density evolves from 3 × 10 3 cm −3 to the (uniform) starting density of the dense core of ∼ 6.3 × 10 4 cm −3 over a period of approximately 1 Myr (see e.g.Nejad et al. 1990).The visual extinction at each position is scaled during the pre-model evolution according to A V ∝ n 2/3 H , in keeping with past implementations of the collapse treatment (e.g.Garrod et al. 2022), such that the visual extinction profile that is present at the initiation of the RHD simulations is attained at the end of the pre-model.The additional background extinction is added to all calculated values, as described above.
Gas temperatures in the pre-model and main chemical model are set to 10 K or to the RHD model-calculated value, whichever is greater.The dust temperature used in the chemical models is chosen as the greater of the RHD-model value and the visual extinction-dependent value obtained from the relation provided by Garrod & Pauly (2011).Both dust temperature treatments assume a minimum value of 8 K.The details of the RHD simulation are described in § 2.2.A more detailed description of the chemistry used in the model is provided in § 2.1.

Chemical model
The chemical simulations conducted for each of the dynamical trajectories use the astrochemical gas-grain kinetics code MAGICKAL.The coupled gas-phase, grain-surface, and icy mantle chemistry are simulated by solving a set of rate equations (Garrod & Pauly 2011;Garrod 2013), supplemented by the modified rate-equation treatment presented by Garrod (2008).
Following Jin & Garrod (2020) and Garrod et al. (2022), the chemical model includes both the typical diffusive grain-surface/ice chemistry and a range of nondiffusive mechanisms by which grain/ice surface or bulk-ice species may meet and thence react.These mechanisms are important to treat accurately the production of COMs at low temperatures.The model includes the following nondiffusive mechanisms: three-body (3-B) reactions, in which a preceding surface/bulk reaction is followed by the immediate reaction of the product with some pre-existing nearby species; photodissociation-induced (PDI) reactions, in which a radical produced through photodissociation of a precursor instantly meets a reaction partner in its immediate vicinity; and the Eley-Rideal process, in which gas-phase species adsorb directly onto a grain-surface reaction partner.The "PDI2" treatment of Garrod et al. (2022) is employed, in which (non-hydrogenic) photodissociation products in the bulk ice that do not find an immediate reaction partner are allowed to recombine with each other.The three-body excited-formation (3-BEF) mechanism proposed by Jin & Garrod (2020), in which the chemical energy released by a preceding reaction allows the energized product to overcome the activation energy barrier to the follow-on reaction, is also used, adopting the generalized treatment of Garrod et al. (2022).All chemical reactions included in the surface/bulk network are able to occur through both diffusive and nondiffusive mechanisms (excluding the 3-BEF mechanism, which occurs only in cases where the initiating reaction produces sufficient energy to overcome the barrier to the subsequent reaction).A discussion of the effects of nondiffusive processes in the model is provided in § 3.8.Importantly, the present model (following Garrod et al. 2022) eliminates bulk diffusion in the ice for all species but atomic and molecular hydrogen.Although this means that thermally-driven, diffusive bulk-ice reactions between radicals do not occur, those radicals may still react through the nondiffusive mechanisms described above.H atoms may also react diffusively with those radicals, as well as with more stable species such as CO via activation energy barrier-mediated reactions.Those diffusive hydrogen atoms are typically produced via photodissociation of molecules.Thus, both the direct photodissociation-driven "PDI" chemistry and the indirect reactions between atomic H and other species are ultimately the result of UV-induced photodissociation of bulk-ice molecules.
The photodissociation rates of grain-surface/bulk species is assumed to be a factor 3 lower than the corresponding gas-phase processes, although the other rate coefficients and the photoproducts are the same.This is based on the modeling study of Kalvāns (2018) that theoretically determines an average, general ratio between the photodissociation coefficients for molecules in ice and gas.
Initial elemental abundances used in the chemical model are shown in Table 1.These values are partially based on those used by Wakelam & Herbst (2008) and Jin & Garrod (2020).However, here we adopt the "low-metal" value (Graedel et al. 1982) for nitrogen.In the dark clouds that should provide the starting point for low-mass star formation, metals (all elements other than hydrogen and helium) are depleted in the gas compared to diffuse clouds, because heavy elements have already been incorporated into solids as the density evolves.Jenkins (2009) suggests a relationship to describe the degree of depletion of each metal element with density, and Wakelam & Herbst (2008) define the initial elemental abundances for chemical models by applying Jenkins' relation to a typical density of dense clouds (2×10 4 cm −3 ).Previously, we did not use the low-metal value only for nitrogen, because nitrogen is the element that does not show distinctive depletion with density evolution (Jenkins 2009).Here, we assume some depletion for nitrogen in dense clouds, consistent with other elements; as argued by Jenkins, observational bias could be possibly involved in measuring the nitrogen depletion.The newly adopted initial abundance for nitrogen is around 3 times lower than the values found in the ζ Oph diffuse cloud (Hincelin et al. 2011).In the models, this change has the effect of reducing the ultimate solid-phase ammonia abundance to a value in keeping with observations of interstellar ice features, ensuring that COMs whose production derives from ammonia also will not be overproduced for this reason.The underlying chemical network used in these models is based on that presented by Garrod et al. (2022), with a number of additions and changes made to the grain-surface/bulk-ice portion of the network.Particular adjustments are made to the reaction schemes involving the formation of formic acid (HCOOH, FA) and ethanol (C 2 H 5 OH, EtOH), which are summarized in Table 2.
Two reactions on grain surfaces mainly contribute to the formation of FA: Reaction (2) has two equally weighted branches resulting in stable products, one of which is formic acid.Reaction (3), which is mediated by a modest activation energy barrier, typically leads to carbon dioxide production, but may also produce the radical COOH, which can easily be transformed into FA through hydrogenation by mobile atomic H.The branching ratio (BR) to form COOH versus CO 2 + H is not well constrained.Garrod et al. (2022) proposed a 1% conversion, although the resulting gas-phase formic acid in those hot core models, derived from grain-surface production, was overproduced by a factor 2 compared with well-observed sources Sgr B2(N2) and IRAS 16293-2422 (high-and low-mass hot core/corino sources, respectively).We here adjust the ratio for COOH production downward by a factor 2, to 0.5%.The production of ethanol on the grains in this model is strong, and can occur early in the evolution.As a result, we paid particular attention to ensuring that the grain-surface network is as complete as possible.Two reactions mainly contribute to the formation of EtOH in our chemical model: The first of these involves the barrier-mediated abstraction of an H-atom from the methyl end of the methanol molecule, by methylene, CH 2 .This is allowed either to form simply two radicals, or for those radicals to immediately recombine to produce ethanol.A branching ratio of 50% is assumed, following the approach of Garrod et al. (2022).
In reaction (5), the C 2 H 4 OH and C 2 H 5 O radicals are formed through ongoing hydrogenation of surface and bulk-ice hydrocarbons However, in the previous chemical network, only addition reactions were considered, rendering the efficient hydrogenation of unsaturated hydrocarbons inevitable and providing a bias toward EtOH production.To avoid overproduction, a number of backward reactions were also included, some of which have activation energy barriers (Table 2).

Radiative-hydrodynamics simulations
The radiative-hydrodynamic equations are solved as described by Barger et al. (2021) with an additional cosmic-ray heating term.For reference, the equations are Here, g is the gravitational acceleration, and ǫ CR is the cosmic-ray heating rate per unit mass of gas.G r and G 0 r are respectively the radiative force and heating/cooling term that couple the hydrodynamics with the radiative transfer equation.The total energy density E is related to the gas pressure P by where γ = 5/3 is the adiabatic index.The source term S (I, n) in the radiative transfer equation ( 9) is defined by where κ a,P and κ a,R are the Planck and Rosseland mean opacity, T = ρk B T /µ is the gas temperature calculated with the mean molecular weight µ = 2.33 u, and J = IdΩ/4π is the mean intensity.The rest of the symbols carry their usual meanings.
The set of equations is solved in spherical coordinates assuming 1-D symmetry in the Athena++ framework (Stone et al. 2020).The radiative transfer equation is solved as described by Jiang et al. (2019) except that the equation is solved covariantly as described by Chang et al. (2020) with additional geometric terms to account for spherical geometry following Davis & Gammie (2020).
The RHD simulation starts with a stationary cloud core with uniform density and temperature.A setup similar to model M1a from Masunaga et al. (1998) is adopted.The computational domain from r = 0.1 to 10 4 AU is resolved by 128 logarithmically spaced cells.A total mass of 1 M ⊙ with an initial temperature of T 0 = 8 K (as opposed to 10K in Masunaga et al. 1998) is enclosed in the domain.To maintain constant temperature in the initial isothermal phase, the cosmic-ray heating rate ǫ CR is set to cκ a,P a r T 4 0 where the opacity κ a,P is evaluated at the initial density ρ 0 and temperature T 0 to account for the radiative cooling.We adopted the same density-and temperature-dependent prescription for frequency-averaged opacity as in Kuiper et al. (2010).To speed up the simulation without sacrificing accuracy, the reduced speed of light approximation (see, e.g., Chang et al. 2020) with a reduction factor of 10 4 is used.Reflecting boundary conditions are employed at both the inner and outer radial boundaries for hydrodynamic variables to prevent gas from exiting the computational domain and conserve the total mass.This is a good approximation at the inner boundary because we are interested in the early phase of the star formation up to the thermally supported FHSC, before the formation of a central protostellar object.Reflecting and vacuum boundary conditions are implemented at the inner and outer boundaries respectively for radiation to both prevent radiation from leaking through the inner boundary while allowing radiation to exit through the outer boundary.The system is evolved for 255.1 kyr, up to the moment when the central temperature reaches 2000 K, at which point the dissociation of molecular hydrogen that triggers the collapse of the FHSC is expected.However, it should be noted that the chemical modeling is performed for infalling parcels for which final temperatures do not exceed 400 K, beyond which the description of the chemistry could become inaccurate, due in part to some of the measured gas-phase reaction rate coefficients used in the network exceeding their recommended temperature range.Thus, 1 AU is the minimum radius resolved within the chemodynamical model at the end of the dynamical evolution (which corresponds to a parcel of gas with a minimum radius of 2700 AU when the collapse begins).This is small enough to resolve the effective radius of the FHSC (∼5 AU).Overall, the dynamical results are consistent with Masunaga et al. (1998).Figure 1 shows snapshots of the gas density and temperature as functions of radius at four different evolutionary times in the RHD model of FHSC formation; these correspond to a moment early in the simulation, to the end-time of the simulation, and to two intermediate moments during the collapse, which becomes rapid at late times.The collapsing cloud core is initially optically thin to the thermal emission from dust grains, and it collapses isothermally for a long evolutionary time because the compressional heating rate is much smaller than the radiative cooling rate.During this time, the density develops an r −2 profile, which resembles the Larson-Penston solution for the isothermal infalling envelope (Larson 1969;Penston 1969).As the density evolves, the collapsing cloud becomes optically thick, forming the adiabatic FHSC in the central region.Once the adiabatic core is formed, the warm-up timescale is very short, showing a step-like profile for the temperature as a function of time (see Figure 2, which indicates the temperature behavior for an individual trajectory).Considering that the second collapse (which terminates the FHSC phase) begins when the central temperature reaches about 2000 K, such rapid evolution of temperature explains well the short lifetime of FHSC.This process is discussed in further detail by Masunaga et al. (1998).

Main ice constituents
Figure 3 shows the four snapshots of the radial distribution of main ice constituents at the same evolutionary times presented in figure 1, i.e. 1, 192,400, 254,800, and 255,100 yr.Abundances are shown as a fraction of total hydrogen in the gas.Although ice species are formed throughout the whole collapsing core, the chemical evolution is presented only for the infalling parcels with initial radius r init > 2700 AU; this ensures that their final radii are no smaller than ∼1 AU, and thus that their final temperatures do not exceed the nominal 400 K maximum temperature of the chemical model.
Figure 4 again shows the final spatial distribution of the main ice components, here with the gas and dust temperature profiles shown.The dust temperature in the model is consistent with the estimate toward Cha-MMS1 provided by Busch et al. ( 2020) (∼30K at 55 AU) based on ALMA dust continuum observations.The divergence between the gas and dust temperatures at radii >200 AU occurs because T gas is determined directly by the RHD simulation, in which only protostellar heating is considered, in combination with a fixed minimum temperature of 10 K; meanwhile T dust is governed by external heating, as determined by the visual extinction.The two temperatures converge when the RHD-derived temperature exceeds, and thus overrides, the extinction-derived dust temperature, at which point it is assumed that the gas density is already high enough to ensure that the two would be well coupled.In fact, in the model, this convergence occurs at a rather high density, around 10 8 cm −3 , while an explicit treatment of the gas-dust thermal coupling would be expected to achieve convergence at a substantially lower value ( 10 4.5 cm −3 ; Goldsmith 2001).However, the difference between T gas and T dust is generally small during this period of divergence, and is due to the use of only an internal heat source in the RHD model (which controls T gas ).This has a negligible effect on the dynamics at this temperature (i.e.10-15 K), while the effect on the gas-phase chemistry is also minor (as reactions with strong temperature dependencies are, in general, already negligibly slow).The focus of our model is the grain-surface and ice chemistry, thus this simplification in the gas temperature treatment at early times/low extinctions is considered acceptable.
The ice abundances within the effective radius of the FHSC (∼ 5 AU; Masunaga et al. 1998) significantly decrease due to active thermal desorption of ice species.The distribution of main ice constituents generally shows a locallypeaked structure, peaking between 10 AU and 1000 AU.This can be explained by a centrally condensed density profile combined with high temperature toward the core center.However, this locally-peaked structure is not always seen.For example, the ice abundance of CO 2 increases with radius, and the abundances of NH 3 and CH 4 in the ice mantles change little outside of the 5 AU effective radius of the FHSC.This implies the involvement of chemistry depending on ice species.
From the abundance data, ice column densities are calculated by integrating the number density (cm −3 ) of each ice species through a line of sight toward the core center (technically, 1AU offset from the core center).In this case, the protostar itself would be the background source; the column densities are thus derived by integrating the number density from the core center to the edge of the core.
Table 3 summarizes ice abundance with respect to water ice column density toward the core center.The relative abundances reproduce well the main ice composition of clearly identified species from past observations towards lowmass YSOs (Boogert et al. 2015).The relative abundances of all main ice species are within the observational range.The CO:CO 2 abundance ratio from the model (0.81) shows a good match with the observations based on median values of CO and CO 2 abundances (0.75), although the relative abundance of CO 2 with respect to water exceeds the upper quartile value from the observations.The broad consistency suggests that the Cha-MMS1 model presented here should provide a reliable basis for the chemical investigation.The peak water ice abundances with respect to total hydrogen (∼ 10 −4 ) in the model (see figure 3) are higher than required to reproduce the highest abundances seen in the observations, by a factor of around 2. This suggests that column densities of other ice species could be overproduced by a similar factor.For this reason, we further investigate the model results with a focus on the relative abundances with respect to water ice rather than on the ice abundances relative to atomic hydrogen; the possible causes of the higher ice abundances produced by the model are discussed in § 3.7.The influence of this effect on synthesized observational spectra is also tested in § 3.5.
Figure 5 shows the radial distribution of abundances for the main ice constituents (with respect to water ice) at the end of the model run, i.e., once the FHSC is formed.The following three features may be identified in the distribution: (1) An anti-correlation between CO and CO 2 at r > 2000 AU -the CO:CO 2 ratio is almost constant (∼ 0.8) at r < 100 AU, while the ratio rapidly decreases in the outer regions.This also results in the CO 2 /H 2 O ratio exceeding 1 at r 3000 AU (which is improbable given that no observational sightline has ever revealed such a ratio).The latter effect is related to the efficient conversion of CO to CO 2 (OH + CO → CO 2 + H) on grain-surfaces and within ice mantles.On the grain surface, this is a diffusive reaction due to the high dust temperatures at outer regions while nondiffusive reaction dominates within the mantles.Thus, OH radicals produced as a result of PD of water ices quickly react with nearby CO (which is one of the most abundant ice species after water).The rate at which  this instantaneous reaction mechanism occurs depends on the photodissociation rate of water ice, thus more efficient conversion appears at outer regions.Although photodissociation of CO 2 acts as a backward reaction to form CO (CO 2 → CO + H), the formation of CO 2 through the nondiffusive process dominates over its destruction.Thus, the formation of CO 2 is sensitive to the description of visual extinction in the model, as this parameter governs dust temperature and photodissociation rates at the same time.Our simple assumption for the single isolated core could cause the underestimation of visual extinctions in the outer regions, resulting in the conversion of CO to CO 2 which might be too efficient (see the discussion in § 3.4).However, the effect will have only a very minor influence on the column densities calculated below.
(2) The correlation between CO and CH 3 OH ice distribution; this indicates that CH 3 OH ice is mainly formed by ongoing CO hydrogenation on grain surfaces.
(3) Minor spatial changes in NH 3 and CH 4 ice abundances; the chemistry related to these two saturated species is fairly insensitive to the physical variations with radius.The gradual increases of their ice abundances with respect to water ice seen at outer radii is caused by the faster accretion of NH 3 and CH 4 than H 2 O, as the gas-phase budget of C and N has not run down yet at those outer radii (Garrod et al. 2022).We have investigated other volatiles expected to be abundant in the ice mantles, including CO-hydrogenation intermediates and well known O-bearing hot-core species; species of interest are listed in Table 4, along with their abbreviations.Figure 6 shows the spatial distribution of these ice species with respect to the local fractional abundance of water ice, X(H 2 O), at the end of the evolution.The two upper panels show large (n atom ≥ 5) O-bearing species that have been detected toward hot cores while the lower panel shows the intermediates related to CO hydrogenation.

Other volatiles expected to be abundant in the ice mantles
The results from this model predict high ice abundances for the CO-hydrogenation intermediates, including radicals ( 10 −3 ×X(H 2 O)) although this value is considerably lower than the upper limit of the storage capacity of the ice for reactive radical species (< 10%; Golden 1958;Jackson & Montroll 1958;Jackson 1959a,b).Garrod et al. (2022) further discuss the abundances of radicals in the ices.For more complex O-bearing species, their solid-phase abundances are typically on the order of 10 −5 -10 −4 ×X(H 2 O), while EtOH shows distinctively high solid-phase abundances of 10 −3 ×X(H 2 O).The five O-bearing species presented in the upper-left panel of figure 6 also show relatively high abundances.Observations of these species have, either directly or indirectly, implied their presence within the ice mantles.Three of them (ACTD, MF, and DME) have been identified in the gas phase in cold, prestellar environments (Bacmann et al. 2012;Vastel et al. 2014;Cernicharo et al. 2012;Scibelli & Shirley 2020), with grain-surface formation being suggested as an explanation.In addition, the analysis of the mid-IR spectrum of the massive protostar W33A -one of the few sources for which a high quality mid-IR spectrum is currently available -identified three prominent features that have been attributed to EtOH, FA, and ACTD, although their identification in the ices remains uncertain (Schutte et al. 1999;Boogert et al. 2008;Öberg et al. 2011).Terwisscha van Scheltinga et al. ( 2018) re-analyzed the same spectrum with more accurate spectroscopic laboratory data, deriving upper limits for ice abundances with respect to water of 1.9 % and 2.3 % for EtOH and ACTD, respectively.This is 1-2 orders of magnitude more than predicted by our model.AA and MM have not been identified from the ice observations so far, but their large abundances in the model suggests them to be plausible candidates for solid-phase detection, their spectral characteristics notwithstanding.In particular, those species could be enhanced at the outer radius (in low-density regions) by almost one order of magnitude compared to the core center.
We categorize the spatial distribution of the species shown in figure 6 into three types: 1) pudding shape, i.e. flattened center, decreasing at large radii; 2) bowl shape, i.e. flattened center, increasing at larger radii; and 3) locallypeaked, with X(i)/X(H 2 O) peaking at a few thousand AU.Table 5 shows the shape with which we identify each of the ice species shown in the figure.Interestingly, aside from the three aforementioned COMs (ACTD, MF, and DME) detected in cold environments, the chemical species following each type of distribution can be linked back to a specific radical or molecule involving their formation.For example, the CO ice distribution is a pudding shape, and most of the ice species whose formation is directly related either to the CO→FA conversion pathway (see equations 1 and 2) or to the hydrogenation products of CO show similar distributions.
As for the larger species, whose formation pathways are more complex, the abundance profiles are nevertheless found to depend on which functional groups are present.All species with bowl-shaped distributions contain a methyl  group ( -CH 3 ), while all species containing the hydroxymethyl radical ( -CH 2 OH) exhibit locally-peaked distributions peaking at a few thousand AU.ACTD, MF, and DME show locally-peaked distributions despite each containing a methyl group, but in the test model of a moderately higher extinction environment their distributions converge to the bowl-shape distribution.This result shows that -CH 3 (which perhaps may include -CH 3 O) and -CH 2 OH might be the key functional groups determining the spatial profiles of some solid-phase species.A recent high-resolution observation indirectly implies this; gas-phase observation of O-bearing COMs in Orion KL show different morphology depending on which functional group is present in the molecule (Tercero et al. 2018).To explain this, Tercero et al. (2018) propose that different radicals, methoxy ( -CH 3 O) and hydroxymethyl ( -CH 2 OH), may dominate, driving the different chemical complexity of the regions.
Another main factor that determines the radial distribution of ice species is the A V distribution, dominantly during the pre-model.At the beginning, the visual extinction ranges from 3.0 mag (innermost) to 2.5 mag (outermost), and it evolves to the ranges from 6.2 mag to 2.7 mag with the density evolution during the pre-model.The chemical differentiation of many species (e.g.EtOH) are determined by the balance between formation and destruction processes that both become more efficient further out from the central core.In outer regions, low A V and/or higher dust temperatures can drive the formation of these species through PDI nondiffusive reactions and/or diffusive reactions, while destruction via PD also becomes active.However, the formation of some species such as MF becomes rather efficient when the density is higher, indicating that the A V has increased as well.The formation of this type of species is not strongly dependent on the PDI processes in the bulk ice, but on the build up through direct surface production via 3-B nondiffusive reactions, when radicals become more abundant on the grain surfaces (see figure 13 in Garrod et al. 2021).

Calculated ice column densities
Figure 7 shows the (absolute) column density profiles produced by our model of Cha-MMS1 for large solid-phase molecules.These are calculated by integrating the number density of each species along parallel lines of sight intercepting different offsets from the core center.For comparison, the gray dashed lines indicate the observed radial offsets from the center of Cha-MMS1 of the lines of sight toward a selection of background (BG) stars.The information of these BG stars is listed in table 6.These lines of sight are showing high-A V in the Ice Age field.Considering the high A V threshold expected for the detection of complex ice species, these lines of sight would be the most probable to identify icy COMs, of those BG stars yet known that lie near on the sky to Cha-MMS1.The two outer BG stars (orange crosses in figure 8 -NIR 38 and J110621.63-772354.1 -will be observed both with NIRSPEC's fixed slit mode (0.6 to 5 microns) and MIRI's low resolution spectroscopy (LRS) mode (5-14 microns) which operate in the near-IR and mid-IR, respectively.The core center and the other two inner positions (pink crosses in figure 8)-J110637.01-772323.0 and J110637.34-772352.7 -will be observed with NIRCAM / Wide Field Slitless Spectroscopy (WFSS) only (2.4 to 5 µm).The combination of these five lines of sight will provide a spatial profile of ice composition possibly including COMs.
The extinction towards the two MIRI LRS targets, NIR 38 and SSTSL2J110621.63-772354.1, was calculated for the Ice Age ERS program using the observed photometric color taken between the K bandpass at 2.2 µm and Spitzer IRAC band 1 at 3.6 µm taken from Persi et al. (2001) and the IPAC Gator database.The calculation assumed an intrinsic K-IR1 color for background giant stars of 0.03 from Majewski et al. (2011) and the highest extinction A λ /A V curve for molecular clouds from McClure (2009).The other two background stars were only observed at WISE bands 1 and 2. Therefore we assumed an intrinsic background giant W1-W2 color of -0.08 from Schlafly et al. (2016) and the A W1 /A W2 extinction of 0.463 from Xue et al. (2016).From comparing the variation in A V derived for other, less reddened stars from both this set of colors and other colors, we estimate a conservative 20% systematic uncertainty in our final A V measurements.
To see how well the density structure described in the model matches with the observations, an observational column density of total hydrogen toward each background star is compared with a calculation of the column density along a line of sight in the model that lies at the appropriate radial separation from the core center.A similar comparison is made between the observed column density of total H toward the center of Cha-MMS1 and the same line of sight toward the model of the source.
For each line of sight, the column density of total hydrogen, N H,obs , is calculated based on the observational visual extinction information, using the relation N H,obs =1.8 × 10 21 A V cm −2 .To take some account of observational effects in the comparison, the corresponding model values, N H,model , are obtained by convolving the column density profiles with the spatial resolution of Spitzer IRAC (2 ′′ ), using a Gaussian beam profile.
For the core center, N H,obs derived from 870 µm dust continuum is also available from the literature (Belloche et al. 2011).This value is compared with the model by convolving the modeled column density profile with the beam size of the dust continuum map (θ HPBW =21.2 ′′ ).As seen in table 6, the column density derived from the model is roughly consistent with the observation at the core center.This result remains consistent even when an alternative estimate of the A V /N H factor is employed.Based on the study by Chapman et al. (2009), around two times lower conversion factor (N H,obs =1.1 × 10 21 A V cm −2 ) might be appropriate for the high density molecular gas.The employment of this alternative conversion factor results in around two times lower N H,obs , which in fact agrees even better with the model at the core center.However, the disparity in total column density between model and observation grows at larger radii, where N H,model falls more than two orders of magnitude below the observational value for the outermost BG position.
The main reason for this discrepancy is that the physical evolution of Cha-MMS1 in the model is described as a single isolated object with the RHD simulation, while Cha-MMS1 is actually one of nine young, low-mass star-forming cores within a distance of 0.2 pc in the plane of the sky (Maureira et al. 2020).For example, another Ice Age targetprestellar core C2 (right white cross in figure 8)-exists near the outermost BG position, and it is even closer to this BG position than Cha-MMS1.This explains the much higher column density toward this line of sight in the observation than the model.Thus, the model is unlikely to produce meaningful predictions for ice abundances for the nearest BG stars yet detected in the vicinity of Cha-MMS1.However, it is expected that the JWST observations will reveal new background stars that are closer to the core center and therefore may provide a better constraint on the model.
In addition to the density structure, we also compare the modeled water-ice column density with values previously reported for observations.There is no direct measurement of water ice column density for this source, but values toward various other low-mass YSOs are available for comparison.This allows the compatibility of the modeled ice compositions with the observations to be determined.Considering the correlation between water ice column density and visual extinction, which has been found for a range of visual extinctions 0 < A V < 30 mag (reviewed in Boogert et al. 2015), water-ice column densities on the order of 10 19 cm −2 would be expected, given the high visual extinction of Cha-MMS1.The corresponding value from the model is simply derived by integrating the number density of water ice along a pencil beam toward the core center (In observations, the water ice column density is generally measured from the optical depth of the 3 µm band relative to the dust continuum, which is expected to be confined within a very small radius for a deeply embedded object).This calculation results in the value of 1.2×10 21 cm −2 , which Note-The SynthIceSpec are assuming Gaussian for each absorption feature for simplicity.We assumed the center and the half of the integrated interval for the derivation of the integrated absorbance as peak position and FWHM of the absorption band, respectively.
is around two orders of magnitude higher than the expected observational one.This discrepancy could be caused by the model having too high a local density near the core center in particular; the N H,model value obtained without beam convolution is overproduced by a similar factor as the water column density, while the maximum water ice fractional abundance of ∼ 10 −4 in the model (see Figure 4) is just a factor of 2 higher than the measured abundances previously reported for observations (Boogert et al. 2013).
Figure 4 indicates also that the modeled dust-grain ice mantles are present to an inner radius of just a few AU, with the temperatures inside this radius being too great for ices to be maintained on the dust grains.A more extended final temperature profile would push the ices out to a region of lower density, reducing the column densities of solid-phase species.Test calculations that curtail the line-of-sight integration of ice column densities to a minimum inner radius on the order of 100 AU indeed reduces the calculated values by an appropriate factor.
These proposed influences on density and temperature could be explained by the presence of a protostellar outflow carving out the central region of the core, as long as the sightline passed through the outflow cavity before reaching the central star.Furthermore, much of the material settles in a disk even at this early stage.The observations typically don't cover the high disk columns because edge-on disks have such high columns that we don't detect the central star anymore.Such a disk and an outflow has actually been identified toward Cha-MMS1 (Busch et al. 2020), but is not included in the 1-D dynamical model used here.
The possibility that the chemical model itself is still overproducing the ice abundances in general cannot be completely ruled out, and the simple method of column density calculation used here could be related to this discrepancy as well.Further discussion of this topic is provided in § 3.7.Here, for the purposes of further comparison of the modeled ices with observations, we make a further uniform, empirical adjustment to the modeled ice abundances, scaling all the modeled ice column densities down by a factor of 100, i.e. f = 0.01.The water-ice column density after scaling roughly corresponds to a value that can be derived via extrapolation of Figure 7 in Boogert et al. (2015) to 100 mag (1.9×10 19 cm −2 ).

Synthesis of model ice spectra and diagnostic detectability of new ice species with JWST
Ice absorption spectra will be obtained toward the core center of Cha-MMS1 with the NIRCAM/WFSS instrument, as a part of the Ice Age program.In preparation for these observations, the near-IR spectrum expected from the core center is synthesized based on the modeled ice composition.By doing so, we are able to provide a practical diagnosis of the detectability of new ice species with JWST, as well as to get a sense of the order of magnitude of solid-phase COM abundances that would be required in order for them to be identified.A simple synthetic ice-spectrum simulation script has been developed for this purpose, with input data provided by members of the Ice Age team.The main inputs for the simulation are the retrieved modeled ice column densities and three spectral parameters for the characterization of each absorption feature (i.e.frequency, width, and apparent band strength A ′ ).
Until recently, laboratory studies for the quantitative measurement of IR spectra of interstellar ice COMs had scarcely kept pace with astronomical discoveries.Fortunately, more and more effort has been expended on this type of study in recent years, although it has typically been focused on the mid-IR regime as this is considered the fingerprint region for solid-state vibrational spectroscopy (see e.g., Terwisscha van Scheltinga et al. 2018;Hudson et al. 2018;Rachid et al. 2020;Hudson & Ferrante 2020;Hudson et al. 2020;Terwisscha van Scheltinga et al. 2021).For example, Terwisscha van Scheltinga et al. ( 2018) present ice spectra (2.5µm -20 µm) for ACTD, EtOH, and DME in astronomically relevant environments, and they characterize the key absorbance features at λ > 5µm from these reference spectra.Hudson et al. (2018), Hudson & Ferrante (2020) and Hudson et al. (2020) provide direct measurements of apparent band strengths and absorption coefficients from the mid-IR spectra for amorphous ices of ACTN, DME and ACTD, respectively.From the aforementioned literature, we could directly collect the absorption feature characteristics relevant to NIRCAM/WFSS coverage (2.4µm -5 µm) for ACTN, ACTD and DME.It should be noted that ice features vary considerably depending on the chemical composition of the ice mixtures.However, there is no clear consensus on the representative ice mixture for the interstellar icy mantles, and the experimental setup for ice mixtures in the past varies significantly.For consistency, we chose the measurements from pure ice spectra as our standard.In the case of EtOH, only the band information of absorption features at λ > 5µm are readily available.However, the reference spectrum of pure EtOH (2.5µm to 20µm) ice is publicly available on the Leiden Ice Database 1 .This allowed us to derive the spectral parameters for the shorter wavelength region (2.5µm -5µm).The unknown apparent band strengths were estimated by using the ratio between the integrated areas of an absorption feature with different absorbance within a single IR spectrum.
In general the column density of a particular ice species, N species , is determined according to where band I λ is the integrated absorbance of the IR band under consideration and A ′ λ the apparent band strength.
The ice column density derived from different absorption features in a single IR should correspond to each other.Thus the unknown band strengths in the near-IR region for EtOH can be constrained by using the representative absorption band for EtOH at 9.514 µm as a reference (Terwisscha van Scheltinga et al. 2018).Table 7 summarizes the spectral parameters used in this simulation.
Ice spectra derived from the chemodynamical models are calculated for the main simple ice constituents and the four aforementioned COMs for which IR bands are characterized.Two sets of ice column densities are tested: the regular ice composition (setup A) which is derived by directly integrating the spatial distribution of the species from the core edge to the core center, and the ice composition scaled down by factor f = 0.01 (as discussed in § 3.4; setup B.1). Table 8 shows the two sets of ice composition for the synthesis of spectra.
13 CO and 13 CO 2 are known as two rarer isotopologues of which ice absorption features are clearly detected in interstellar environments.As the chemical model considers the chemistry involving the main isotopologues only, the contribution of the two aforementioned species to the IR absorption features is considered by scaling the column densities of 12 CO and 12 CO 2 with the ratios as follow; 12 CO/ 13 CO=70 and 12 CO 2 / 13 CO 2 =86.The former is adopted from Boogert et al. (2002), and the latter is derived from the solid state 13 CO 2 / 12 CO 2 ratio as a function of Galactocentric radius (Boogert et al. 2000).
Figure 9 shows the IR spectra for the NIRCAM/WFSS coverage for the two different setups of the ice composition.For each setup, a control spectrum is additionally synthesized (blue dashed line) assuming zero contribution of COMs larger than CH 3 OH to the spectrum for comparison.The IR band at NIRCAM/WFSS coverage is dominated by water features and a mixture of methyl (CH 3 ) functional group vibrational stretching modes from the different ice constituents in the simulation.In particular, the spectrum of regular ice composition without scaling (setup A -left panel) is highly saturated by these features, making it hard to identify IR absorption peaks around 3.4 and 3.6 µm, which are distinctively seen in setup B.1 (right panel).However, regardless of whether or not the ice composition is scaled down, the contribution of icy COMs to the IR absorption is either negligible or would not be distinguishable in the case of real observations.To determine the approximate COM ice abundances required for them to be detected observationally, we again simulate ice spectra based on setup B.1, but this time with the ice abundances of the four COMs increased with respect to the water-ice column density by a factor of 2.5 (setup B.2). Figure 10 shows the near-IR absorption spectra for setup B.2.The absorption from COM ices is more notable in this case, adding up to three absorption peaks at 2.8, 3.4 and 3.6 µm to the control IR absorption spectra.However, even if the ice abundances of EtOH are increased to as high as 1% with respect to icy H 2 O in this case, the COM absorption features are still minor.Furthermore, the absorption peaks overlap with the absorption IR band of the main ice constituents.This would introduce additional uncertainty in the identification of COM features in actual observational spectra.
Although only near-IR spectra will be obtained toward the core center as part of the Ice Age project, we may extend our investigation to the mid-IR regime (5µm -14 µm; MIRI/LRS coverage) with the same ice composition to seek any distinctive ice feature from COMs; MIRI observations will be performed towards the two high-A V lines of sight located in the outer envelope of Cha-MMS1.
Figure 11 shows the resulting mid-IR absorption spectra with setup B.2 (COM-rich).The IR band at 9.2 µm, which is mainly attributed to the CH 3 rocking mode of EtOH, could be a prospective IR band for COM detection (Terwisscha van Scheltinga et al. 2018).Although it is unclear whether such high abundances of COMs might be achieved in typical star-forming environments, the production rates of COMs can change sensitively with UV flux.For example, in a test model with a higher extinction environment (A V,bac = 5 mag; not shown) EtOH ice abundance decreases by around a factor of four compared to the standard model.When the additional background visual extinction is not included at all (A V,bac = 0 mag; not shown), the EtOH ice abundance significantly decreases again as photodissociation overwhelms the synthesis of the species via PDI.Considering all this, the identification of COM ice might be possible with MIRI observations in an environment with a modest UV.Similar variations in COM abundances might be expected if an alternative N H :A V ratio than the one used here were adopted (see also Sec. 3.4).
It should be noted that dust features were not considered for our calculations but only absorption features from species within ice mantles.Silicates will produce deep dust-continuum baselines around 9.7µm (see L1014 IRS in figure 1; Öberg et al. 2011), and this could make the absorption features of COMs from 8-11 microns more difficult to detect.Here we assume that the dust continuum baseline will be constrained well and already subtracted.By observing multi lines-of-sight within a single source with the Ice Age project, we expect the dust continuum baseline would be determined with unprecedented accuracy.MIRI is expected to have at most a signal-to-noise ratio of 300.This corresponds to the detection of an absorption feature at 1% with respect to the continuum flux at 3σ.

Gas-phase abundances
Gas-phase molecules have been detected in many single-dish observational studies of Cha-MMS1.For example, Kontinen et al. (2000) observed a selection of simple species, carbon-chain molecules, and methanol using the Swedish ESO Submillimeter Telescope (SEST), while Tennekes et al. (2006) determined the HCN/HNC relative abundance ratio by mapping the HCN and HNC emission with the same instrument.Belloche et al. (2006) measured the deuterium fractionation for molecular ions with the APEX telescope, and Cordiner et al. (2012) performed 7mm observation with the ATNF Mopra telescope for a selection of species including polyynes, sulphuretted carbon chains and methanol.
Although the chemical modeling presented here is focused on solid-phase abundances, the models indeed produce fractional abundance profiles for gas-phase species.In order to compare these with the observational values, column densities must be produced.Molecular column densities through the model core are produced for various offsets from the core center.These raw values are then convolved with a beam size appropriate to the observations, with the beam centered on the core.Cordiner et al. (2012) provide the most extensive dataset of observed molecular lines.However, the column densities derived from our model are significantly lower than those in the Cordiner et al. dataset, by more than one order of magnitude.The large observing beam size of the Mopra telescope is attributed to this discrepancy; the observational beam size (96 ′′ at 36 GHz and 77 ′′ at 45 GHz) is comparable to the final radius of the outermost trajectory (r≃8660 AU).However, the two nearby sources -Ced 110 IRS4 (class I protostar) and C2 (prestellar core) -contaminate the observing field (see figure 8), contributing to the higher molecular column densities, while the model assumes an isolated single FHSC.This may be a major factor in the poor comparison between models and observations with such a large observing beam size.
Therefore, the abundances measured with a smaller beam should be better reproduced by the chemical model.For example, the SEST observation (FWHM=55 ′′ at 91GHz) shows that the abundances of some key species are roughly consistent with the modeled abundances within one order of magnitude (Table 9).Methanol is in fact well reproduced, and HCN and HNC are quite adequate.Sulfuretted carbon and cyanopolyyne HC 3 N abundances are poorly reproduced; however, the overall abundance of sulfur (in whatever form) is poorly constrained in astrochemical models, as discussed by Laas & Caselli (2019).Furthermore, HC 3 N is traditionally considered an "early-time" molecule that becomes prevalent while substantial amounts of atomic carbon are still available in the gas phase; conversely, our model specifically corresponds to a more chemically-evolved state of the gas; thus the substantial abundance of this molecule could in practice originate from a spatial scale beyond that included in our model.Ultimately, an accurate comparison of gas-phase species toward the Cha-MMS1 source may require high-resolution observations.Although there have been a few studies performing ALMA observation toward this source, it is more focused on the dynamics of the outflow, and/or quantitative measurements of the chemical composition were not fully available (Busch et al. 2020;Allen et al. 2020;Maureira et al. 2020).
Using the gas-phase CO abundance reported by Kontinen et al. (2000) for Cha-MMS1, Cordiner et al. (2012) derived an observational depletion factor of 2, assuming an undepleted CO abundance of 9.5 × 10 −5 (Crapsi et al. 2005); this value is rather low compared with the typical CO depletion factor reported toward prestellar cores (Crapsi et al. 2005).By contrast, the CO depletion factor suggested by our model is ∼10, assuming a convolving beam size of 50 ′′ (the same beam size used in the measurement of N (H 2 ) for the derivation of CO abundance; Tennekes et al. 2006), which is much more consistent with typical observational values for other sources.The discrepancy with the Cha-MMS1 result is likely due to the assumed beam size used in the CO observations, which is not available in the literature.The uncertainty in this value renders a comparison of observational and modeled CO depletion factors of only limited value.

High ice column density of H 2 O
As discussed in § 3.4, the high density regions at the very center of the core may be particularly poorly represented using our simple assumption for the density structure and integration technique.The difference of ice column densities between the model and the observation could indicate the presence of an outflow cavity or other structural differences that are poorly represented by a one-dimensional dynamical treatment.Furthermore, this divergence between the model and generally-observed values may also be related to the longstanding mystery of the oxygen budget.The water abundance (gas and ice) estimated from observations of various starforming regions is much lower than the one expected from the known oxygen elemental abundance.This conundrum is reviewed by van Dishoeck et al. (2021), who discuss two possible explanations: (i) The missing oxygen is locked up in an unidentified refractory component.If this is the case, the high abundance of water in the model is caused by a lack of gas-grain chemical network for the unidentified refractory component, which results in the oxygen budget converging to the formation of water ice.If one assumes that all volatile oxygen is locked in water, the water abundance is expected to be around 4.0×10 −4 (van Dishoeck et al. 2021).This value is just a few times higher than the water-ice abundance in our model in cold environments (∼1.0×10 −4 ).(ii) Water ice abundance estimated from ice observations can be underestimated if the missing oxygen is locked up in larger µm-sized grains that do not contribute to infrared ice absorption.Much observational evidence for grain growth to µm sizes has been published (e.g.Boogert et al. 2013); the toy model with two dust populations (0.1 µm-and 10 µm-sized grain) shows that a fraction f large = 0.01% of the grain population by number in 10 µm-sized grains would catch 50% of the atomic oxygen during freeze-out and make it invisible (van Dishoeck et al. 2021).However, one argument against water ice being hidden on very large grains is that after sublimation in hot cores or grain sputtering in a shocked region, the gas phase H 2 O abundance along with the abundance of other O-bearing species does still not add up to the cosmic O abundance (van Dishoeck et al. 2021).Possibly the 40 µm H 2 O lattice mode can shed some light on this, though it lies beyond wavelength range available to JWST.

The effect of nondiffusive chemical mechanisms
Past astrochemical models have relied on the thermal diffusion of surface or bulk-ice radicals to drive reactions that form COMs (Garrod & Herbst 2006;Garrod et al. 2008;Garrod 2013).However, recent detections of COMs in prestellar environments imply that COMs can be efficiently formed before star-forming regions heat up.In such colder environments, alternative mechanisms not involving thermal diffusion can become important.For example, COMs may be formed nondiffusively as the result of surface radical production in close proximity to pre-existing radicals (Chang & Herbst 2016;Jin & Garrod 2020;Garrod 2019;Garrod et al. 2022).The chemo-dynamical model used here employs the two nondiffusive mechanisms -three-body (3-B) (Jin & Garrod 2020) and photodissociationinduced (PDI) mechanisms (Garrod 2019).In this section, we discuss the effect of those nondiffusive mechanisms on the composition of icy mantles by turning off either or both mechanisms.
Table 10 summarizes the results of the analysis for the main ice constituents.The relative abundances in this table are derived from the ratios of column densities, which are calculated by integrating the number density of each ice species through a line of sight toward the core center.Both nondiffusive mechanisms transfer more of the oxygen budget from CO to CO 2 , through the reaction OH + CO → CO 2 + H, which results in a better reproduction of generic observations.The modeled ice abundance of CH 4 significantly increases with the PDI mechanism switched off.This implies that with the PDI mechanism, a significant fraction of PD products of methane ice such as CH 3 on the grain consecutively react with nearby ice species to form larger species rather than reforming CH 4 , while the effect of the three-body mechanism is negligible in this case.This PDI production of methyl-group bearing molecules occurs mainly in the ice mantle, due to the much larger quantities of material stored beneath the ice surface.The solid-phase abundances of CH 3 OH and NH 3 barely change regardless of the choice of nondiffusive mechanism, reproducing the observed abundances well.Aside from the above-mentioned effects, the influence of nondiffusive mechanisms on the main, simple ice constituents is very minor.
Table 11 shows the results for solid-phase O-bearing COMs.As seen in the last row of this table, the inclusion of the nondiffusive mechanisms increases the abundances of COMs in general, as it provides additional pathways to form large compounds without a temperature dependence.In particular, the abundances of the three structural isomers of C 2 H 4 O 2 (MF, GA, and AA) and two CH 2 OH-bearing COMs (EG and MM) increase by more than three orders of magnitudes with the inclusion of any of the nondiffusive mechanisms.It is notable that the change in chemical behavior produced by such changes can be non-linear; when one strong mechanism is removed, a weaker one may take over to replace (partially) what would otherwise be done by the stronger mechanism.Simple comparison of these models therefore provides only a limited indication of which mechanism is most influential in the development of chemical complexity in general, as the chemical behavior of the species in the various nondiffusive setups is different.For example, MF is more efficiently formed via 3-B excited formation (see Jin & Garrod 2020, for further detail) than via the PDI mechanism, although either mechanism will promote the formation of this species.The ice abundance of EtOH is higher when only the PDI mechanism switched on rather than the both, while FA and ACTD are the opposite cases; the ice abundances of the two COMs decrease with the PDI mechanism switched on.This is because OH and CH 2 radicals -which are the photodissociation products of H 2 O and CH 4 , respectively, at early times in the model evolution -are consumed in various nondiffusive reactions, to form larger species.The main formation routes for FA and ACTD involve OH and CH 2 radicals, respectively: FA is synthesized either by the reaction CO + OH → COOH, followed by hydrogenation, or by the addition of OH and HCO radicals.The two-stage reaction between CH 2 and H 2 CO, in which CH 2 abstracts an H-atom from formaldehyde with the product radicals then immediately recombining, is one of the main formation routes for CH 3 CHO.Thus, the PDI mechanism is critical to transfer the CNO budget to various species at early times, developing the chemical complexity of the bulk ice.
3.9.Implications of this study and future perspectives Motivated by the imminent observational improvements to be provided by JWST, laboratory experiments have recently been undertaken to characterize the spectroscopic behavior of larger species in the ice (Terwisscha van Scheltinga et al. 2018;Rachid et al. 2020;Hudson & Ferrante 2020;Hudson et al. 2020;Terwisscha van Scheltinga et al. 2021).As these quantitative measurements have accumulated, the demand has increased for the application of the experimental outputs into a chemical model; this combined method can more practically probe the strategy needed to identify new ice species with JWST.As ice spectra are influenced by the ice composition as well as the physical structure of the solid (ices and dust grains), it is important to obtain a plausible estimate for the ice composition with a chemical model.This initial study begins to meet this demand.However, our results show that even if ice abundances of COMs are predicted to be as high as 0.5% with respect to water ice, the absorption features of these species are not distinguishable from the water absorption features.This is a good example highlighting the importance of a comprehensive study: the detectability of the ice species can be securely diagnosed with careful consideration of spectroscopic characteristics as well as molecular content before observations follow.
Based on the combination of outstanding features such as a relatively strong IR band, high ice abundance, and the potential to be enhanced in high UV-flux environments, this study indicates EtOH as a good candidate for a new ice detection with MIRI observations, given an appropriate ice abundance.So far, EtOH has not been securely identified in ice observations, but this species is one of the three COMs (MF, ACTD, and EtOH) that is considered to be a tentatively identified species (e.g.Terwisscha van Scheltinga et al. 2018Scheltinga et al. , 2021)), which supports the finding of this study.Only gas-phase detections of EtOH in hot core regions has been reported in the Milky Way, indirectly implying a possible origin for this species in the sublimated ices.However, if EtOH were to be detected in cold gas at a cloud edge, perhaps as the result of photodesorption from the icy grain surfaces, this would separately indicate the richness of EtOH in the ice.Such an observational probe could be a reasonable diagnostic to find candidate regions for new EtOH ice detection with JWST in the future.
The results of this study show that the detection of new ice species may be demanding even with JWST.However, it is noteworthy that the IR spectroscopic considerations in this study were only possible for a handful of COMs, due to a lack of quantitative measurements of IR data.Spectroscopic behavior can be significantly different between species, and some molecules that are not covered in this study yet could be "IR-bright" (strong IR absorption band).For example, a recent study found that the infrared absorbance of propynal's alkyne bond is about 30,000% stronger than the corresponding feature in acetylene (Hudson & Gerakines 2019).Although propynal's ice abundance is not high enough in our FHSC chemical model, further IR investigation of other larger interstellar species may unveil new IR bright COMs that are also abundant in the chemical model.With the interplay between laboratory work and chemical modeling studies, physical environments where such IR-bright COMs are efficiently synthesized can be also examined.This comprehensive study provides firm ground to further explore new ice species in the future.

CONCLUSIONS
To timely diagnose the detectability of new ice species with JWST, we predict the chemical compositions within the dust-grain ice mantles in the young star-forming core Cha-MMS1 using a 1-D chemo-dynamical model.The expected IR spectrum based on the given ice composition is simulated.The main conclusions of this study are enumerated below: 1.In the chemo-dynamical model, the relative abundances of the main ice constituents with respect to water toward the core center match well with generic observational values, providing a firm basis to further explore the ice chemistry.
2. The model predicts relatively high solid-phase abundances (i.e.> 0.01% with respect to water ice column density) of seven large oxygen-bearing species (ethanol, acetaldehyde, methyl formate, dimethyl ether, formic acid, acetic acid, and methoxy methanol).All but acetic acid and methoxy methanol have, either directly or indirectly, implied their presence within the ice mantles through the past observations.
3. For the large molecules, the abundance profiles in the model are found to depend on which functional groups are present; -CH 3 and -CH 2 OH might be the key functional groups delineating the spatial distributions of species.
4. The IR spectrum is synthesized based on the modeled ice column densities of four COMs where IR band information is known, as well as the main ice constituents.The contribution of COMs to IR absorption is minor compared to the main ice constituents, making the identification of COMs with the NIRCAM/WFSS instrument unlikely.
5. Mid-IR observations of COM-rich (ice abundances of COMs w.r.t. the water ice column density > 1%) environments could reveal distinctive ice features of COMs.As this ice feature overlaps with silicate features from dust grains at around 9.7µm, accurate treatment for the continuum subtraction is critical.
6.The results of this study show that the detection of new icy COMs with JWST could still be challenging.However, Further IR investigation of other larger interstellar species may unveil new IR bright COMs that are also abundant enough within the ice mantles to be identified.

Figure 1 .
Figure 1.Snapshots of the physical conditions as functions of radius at four different evolutionary times in the RHD model ; 1 yr (blue), 192,400 yr (green), 254,800 yr (yellow), and 255,100 yr (red)

Figure 2 .
Figure 2. Physical evolution of the innermost trajectory as a function of time during the RHD model.The black solid line and the yellow dashed line indicates the evolution of density and temperature, respectively.

Figure 3 .
Figure3.Radial distribution of the main ice constituent abundances during the evolution of Cha-MMS1; H2O (solid black), CO (solid green), CO2 (dashed green), NH3 (solid red), CH4 (dashed red) and CH3OH (solid blue).The yellow lines indicate density profile within the source obtained from the RHD simulation.

Figure 5 .
Figure5.Final distribution of the main ice constituent abundances w.r.t.water ice within Cha-MMS1; H2O (solid black), CO (solid green), CO2 (dashed green), NH3 (solid red), CH4 (dashed red) and CH3OH (solid blue).The yellow lines indicate the density profile within the source determined by solving the RHD equations.

Figure 6 .
Figure 6.Radial distribution of the ice abundances with respect to water ice abundances at the end of the evolution.The yellow lines indicate the density profile within the source determined by solving the RHD equations.

Figure 7 .
Figure 7. Radial distribution of the column density produced by our model of Cha-MMS1 for large O-bearing ice species.Gray dashed lines indicate the observed radial offsets from the center of Cha-MMS1 for the lines of sight toward a selection of background stars.

Figure 8 .
Figure 8. Map field of ice mapping in the ice age project.Green contours indicates high-AV 870 µm dust continuum emission (Belloche et al. 2011).The targets of ice age are denoted as large white crosses.From the left, Ced 110 IRS4 (Class I protostar), Cha-MMS1 (Class 0 protostar), and C2 (prestellar core).

Figure 9 .
Figure 9. Simulated IR spectrum at NIRCAM/WFSS coverage.The left panel shows the one based on the regular ice composition without scaling (setup A) while the right panel shows the one from scaled down ice composition(0.01×Ni;setup B.1)

Figure 11 .
Figure 11.Simulated IR spectrum at MIRI coverage when 2.5 times higher ice abundances for COMs are assumed (setup B.2) The left panel shows the spectrum including all ice species, while the right panel has the contribution by the main ice constituents subtracted.The grey dashed line indicates the position of 9.7µm silicate dust absorption band which is assumed to be already subtracted.
model development was funded by the National Science Foundation through the Astronomy & Astrophysics program (grant No. AST 19-06489).Development of the coupled chemical and dynamical modeling treatment was funded by the NASA Astrophysics Theory program (grant number 80NSSC18K0558) and is supported in part by NASA NSSC20K0533, NSF AST-1815784, and NSF AST-1716259.

Table 2 .
Updated and/or newly-included grain-surface and ice-mantle reactions

Table 3 .
The relative abundances of the main ice constituents with respect to water ice column density a Öberg et al. (2011); For each molecule, the second row gives the median and lower and upper quartile values of the detections, and in brackets the median including upper limits.The third row gives the full range of abundances.

Table 4 .
Ice species investigated in this study

Table 5 .
Spatial distribution of ice abundance of COMs w.r.t.water ice within Cha-MMS1

Table 6 .
The information on background stars near to Cha-MMS1 This estimation is based on the A V towards the surrounding background stars.There is a systematic uncertainty of ∼ 20% on this estiamtion.See § 3.4 for more details.b Belloche et al. (2011) -N H,obs is derived from the 870 µm dust continuum emission map.c only two WISE photometric bands are available for the estimation of A V .

Table 7 .
Ice spectral parameters assumed for the synthesis of ice spectrum

Table 8 .
Ice compositions for the synthesis of IR spectra

Table 9 .
Comparison of the abundances of gas-phase species between the model and the observation

Table 10 .
The effect of the nondiffusive mechanisms on the main ice constituents

Table 11 .
The effect of the nondiffusive mechanisms on the O-bearing ice species