Triple Stellar Evolution Code

Updated 6 December 2025

Triple Stellar Evolution (TSE) code is a modular framework that simulates the coupled stellar and dynamical evolution of hierarchical triple star systems.
It integrates orbit-averaged secular dynamics, direct N-body integration, and detailed prescriptions for processes like mass transfer, common-envelope phases, and gravitational radiation.
The code employs an event-driven architecture with SSE/BSE routines, enabling efficient statistical studies of triple-induced mergers and gravitational wave sources.

The Triple Stellar Evolution (TSE) code is a modular population-synthesis framework for simulating the coupled stellar and dynamical evolution of hierarchical triple star systems. TSE implements orbit-averaged secular dynamics, direct N-body integration for non-hierarchical or chaotic regimes, detailed prescriptions for mass transfer and common-envelope phases, and full integration with standard stellar evolution recipes. This enables modeling of a diverse range of evolutionary pathways, including the formation of compact object binaries and their potential as sources of gravitational waves, over large synthetic populations with comprehensive event tracking and output diagnostics (Hamers et al., 2020, Stegmann et al., 2021).

1. Software Architecture and Data Structures

TSE is structured as an event-driven, modular code, with a core written in C++ (or C/Fortran style in some variants) and interfaces in Python/Cython for setup and analysis. The software links to:

SSE/BSE routines: for fast single- and binary-star evolution tracks (mass, radius, metallicity, core/envelope state).
Orbit-averaged secular integrator: double-averaged Hamiltonian dynamics up to octupole order, encapsulating von Zeipel–Lidov–Kozai (vZLK) cycles, with additional 1PN, 2PN, and 2.5PN corrections.
Direct N-body engine: for dynamical instabilities, using regularized N-body codes such as MSTAR.
Interaction modules: Roche lobe overflow (RLOF), common envelope (CE), tidal friction, mass transfer, supernova/natal kicks.

The triple system is represented as a recursive hierarchy of “Particle” objects, each holding physical parameters and pointers to child components (stars or binaries). Typical workflow proceeds by time-stepping all modules synchronously, coordinated by a driver process that selects the shortest relevant physical timescale at each step and processes corresponding events (Hamers et al., 2020, Perets, 3 Apr 2025).

2. Dynamical Regimes and Secular Evolution

The core of triple dynamics in TSE is the secular approximation, valid for strongly hierarchical systems (a_out ≫ a_in), where the total Hamiltonian is expanded in powers of the small parameter α = a_in/a_out:

Quadrupole term: governs classical vZLK oscillations, leading to periodic exchanges of inner eccentricity and inclination.
Octupole term: arises when mass asymmetry (m1 ≠ m2) and outer eccentricity are nonzero, introducing extreme eccentricity excursions and potential orbital flips.

Explicitly, the octupole strength parameter is

$\epsilon_{\rm oct} = \frac{m_1 - m_2}{m_1 + m_2}\; \frac{a_{\rm in}}{a_{\rm out}}\; \frac{e_{\rm out}}{1-e_{\rm out}^2}$

Secular equations of motion are evolved for orbital angular momentum and eccentricity vectors:

$\frac{d\mathbf{e}_k}{dt},\qquad \frac{d\mathbf{h}_k}{dt}$

with timescales (e.g., for quadrupole vZLK)

$\tau_{\rm ZLK} \sim \frac{P_{\rm out}^2}{P_{\rm in}} \frac{m_1+m_2+m_3}{m_3} (1-e_{\rm out}^2)^{3/2}$

Breakdown of the secular approximation (e.g., in the semisecular or chaotic regime) is detected by comparing angular-momentum change timescales or by violation of the hierarchical stability criterion, triggering a switch to direct N-body integration (Hamers et al., 2020, Antonini et al., 2017, Perets, 3 Apr 2025).

3. Physical Processes: Stellar and Binary Evolution

TSE incorporates detailed prescriptions for all relevant physical processes in triple evolution:

Stellar evolution: via analytic SSE/BSE fitting formulae, including metallicity-dependent winds, remnant formation, and evolutionary phase transitions.
Tidal processes: equilibrium tide models (Hut 1981) for both convective and radiative envelopes; coupled spin evolution, circularization, and semi-major axis evolution.
Mass transfer and common envelope:
- Stable mass transfer: RLOF detection using the Eggleton approximation; mass exchange rate and angular momentum loss via analytic or polytropic-stream models, including Hamers & Dosopoulou (2019) for eccentric orbits.
- Common envelope: α–λ energy formalism; post-CE separations and merger criteria based on envelope binding energy and orbital energy change.
Gravitational radiation: Peters (1964) prescription for orbital decay and eccentricity damping, relevant for compact object binaries.
Supernovae and natal kicks: instantaneous mass loss and velocity kick implementation, with resulting recalculation of orbital elements, possible system unbinding, and fallback-modulated kick magnitudes for black holes.

