GalCEM: Galactic Chemical Evolution Model

• GalCEM is an open-source, Python-based framework designed for precision isotopic chemical evolution modeling in galaxies, tracking 86 elements and 451 isotopes.
• It employs a modular architecture that integrates detailed yield data from massive stars, LIMs, and SNIa using advanced numerical methods like a 4th-order Runge–Kutta integrator.
• Its design enables fast, high-fidelity simulations validated against observational benchmarks, with extensible features for future multi-zone and neutron-capture enhancements.

GalCEM (GALactic Chemical Evolution Model) is an open-source, modular, Python-based numerical framework for precision isotopic chemical evolution modeling within galaxies. Designed to solve the full integrodifferential one-zone chemical evolution equations governing the time-dependent masses of individual isotopes in the interstellar medium (ISM), GalCEM integrates detailed yield data from multiple nucleosynthetic channels, leveraging array-based numerical engines to achieve both high speed and high fidelity. The current implementation tracks 86 elements subdivided into 451 isotopes and includes contributions from massive stars, low- and intermediate-mass stars, and Type Ia supernovae. GalCEM is publicly available and extensible, providing researchers with a flexible platform for simulating the chemical enrichment history from light to intermediate elements, with future expansions planned for neutron-capture channels and spatially resolved treatments (Gjergo et al., 2023).

1. Mathematical Formulation and Core Equations

GalCEM solves the canonical one-zone galactic chemical evolution (GCE) equation for each tracked isotope $i$:

$$\frac{\mathrm{d} M_{i,\rm gas}}{\mathrm{d} t} = X_{i,\rm inf} \, \dot{M}_{\rm inf}(t) - (1-\omega) \, X_{i,\rm gas}(t) \, \psi(t) + \sum_P R_{P,i}(t)$$

where:

• $M_{i,\rm gas}(t)$ is the mass of isotope $i$ in the ISM,
• $X_{i,\rm inf}$ and $X_{i,\rm gas}$ are the abundance fractions in infalling and ISM gas,
• $\psi(t)$ is the star formation rate (SFR),
• $\omega=0$ for no outflow,
• $R_{P,i}(t)$ are the rates of isotope $i$ returned via nucleosynthetic channel $$P \in \{\mathrm{LIM},\,\mathrm{SNCC},\,\mathrm{SNIa}\}$$: low-intermediate mass stars, core-collapse supernovae, and Type Ia supernovae respectively.

For SNCC and LIM channels, the birth-time convolution rate is:

$$R_{P,i}(t) = \int_{t_{\rm min}}^{t_{\rm max}} \psi(t') \left[-\frac{\mathrm{d}M_*}{\mathrm{d}\tau} \phi(M_*) Y_{P,i}(M_*,Z(t')) \right]_{M_*=M_*(t-t', Z(t'))} \;\mathrm{d}t'$$

For SNIa, a delay-time distribution (DTD, Greggio 2005) produces:

$$R_{\rm SNIa}(t) = \alpha_{\rm SNIa}\,Y_{\rm SNIa,i}\,\int_{\tau_{\rm i}}^{\tau_{\rm x}} \psi(t-\tau) \mathrm{DTD}_{\rm SNIa}(\tau) \,\mathrm{d}\tau$$

where $Y_{P,i}$ is the stellar yield, $\phi(M)$ the IMF, $\tau(M,Z)$ the lifetime, $\mathrm{d}M_*/\mathrm{d}\tau$ the Jacobian, and $\mathrm{DTD}_{\rm SNIa}$ is computed analytically from binary mass-ratio distributions. The system is advanced using a classical 4th-order Runge–Kutta method with uniform time step $\Delta t = 2$ Myr.

2. Modular Architecture and Numerical Engines

GalCEM features three primary ingredient classes—morphology, yields, and input parameters—tightly coupled with its numerical engines. Morphology encompasses the IMF ($\phi(M)$), SFR ($\psi(t)$), infall, outflow, and stellar lifetime functions ($\tau(M, Z)$). Yield tables for different enrichment channels (FRUITY AGB [Cristallo et al.], Limongi & Chieffi SNCC, Iwamoto SNIa) provide data for 86 elements and 451 isotopes. Input parameters configure solver options, including channel selection, yield choices, and timestep control.

Numerical engines perform:

