NAUTILUS Gas-Grain Chemical Model

• NAUTILUS Gas-Grain Chemical Model is a comprehensive astrochemical framework that integrates gas-phase reactions with surface and mantle chemistry using two-phase and three-phase formulations.
• It employs rate-equation methods with explicit reaction–diffusion–evaporation competition, accurately modeling adsorption, desorption, and surface processes to track molecular evolution.
• The model is flexibly coupled to various physical structures, including cold cores, disks, shocks, and collapsing protostars, providing critical insights into interstellar chemical dynamics.

Searching arXiv for core NAUTILUS gas-grain model papers and recent applications. NAUTILUS is a time-dependent gas-grain astrochemical modeling framework used to evolve molecular abundances in interstellar and circumstellar environments by coupling gas-phase chemistry with adsorption, grain-surface chemistry, desorption, and, in many implementations, bulk-ice chemistry. Across the literature, the name refers not to a single immutable configuration but to a family of related rate-equation models that have been deployed in two-phase and three-phase forms, extended to multi-grain populations, deuteration and spin chemistry, shock chemistry, and specialized reaction networks for complex organic molecules (COMs) and prebiotic precursors (Ruaud et al., 2016, Iqbal et al., 2018, Majumdar et al., 2016). Its scientific role is typically that of a chemical solver attached to an external physical model, ranging from static photodissociation-region structures to infalling parcels, shock slabs, disk columns, and tracer-particle histories extracted from magnetohydrodynamic collapse simulations (Gal et al., 2017, Wakelam et al., 2014, Navarro-Almaida et al., 2024).

1. Phase structure and model identity

The foundational distinction in the NAUTILUS literature is between two-phase and three-phase formulations. In the older two-phase picture, the system is divided into gas and grain material, with the adsorbed reservoir treated as a single chemically active phase. The three-phase extension separates that grain reservoir into a grain surface and a grain mantle, so that gas, surface, and mantle are evolved distinctly (Ruaud et al., 2016). In that formulation, the surface is the outer active ice, while the mantle is the deeper bulk ice; both can be chemically active, and the model includes exchange between them.

In the three-phase cold-core implementation, the outermost two monolayers are treated as the surface, while material below is mantle. Surface species can desorb directly, whereas mantle species are ordinarily trapped unless a process acts on the bulk ice or mantle-to-surface replenishment occurs (Wakelam et al., 2021). Other applications use a simpler two-phase gas + grain-surface description. In high-mass star-forming region simulations of the chemical clock $N(\mathrm{HC_3N})/N(\mathrm{N_2H^+})$, a two-phase model was adopted as the main framework because a tested three-phase model gave poorer agreement, with mantle species described there as largely chemically inert at low temperature (Wang et al., 26 Feb 2025).

This variability is methodologically important. It suggests that “NAUTILUS” denotes a common chemical architecture rather than a single standardized setup. The phase partition, degree of mantle activity, and even whether grain populations are single-size or multi-size are selected according to the target problem (Ruaud et al., 2016, Iqbal et al., 2018).

2. Chemical formalism and grain-mediated processes

NAUTILUS is consistently described as a rate-equation model. In static PDR applications, the abundance of species $i$ is evolved with coupled production and loss terms written in the generic form

$\frac{dn_i}{dt} = P_i - L_i ,$

with $P_i$ and $L_i$ including gas-phase reactions, adsorption, desorption, and grain chemistry (Gal et al., 2017). In three-phase cold-core work, separate coupled equations are written for gas, surface, and mantle abundances, with explicit accretion, desorption, dissociation, and surface–mantle transfer terms (Ruaud et al., 2016).

A defining contribution of the three-phase formalism is reaction–diffusion–evaporation competition. For surface or mantle chemistry, the Langmuir–Hinshelwood rate is expressed as

$$k_{ij}^x = \kappa_{ij}\Bigg(\frac{1}{t_\text{hop}^x(i)} + \frac{1}{t_\text{hop}^x(j)}\Bigg)\frac{1}{N_\text{site}n_\text{dust}},$$