These modules interact self-consistently, allowing the tertiary to induce RLOF or merger in otherwise noninteracting binaries via eccentricity excitation, and permitting tertiary RLOF and dynamically unstable episodes with recourse to N-body integration (Hamers et al., 2020, Rajamuthukumar et al., 13 Feb 2025, Stegmann et al., 2021).

4. Stability Criteria and Dynamical Transitions

The dynamical regime of a triple is governed by the Mardling & Aarseth (2001) criterion, which stipulates that

$\frac{a_{\rm out} (1-e_{\rm out})}{a_{\rm in}} > 2.8 \left[ (1+q_{\rm out}) \frac{1+e_{\rm out}}{\sqrt{1-e_{\rm out}}} \right]^{2/5} \left[1-0.3\,\frac{i}{\pi}\right]$

If violated, the system is flagged as dynamically unstable and the code transitions from secular evolution to direct N-body integration. Additional triggers include the entrance into the semisecular regime (where the angular momentum changes on timescales shorter than the orbital period) and unbinding events.

While in the N-body regime, the code periodically checks whether a stable hierarchy has re-formed and, if so, reconstructs orbital elements and returns to secular evolution provided the fractional change in orbital elements over a test interval falls below threshold (∼1%) (Hamers et al., 2020, Antonini et al., 2017).

5. Input Parameters, Boundary Conditions, and Procedures

TSE adopts initial conditions based on observed properties of massive stars and multiples, e.g.,

Masses from Kroupa or Salpeter IMF, with secondary and tertiary masses from empirical or uniformly sampled mass-ratio distributions.
Initial semi-major axes and eccentricities from log-flat or thermal distributions; mutual inclinations drawn isotropically or following observed correlations.
Metallicity, wind prescription, and natal kick models specified by user or sampled across population.
Systems violating the stability criterion or with initial Roche-lobe overflow are rejected at initialization.

Simulation proceeds via timesteps determined by the minimum of the relevant physics timescales (secular, tidal, gravitational wave, mass transfer, nuclear evolution), and proceeds until a stopping condition is met (e.g., merger, unbinding, core-collapse, or Hubble time) (Rajamuthukumar et al., 13 Feb 2025, Stegmann et al., 2021).

Example pseudocode capturing the population-synthesis workflow is provided in (Rajamuthukumar et al., 13 Feb 2025, Perets, 3 Apr 2025):

for i in 1..N_triples:
    sample ICs
    while t < t_max:
        compute dt = min(ts, t_LK, t_tide, t_GW, t_MT)
        advance all modules by dt
        process discrete events (RLOF, CE, SN, instability)
        if system unbound or merger: break

6. Outputs, Diagnostics, and Validation

TSE tracks and outputs, for each system:

Time series of all relevant orbital and stellar parameters (a, e, i, mass, stellar type, spin, etc.) for each component.
Event log: records RLOF onset, mass transfer onset, common envelope, SN kicks, mergers, collisions, secular/N-body switching.
Evolutionary flags indicating channel (e.g., inner mass transfer, tertiary mass transfer, dynamical instability, compact object formation).
Terminal configuration: binary or triple DCOs, unbound remnants, merged objects.

Validation is performed by comparing against analytic secular test cases (e.g., vZLK cycles, maximum eccentricity, timescales) and, in the N-body and merger regime, by comparison with direct few-body simulations and published benchmark results. Performance benchmarks indicate high efficiency for large populations (∼10^5–10⁶ systems), suitable for statistical studies of triple-induced mergers and transients (Hamers et al., 2020, Stegmann et al., 2021, Perets, 3 Apr 2025).

7. Limitations and Prospects

Present limitations and approximations:

Secular dynamics truncated at octupole order; higher-order resonance and cross-terms neglected except in N-body regime.
Tidal evolution is based on equilibrium models; dynamical tides and complex dissipation physics are excluded.
Eccentric mass transfer employs analytic recipes; full hydrodynamic modeling is not included.
SN kick prescriptions remain heuristic with large empirical uncertainty.
Planet and fly-by modeling is restricted to simple analytic approximations; long-lived external perturbations are not tracked.
Sensitive dependence on initial conditions and numerical rounding can result in run-to-run variation at the single-system level; the code is best suited for statistical population studies.

Possible future improvements include implementation of higher-order multipole dynamics, more sophisticated treatments of tides, magnetic braking, and integration with recent empirical distributions for high-mass triples (Hamers et al., 2020, Stegmann et al., 2021).

References:

(Hamers et al., 2020, Rajamuthukumar et al., 13 Feb 2025, Toonen et al., 2016, Kummer et al., 2023, Antonini et al., 2017, Perets, 3 Apr 2025, Stegmann et al., 2021)