1. Multidimensional yield table preprocessing and spline interpolation;
2. Construction of time-step grid arrays for each channel and isotope;
3. Efficient convolution integral evaluation via Simpson’s rule;
4. ODE integration via a fourth-order Runge–Kutta solver.

The entry point, the OneZone class, orchestrates initialization, isotope extraction, and the iterative evolution loop.

3. Yield Table Preprocessing and Interpolation

To optimize computation, GalCEM employs the standalone GalCemInterpolant tool, which ingests raw stellar yield tables and constructs smooth multidimensional splines for each isotope.

• For LIM (FRUITY), interpolation occurs in $(M, Z)$ space.
• For SNCC yields (Limongi & Chieffi), $(M, Z, \mathrm{IRV})$ space is used.
• SNIa yields (Iwamoto) are similarly treated.

Splines are typically constructed in $(\log M, \sqrt Z)$, and the mass–lifetime $(\tau, M)$ relation and its inverse are also interpolated for analytic Jacobian computation. The framework is dimensionally agnostic and supports extensions (e.g., extra dimensions for He enhancement). In-sample RMSE errors are $<10^{-15}$, yielding machine-precision fidelity in tabulated values.

4. Convolution Grid-Mapping and Computational Performance

At each time step $t_n$ and for each enrichment channel $P$, GalCEM constructs a grid in $\log M$ (with $N_k=200$ points) over relevant mass ranges. Array-based mappings provide immediate evaluation of $\tau(M, Z)$, birth times $t'$, IMF values $\phi(M)$, yields $Y_{P,i}(M, Z(t'))$, and stored past SFR values $\psi(t')$. The convolution integral for each isotope/channel reduces to a dot-product of pre-computed arrays, numerically integrated via Simpson’s rule.

This array-centric approach suppresses expensive Python-level looping, providing order-of-magnitude speed enhancements. A single full run (451 isotopes, 86 elements, typical spiral galaxy) completes in minutes on standard workstation hardware (Gjergo et al., 2023).

5. Validation Against Observational Data

Benchmark simulations using Milky Way-like parameters ($M_{\rm tot} = 5\times 10^{10}\,M_\odot$, $\tau_{\rm inf}=7$ Gyr, SFR $\propto M_{\rm gas}$, Kroupa (2001) IMF, standard yield choices) yield age–metallicity and [Fe/H] tracks that intersect the solar point at $t=4.6$ Gyr and agree with APOKASC stellar-age measurements (Silva Aguirre et al. 2018). Elemental [X/Fe] versus [Fe/H] tracks for C to Zn and Li reproduce key features: the $\alpha$-plateau at low [Fe/H], knee at $-1\lesssim \mathrm{[Fe/H]} \lesssim -0.5$, and iron-peak rises from SNIa. Minor discrepancies at the few tenths dex level suggest areas for further model refinement, including treatment of rotating massive stars and additional enrichment sources.

6. Public Availability and Usage

GalCEM is hosted at https://github.com/egjergo/GalCEM and distributed via PyPI:

pip install galcem

A minimal script for one-zone GCE simulation is:

import galcem as gc inputs = gc.Inputs() oz = gc.OneZone(inputs, outdir='my_output') oz.main()

Channels can be toggled and yield options specified:

inputs.include_channel=['SNCC','LIMs','SNIa'] inputs.yields_LIMs_option='c15' inputs.yields_SNCC_option='lc18' inputs.yields_SNIa_option='i99'

Solver timestep granularity is adjustable:

inputs.ntime_step = 0.002 # in Gyr

Complete configuration documentation is provided within galcem/classes/inputs.py. Planned future expansions include adaptive timestepping, multi-zone capability, neutron-capture channels (e.g., r-process), and expanded yield libraries.

7. Roadmap and Prospective Extensions

Current limitations restrict nucleosynthetic coverage to light and intermediate elements (C to Zn, Li) and one-zone galactic models. Upcoming releases aim to incorporate the full periodic table by adding neutron-capture channels and spatially-dependent galaxy properties. Further anticipated features include adaptive timestep solvers and multi-zone enrichment coupling. The GalCEM architecture is designed for extensibility, facilitating precision studies of galactic chemical evolution across parameter regimes and astrophysical contexts (Gjergo et al., 2023).