with $x=s$ or $m$, and with $\kappa_{ij}$ modified for barriered reactions by competition between barrier crossing, hopping, and evaporation (Ruaud et al., 2016). In this treatment, barriered reactions involving H and even H$_2$ can be strongly enhanced relative to simpler rate prescriptions. One reported consequence is that major N-bearing grain reservoirs shift from NH$i$0 to N$i$1 and HCN when reaction-diffusion competition is included (Ruaud et al., 2016).

Across applications, the recurrent physical processes are gas-phase reactions, adsorption or freeze-out, grain-surface diffusion and reaction, and thermal plus non-thermal desorption. The non-thermal channels commonly included are cosmic-ray-induced desorption, photodesorption, and chemical desorption (Gal et al., 2017, Navarro-Almaida et al., 2024, Wakelam et al., 2016). Several studies further extend this set. Cosmic-ray sputtering has been implemented as an additional first-order loss process from ice phases in Nautilus-1.1, with a sputtering rate coefficient $i$2 added to the total desorption budget (Paulive et al., 2022). In disk chemistry, photoprocess rates have been evaluated from frequency-dependent cross-sections using

$i$3

with shielding handled through a wavelength-dependent local field rather than a single empirical attenuation factor (Gavino et al., 2021).

3. Coupling to physical structure and dynamics

NAUTILUS is usually not the source of the physical model. Instead, it is fed density, temperature, extinction, and radiation fields from external calculations or parameterized histories. In the Horsehead nebula study, the Meudon PDR Code supplies a 1D depth-dependent structure, and NAUTILUS computes time-dependent chemistry independently in each layer of that static background (Gal et al., 2017). In a protoplanetary-disk vertical-column model at 300 au, NAUTILUS is run in 1D with prescribed hydrostatic density, vertical temperature structure, self-shielding, and UV transport (Wakelam et al., 2016).

Other studies attach the model to explicitly evolving trajectories. For IRAS 16293-2422, the chemistry is solved in infalling fluid parcels inside a 1D protostellar core model, allowing the code to follow changing conditions along parcel histories rather than a fixed grid point (Wakelam et al., 2014). In first-hydrostatic-core work, $i$4 tracer particles from a $i$5 non-ideal MHD simulation are postprocessed with NAUTILUS, so that the code receives the evolving density and temperature histories of individual parcels during collapse (Navarro-Almaida et al., 2024).

Shock modeling provides another coupling mode. In C-shock chemistry for L1157, an adapted three-phase NAUTILUS is embedded in a 1-D plane-parallel slab using parametric fits for gas velocity, grain velocity, density, and temperature, plus dynamic dust heating and sputtering (Burkhardt et al., 2019). In high-mass star-forming region models, the code is coupled to 1D spherically symmetric evolutionary structures spanning HMSC, HMPO, HMC, and UCHII stages, with density and temperature represented by modified power laws and physical jumps between stages (Wang et al., 26 Feb 2025).

This range of couplings shows that the model’s core competence is chemical time integration under externally supplied microphysical conditions. A plausible implication is that NAUTILUS serves as a modular chemical back end for environments whose dynamics, radiative transfer, and dust physics are computed elsewhere.

4. Major extensions of the framework

The literature contains several substantive extensions beyond the baseline gas–grain network. One of the earliest broad chemical generalizations is the public network for multiply deuterated species and spin chemistry of H$i$6, H$i$7, H$i$8, and their isotopologues. In that implementation, the starting gas-phase network has 489 species and 7509 reactions, while the extended network reaches about 111,000 gas-phase reactions and about 7700 surface reactions; Oka-based nuclear-spin branching is used for the light hydrogenic system, and the initial ortho:para H$i$9 ratio is taken as 3:1 (Majumdar et al., 2016).

