Irregularly Sampled Time Series Interpolation for Binary Evolution Simulations Using Dynamic Time Warping

1. Problem Definition

Binary stellar evolution simulations (e.g., using the MESA code within the POSYDON framework) are computationally expensive, often requiring hundreds of CPU hours per system. To perform large-scale population synthesis, researchers rely on interpolating pre-simulated evolutionary tracks for arbitrary initial conditions. However, traditional interpolation methods face two critical challenges in binary systems:

1. Irregular Sampling: Evolutionary codes use adaptive time-stepping, concentrating resolution on rapid events (like mass transfer) and using sparse sampling during quiescent phases. This results in tracks of varying lengths and non-aligned time steps.
2. Morphological Complexity & Discontinuities: Unlike single stars, binary systems involve mutual interactions (mass transfer, tidal forces) that cause abrupt, non-linear changes in evolutionary trajectories. These interactions create discontinuities that standard interpolation methods cannot align correctly.
3. Physical Inconsistency: Existing methods (e.g., the "Change Point" algorithm) often align parameters independently or reduce tracks to uniform lengths by selecting representative points. This can decouple physically linked parameters (e.g., Luminosity $L$, Radius $R$, and Effective Temperature $T_{eff}$), violating fundamental laws like the Stefan-Boltzmann law ($L \propto R^2 T_{eff}^4$).

2. Methodology

The authors propose a novel framework combining K-Nearest Neighbor (KNN) retrieval with Dynamic Time Warping (DTW) and Iterative Weighted Averaging.

A. Data Representation

• Input: A 3D grid of initial conditions defined by Primary Mass ($M_1$), Secondary Mass ($M_2$), and Orbital Period ($P$).
• Tracks: Each grid point contains a multivariate time series of physical parameters (e.g., $M_1, M_2, P, \dot{M}_{transfer}, T_{eff}, L, R$).
• Classification: Tracks are classified by mass transfer behavior (No Mass Transfer, Stable, Unstable, Initial). The method uses Class-Aware Interpolation, restricting neighbor searches to tracks within the same evolutionary class to avoid averaging physically distinct pathways.

B. K-Nearest Neighbor Retrieval

For a target initial condition $i^*$, the algorithm identifies $K$ nearest neighbors in the grid.

• Weighting: Instead of simple Euclidean distance, the method employs barycentric weights derived from Delaunay triangulation of the 3D parameter space. This provides a more geometrically accurate interpolation weight than simple distance decay.

C. Dynamic Time Warping (DTW) Alignment

To handle irregular sampling and morphological shifts, the method uses DTW to align neighboring tracks before averaging.

• Joint Alignment: A critical innovation is the computation of a single shared warping path across a subset of alignment parameters (e.g., orbital period and luminosity). This path is then applied simultaneously to all physical parameters ($D$ dimensions).
• Preprocessing: To prevent the "age" variable (millions of years) from dominating the distance metric over physical parameters (tens/hundreds), both time and parameter values are normalized to $[0, 1]$.
• Variants: The authors test vanilla DTW, Constrained DTW (CDTW), Derivative DTW (DDTW), and Shape DTW (SDTW), finding that standard DTW with 0-1 normalization and class-aware constraints performs best.

D. Iterative Weighted Averaging

Since DTW aligns two tracks at a time, the authors implement an iterative strategy for $K$ neighbors:

1. Sort neighbors by barycentric weight (most similar first).
2. Align the first two tracks and compute their weighted average.
3. Iteratively align the next neighbor to the current averaged track and update the average.
4. This process preserves the temporal consistency of the final interpolated track.

3. Key Contributions

1. Joint-Alignment Strategy: The paper introduces a method to align all physical parameters of a binary system using a single temporal grid. This ensures that coupled quantities (like $L, R, T_{eff}$) remain causally consistent, preventing physical anomalies.
2. Theoretical Proof of Physical Consistency: The authors prove mathematically that if neighboring tracks satisfy the Stefan-Boltzmann law and parameters are in logarithmic scale, the weighted average of aligned tracks preserves this law exactly. This is a direct consequence of the linearity of averaging in log-space combined with shared temporal alignment.
3. Class-Aware Interpolation: By segregating the grid into evolutionary classes (NMT, SMT, UMT) before neighbor selection, the method avoids averaging tracks with fundamentally different morphologies.
4. Comprehensive Evaluation: The study introduces Area Error as a metric to evaluate the total deviation between curves, addressing the limitations of point-wise error metrics (RAE) which fail to capture morphological fidelity in irregularly sampled data.

4. Results

The method was evaluated on three POSYDON grids: HMS-HMS, CO-HeMS, and CO-HMS.

• Performance vs. Baselines: The DTW-based method consistently outperformed the "Change Point" algorithm and simple Nearest Neighbor (NN) baselines.
  • Median Relative Absolute Error (RAE): Reduced significantly across most parameters. For example, in the HMS-HMS grid, the median error for Mass Transfer Rate ($\dot{M}$) dropped from ~22% (Change Point) to ~6% (DTW).
  • Area Error: The DTW method achieved the lowest area errors, indicating superior preservation of the overall track shape and evolutionary history.
• Visual Fidelity: Visual comparisons showed that the DTW method successfully captured sharp spikes and rapid transitions (e.g., mass transfer onset) that the Change Point algorithm smoothed over or misaligned.
• Physical Consistency:
  • HR Diagram: Interpolated tracks maintained coherent trajectories in the Hertzsprung-Russell diagram.
  • Stefan-Boltzmann Law: While the Change Point algorithm produced spurious spikes in the Stefan-Boltzmann ratio ($L / (R^2 T^4)$) due to temporal misalignment of $L, R,$ and $T$, the DTW method maintained the ratio close to unity throughout the evolution.

5. Significance and Impact

• Computational Efficiency: The method reduces the computational cost of generating binary population samples by six orders of magnitude. While full MESA simulations take ~10–100 CPU hours per system, the interpolation method generates a full track in ~0.01 seconds.
• Astrophysical Utility: It enables the generation of large-scale, statistically robust binary populations for gravitational wave astronomy and stellar population studies, which were previously restricted by computational bottlenecks.
• Generalizability: While applied to stellar evolution, the framework of geometric neighbor selection combined with sequence alignment (DTW) offers a general solution for interpolating irregularly sampled time series in other scientific domains (e.g., climate science, fluid dynamics) where morphological preservation is critical.

In conclusion, this work bridges the gap between high-fidelity physics simulations and large-scale statistical analysis by solving the fundamental problem of aligning irregular, interaction-driven time series while rigorously preserving physical laws.