Dust-grain generalization is another major direction. The multi-grain Nautilus model divides the grain population into size bins from $\frac{dn_i}{dt} = P_i - L_i ,$0 to $\frac{dn_i}{dt} = P_i - L_i ,$1, typically using 10, 30, or 60 logarithmically spaced bins and MRN or WD size distributions (Iqbal et al., 2018). Each bin has its own surface chemistry, temperature, and cosmic-ray-induced desorption rate, while all bins are coupled through a shared gas phase. The smallest grains can reach $\frac{dn_i}{dt} = P_i - L_i ,$2 after cosmic-ray hits, whereas $\frac{dn_i}{dt} = P_i - L_i ,$3 grains reach only about $\frac{dn_i}{dt} = P_i - L_i ,$4, making grain-size-dependent desorption a major chemical control (Iqbal et al., 2018). In disks, a separate multi-grain implementation uses 16 grain sizes from $\frac{dn_i}{dt} = P_i - L_i ,$5 nm to $\frac{dn_i}{dt} = P_i - L_i ,$6 mm, with temperatures and UV fluxes computed by POLARIS and passed to the three-phase chemistry (Gavino et al., 2021).

Specialized non-thermal desorption physics has also been added. One study tests chemical desorption, photodesorption, whole-grain cosmic-ray heating, and cosmic-ray sputtering on a 1D TMC1 structure, emphasizing that sputtering can act on both surface and mantle species simultaneously (Wakelam et al., 2021). Another implements sputtering in Nautilus-1.1 as a first-order mantle-loss term, with rates derived from estimated sputtering cross sections and interstellar cosmic-ray fluxes for water ice, CO$\frac{dn_i}{dt} = P_i - L_i ,$7 ice, and a simple mixed ice (Paulive et al., 2022).

Chemical-network specialization is a further hallmark. A dedicated CH$\frac{dn_i}{dt} = P_i - L_i ,$8NCO network was added to the NAUTILUS 3-phase model after quantum-chemical screening of candidate routes, with the conclusion that grain-surface or ice chemistry is the dominant origin of methyl isocyanate (Majumdar et al., 2017). Hot-core glycine modeling added nearly 4500 grain-surface reactions and about 9500 gas-phase reactions, with a revised glycine binding energy of $\frac{dn_i}{dt} = P_i - L_i ,$9 K shifting emphasis toward gas-phase production from precursors rather than simple thermal release (Suzuki et al., 2018). Aminoacetonitrile chemistry added over 300 reactions for AAN and its isomers to a three-phase hot-core model (Zhang et al., 2024). A TMC-1 propargyl study updated the network around resonance-stabilized radicals, added grain-surface hydrogenation routes, and improved the modeled propargyl abundance by about half an order of magnitude relative to an earlier network, though the radical remained underpredicted by two orders of magnitude (Byrne et al., 2023).

5. Representative scientific applications

In cold cores, the three-phase model has been used to compare gas and ice predictions with observations of dark clouds and young stellar object sight lines. The active-mantle, reaction-diffusion-competition formulation reproduces observed main ice species well and leads to a strong late-time decrease in gas abundances, which was interpreted as constraining cold-core chemical ages to a few $P_i$0 yr rather than $P_i$1 yr (Ruaud et al., 2016). In TMC-1-like conditions, the deuteration and spin implementation reproduces many observed D/H ratios reasonably and shows that explicit spin chemistry lowers deuterium fractionation relative to models without it (Majumdar et al., 2016).

In UV-irradiated regions, the Horsehead nebula study uses Nautilus on Meudon PDR structures and concludes that grain-surface chemistry is essential for reproducing methanol and other small organics; FUV photodesorption is singled out as crucial near the illuminated edge (Gal et al., 2017). In disks, the treatment of H$P_i$2 physisorption alone can alter vertically integrated gas-phase column densities of CO, CS, CN, and HCN by factors of several to tens, with the encounter-desorption prescription providing an intermediate behavior between almost no H$P_i$3 sticking and effectively complete trapping of H$P_i$4 on grains (Wakelam et al., 2016). Grain-size-dependent dust temperatures in multi-grain disk chemistry produce broadened snowline behavior, enhanced COM production on warmer grains, and distinct partitioning of C and O into H$P_i$5O, CH$P_i$6, and CO$P_i$7 ices (Gavino et al., 2021).

In protostellar and hot-core environments, NAUTILUS has been used to analyze inheritance and warm-up chemistry. In first-hydrostatic-core collapse postprocessing, most major C, O, N, and S reservoirs in the hot corino are reported to reflect prestellar ice composition, while methyl formate, acetaldehyde, and formamide are formed during the warm-up phase (Navarro-Almaida et al., 2024). For CH$P_i$8NCO, the model indicates an ice-phase origin and thermal release in warm inner regions (Majumdar et al., 2017). In glycine chemistry, the high binding energy adopted from TPD experiments implies substantial grain-surface destruction before full evaporation, making gas-phase formation from precursor chemistry more important than direct thermal desorption (Suzuki et al., 2018). For aminoacetonitrile in hot cores, the three-phase model predicts peak gas-phase abundances of order $P_i$9, consistent with observations toward Sgr B2(N1), and attributes AAN formation primarily to free-radical reactions on grain surfaces during early stages followed by thermal desorption during warm-up (Zhang et al., 2024).

In shocks and dynamical collapse, the framework has been used to move beyond static freeze-out–desorption balances. The adapted C-shock model for L1157 includes sputtering, high-temperature gas-phase reactions, and collisional dust heating, and finds that NH$L_i$0CHO, HCOOCH$L_i$1, and CH$L_i$2CHO exhibit significant post-shock formation pathways rather than tracing sputtering alone (Burkhardt et al., 2019). In high-mass star-forming regions, 1D stage-dependent models reproduce the observed $L_i$3 ratio only when large spatial averages and envelope contributions are included; the best model matches the observed medians at late HMSC, early HMPO, and early UCHII times and identifies 157 observable candidate chemical-clock ratios out of 178 evolutionary trends found among 350 tested ratios (Wang et al., 26 Feb 2025).

6. Limitations, sensitivities, and recurrent methodological issues

The literature repeatedly emphasizes that NAUTILUS results are controlled as much by the supplied physical model and microphysical assumptions as by the chemical solver itself. Static-structure applications inherit uncertainties in density, temperature, and radiation fields from the external physical calculation; the Horsehead study explicitly notes that the chemistry is not dynamically coupled and therefore cannot capture evolving density, turbulent mixing, or grain processing beyond the parameterized chemistry (Gal et al., 2017). In high-mass star-forming region modeling, the authors state that 1D jump models are reaching their limits and that smoother physical evolution is needed (Wang et al., 26 Feb 2025).

Grain physics is a major uncertainty. Single-size descriptions remain common, but multi-grain work shows that small grains dominate cosmic-ray-induced desorption and can create net mass transfer from small to large grains through size-dependent heating (Iqbal et al., 2018). In protoplanetary disks, single-grain models are said to fail to reproduce the complexity introduced by a spread of grain temperatures, surface areas, settling heights, and H$L_i$4-formation efficiencies (Gavino et al., 2021). The H$L_i$5 physisorption problem is especially sensitive: almost-no-sticking and strong-sticking limits were both judged physically unsatisfactory in one disk study, with encounter desorption preferred as an intermediate prescription (Wakelam et al., 2016).

Non-thermal desorption remains another unsettled area. In cold-core TMC1 modeling, chemical desorption is found essential below $L_i$6 cm$L_i$7, while sputtering becomes essential above that density, yet the standard models still underproduce methanol and require consideration of more efficient chemical desorption or more efficient sputtering, especially for CO$L_i$8-rich ices (Wakelam et al., 2021). The sputtering study based on Nautilus-1.1 likewise concludes that quantitative impacts depend strongly on assumed ice composition and that mixed-ice sputtering physics needs improvement (Paulive et al., 2022).

Finally, some NAUTILUS-based conclusions remain explicitly speculative. The aminoacetonitrile study proposes that AAN and NCCN might react toward adenine formation, but states that the efficiency of this route on grain surfaces requires future theoretical or laboratory validation (Zhang et al., 2024). The glycine study suggests that suprathermal H atoms could greatly accelerate COOH formation from CO$L_i$9, but identifies that mechanism as requiring detailed theoretical and experimental investigation (Suzuki et al., 2018). Such cases underline a broader point: the framework is chemically expressive, but its predictions remain contingent on uncertain surface kinetics, binding energies, branching ratios, and desorption physics.