diff --git a/dingo/asimov/dingo.ini b/dingo/asimov/dingo.ini
index 9672df89d..ea5d0eb5e 100644
--- a/dingo/asimov/dingo.ini
+++ b/dingo/asimov/dingo.ini
@@ -203,9 +203,9 @@ prior-dict-updates = {
 {%- assign p = priors['eccentricity'] %}
    eccentricity = {{ p['type'] }}(name='eccentricity', minimum={{ p['minimum'] }}, maximum={{ p['maximum'] }}),
 {%- endif %}
-{%- if priors.keys() contains "mean_anomaly" %}
-{%- assign p = priors['mean_anomaly'] %}
-   mean_anomaly = {{ p['type'] }}(name='mean_anomaly', minimum={{ p['minimum'] }}, maximum={{ p['maximum'] }}, boundary={{ p['boundary'] | default: 'periodic' }} ),
+{%- if priors.keys() contains "mean_per_ano" %}
+{%- assign p = priors['mean_per_ano'] %}
+   mean_per_ano = {{ p['type'] }}(name='mean_per_ano', minimum={{ p['minimum'] }}, maximum={{ p['maximum'] }}, boundary={{ p['boundary'] | default: 'periodic' }} ),
 {%- endif %}
 {%- if priors.keys() contains "luminosity distance" %}
 {%- assign p = priors['luminosity distance'] %}
diff --git a/dingo/gw/likelihood.py b/dingo/gw/likelihood.py
index 49b29941f..f15be844d 100644
--- a/dingo/gw/likelihood.py
+++ b/dingo/gw/likelihood.py
@@ -14,7 +14,7 @@
     DecimateWaveformsAndASDS,
     create_mask_based_on_frequency_update,
 )
-from dingo.gw.waveform_generator import WaveformGenerator
+from dingo.gw.waveform_generator import WaveformGenerator, wfg_utils
 from dingo.gw.domains import (
     Domain,
     UniformFrequencyDomain,
@@ -382,6 +382,18 @@ def _log_likelihood_phase_grid_mode_decomposed(self, theta, phases=None):
         pol_m = self.signal_m({**theta, "phase": 0})
         pol_m = {k: pol["waveform"] for k, pol in pol_m.items()}
 
+        # For TEOB, generate_hplus_hcross_m applies a phase-dependent time shift
+        # bookkept in wfg._deferred_timeshift_data. The analytical precomputation
+        # below assumes mu(phi) = sum_m mu_m(0) * exp(-i m phi), which is only
+        # valid when that timeshift is absent. When it's present we fall back to
+        # a per-phase resum so the phase-grid likelihood stays consistent with
+        # the full direct likelihood.
+        dts = getattr(self.waveform_generator, "_deferred_timeshift_data", None)
+        if dts is not None:
+            return self._log_likelihood_phase_grid_mode_decomposed_with_dts(
+                pol_m, d, phases, dts
+            )
+
         # Step 2: Precompute complex inner products (mu, mu) and (d, mu) for the
         # individual modes m.
         min_idx = self.data_domain.min_idx
@@ -475,6 +487,56 @@ def _log_likelihood_phase_grid_mode_decomposed(self, theta, phases=None):
 
         return log_likelihoods
 
+    def _log_likelihood_phase_grid_mode_decomposed_with_dts(
+        self, pol_m, d, phases, dts
+    ):
+        """Variant of _log_likelihood_phase_grid_mode_decomposed for approximants
+        (currently TEOB) whose mode decomposition carries a phase-dependent time
+        shift recorded in `dts` (= wfg._deferred_timeshift_data).
+
+        We cannot use the analytical precomputation trick because the timeshift
+        depends on the argmax of the TD-resummed h+, which varies with phase.
+        Instead we resum modes per phase and apply the exp(-2 pi i f dt(phi))
+        correction before computing inner products. The geocenter-level time
+        shift commutes with detector projection, so we can apply it directly to
+        the per-detector signals returned by signal_m.
+        """
+        m_vals = sorted(pol_m.keys())
+        ifos = list(d.keys())
+        freqs = self.data_domain()
+
+        log_likelihoods = np.empty(len(phases))
+        for idx, phi in enumerate(phases):
+            # Analytical resum of per-detector mode contributions.
+            mu = {
+                ifo: sum(
+                    pol_m[m][ifo] * np.exp(-1j * m * phi) for m in m_vals
+                )
+                for ifo in ifos
+            }
+            # Deferred-timeshift correction (no-op when phi == 0.0).
+            if phi != 0.0:
+                target_phase = dts["phase_ref"] + phi
+                epoch_target = wfg_utils.compute_epoch_from_resized_td_modes(
+                    dts["resized_td_modes"], dts["iota"], target_phase,
+                    dts["delta_t"],
+                )
+                dt_target = 1.0 / dts["delta_f"] + epoch_target
+                dt_correction = dt_target - dts["dt_ref"]
+                shift = np.exp(-2j * np.pi * dt_correction * freqs)
+                for ifo in ifos:
+                    mu[ifo] = mu[ifo] * shift
+
+            rho2opt = sum(
+                inner_product(mu_ifo, mu_ifo) for mu_ifo in mu.values()
+            )
+            kappa2 = sum(
+                inner_product(d_ifo, mu_ifo)
+                for d_ifo, mu_ifo in zip(d.values(), mu.values())
+            )
+            log_likelihoods[idx] = self.log_Zn + kappa2 - 0.5 * rho2opt
+        return log_likelihoods
+
     def _log_likelihood_phase_marginalized(self, theta):
         """
 
diff --git a/dingo/gw/result.py b/dingo/gw/result.py
index 0dc3b86c9..9a6d658ef 100644
--- a/dingo/gw/result.py
+++ b/dingo/gw/result.py
@@ -701,6 +701,10 @@ def get_pesummary_samples(self, num_processes=1):
             * The spin angles phi_jl and theta_jn are transformed to account for a
             difference in phase definition.
             * Some columns are dropped: delta_log_prob_target, log_prob
+            * Aligned-spin columns ``chi_1``/``chi_2`` are renamed to
+            ``spin_1z``/``spin_2z`` so that PESummary's standard-name
+            dictionary recognises them and derives ``a_1``, ``a_2``,
+            ``chi_eff``, etc.
         """
         if hasattr(self, "_pesummary_samples"):
             return self._pesummary_samples
@@ -727,18 +731,27 @@ def get_pesummary_samples(self, num_processes=1):
             if "time" in col:
                 samples.loc[:, col] += self.t_ref
 
+        # Rename aligned-spin columns to match PESummary's standard names.
+        # PESummary's standard_names.py maps spin_1z/spin_2z but not chi_1/chi_2,
+        # so without this rename pesummary cannot derive a_1, chi_eff, chi_p, etc.
+        samples.rename(
+            columns={"chi_1": "spin_1z", "chi_2": "spin_2z"}, inplace=True
+        )
+
         spin_conversion_phase_old = self.base_metadata["dataset_settings"][
             "waveform_generator"
         ].get("spin_conversion_phase")
 
-        # Redefine phase parameter to be consistent with Bilby. COMMENTED BECAUSE SLOW
-        samples = change_spin_conversion_phase(
-            samples,
-            self.f_ref,
-            spin_conversion_phase_old,
-            None,
-            num_processes=num_processes,
-        )
+        # Redefine phase parameter to be consistent with Bilby.
+        # Only apply spin conversion when precessing spin parameters are present.
+        if "phi_jl" in samples.columns:
+            samples = change_spin_conversion_phase(
+                samples,
+                self.f_ref,
+                spin_conversion_phase_old,
+                None,
+                num_processes=num_processes,
+            )
 
         self._pesummary_samples = samples
 
@@ -750,7 +763,9 @@ def pesummary_prior(self):
 
         By convention, Dingo stores all times *relative* to a reference time, typically
         the trigger time for an event. The prior returned here corrects for that offset to
-        be consistent with other codes.
+        be consistent with other codes. Aligned-spin keys ``chi_1``/``chi_2`` are
+        renamed to ``spin_1z``/``spin_2z`` to match the renaming applied in
+        :meth:`get_pesummary_samples`.
         """
         prior = copy.deepcopy(self.prior)
         for p in prior:
@@ -760,4 +775,7 @@ def pesummary_prior(self):
                     prior[p].minimum += self.t_ref
                 except AttributeError:
                     continue
+        for old, new in (("chi_1", "spin_1z"), ("chi_2", "spin_2z")):
+            if old in prior and new not in prior:
+                prior[new] = prior.pop(old)
         return prior
diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py
index 43f6a6488..2a17aeca8 100644
--- a/dingo/gw/waveform_generator/waveform_generator.py
+++ b/dingo/gw/waveform_generator/waveform_generator.py
@@ -3,7 +3,9 @@
 from math import isclose
 
 import numpy as np
+import astropy.constants as ac
 import astropy.units as u
+import astropy.constants as ac
 from typing import Dict, List, Tuple, Union, Callable
 from numbers import Number
 import warnings
@@ -235,8 +237,12 @@ def generate_hplus_hcross(
             if not catch_waveform_errors:
                 raise
             else:
-                EDOM = e.args[0] == "Internal function call failed: Input domain error"
-                if EDOM:
+                EDOM = (
+                    len(e.args) > 0
+                    and e.args[0] == "Internal function call failed: Input domain error"
+                )
+                recoverable = EDOM or isinstance(e, RuntimeError)
+                if recoverable:
                     warnings.warn(
                         f"Evaluating the waveform failed with error: {e}\n"
                         f"The parameters were {parameters_generator}\n"
@@ -727,7 +733,7 @@ def generate_hplus_hcross_m(
             if LS.SimInspiralImplementedFDApproximants(self.approximant):
                 # Step 1: generate waveform modes in L0 frame in native domain of
                 # approximant (here: FD)
-                hlm_fd, iota = self.generate_FD_modes_LO(parameters)
+                hlm_fd, iota = self.generate_FD_modes_L0(parameters)
 
                 # Step 2: Transform modes to target domain.
                 # Not required here, as approximant domain and target domain are both FD.
@@ -769,7 +775,7 @@ def generate_hplus_hcross_m(
 
                 self._use_base_domain = True
                 self._domain_transform = DecimateAll(self._domain)
-                hlm_fd, iota = self.generate_FD_modes_LO(parameters)
+                hlm_fd, iota = self.generate_FD_modes_L0(parameters)
                 pol_m = wfg_utils.get_polarizations_from_fd_modes_m(
                     hlm_fd, iota, parameters["phase"]
                 )
@@ -800,7 +806,7 @@ def generate_hplus_hcross_m(
         else:
             return pol_m
 
-    def generate_FD_modes_LO(self, parameters):
+    def generate_FD_modes_L0(self, parameters):
         """
         Generate FD modes in the L0 frame.
 
@@ -957,7 +963,8 @@ def __init__(self, **kwargs):
             "antisymmetric_modes_hm",
         }
         extra_wf_kwargs = {k: v for k, v in kwargs.items() if k in allowed_extra_kwargs}
-        extra_wf_kwargs["lmax_nyquist"] = kwargs.get("lmax_nyquist", 2)
+        if "SEOB" in self.approximant_str and not "ROM" in self.approximant_str: 
+            extra_wf_kwargs["lmax_nyquist"] = kwargs.get("lmax_nyquist", 2)
         self.extra_wf_kwargs = extra_wf_kwargs
 
     @property
@@ -1028,6 +1035,10 @@ def _convert_parameters(
             f_min = self.f_start
         else:
             f_min = self.domain.f_min
+            # for SEOBNRv5EHM, the starting frequency must be the same as the reference frequency
+            # note this is the orbit averaged reference frequency
+            if self.approximant_str == "SEOBNRv5EHM" or self.approximant_str == "TEOBResumSDALI":
+                f_min = self.f_ref
         # parameters needed for TD waveforms
         delta_t = 0.5 / self.domain.f_max
 
@@ -1052,8 +1063,39 @@ def _convert_parameters(
             "condition": 1,
         }
 
-        # SEOBNRv5 specific parameters
-        params_gwsignal.update(self.extra_wf_kwargs)
+        if self.approximant_str in ["SEOBNRv5EHM", "TEOBResumSDALI"]:
+            # eccentric parameters
+            log_ecc = p.get("log10_eccentricity")
+            ecc = p.get("eccentricity")
+
+            if log_ecc is not None and ecc is not None:
+                if abs(10**log_ecc - ecc) > 1e-6:
+                    raise ValueError(
+                        f"log10_eccentricity={log_ecc} and eccentricity={ecc} are inconsistent."
+                    )
+
+            eccentricity = 10**log_ecc if log_ecc is not None else p.get("eccentricity", 0.0)
+            longitude_ascending_nodes = p.get("long_asc_nodes", 0.0)
+            if "relativistic_anomaly" in p and "mean_per_ano" in p:
+                raise ValueError(
+                    "Cannot specify both relativistic_anomaly and mean_per_ano."
+                )
+            else:
+                mean_per_ano = p.get("mean_per_ano", p.get("relativistic_anomaly", 0.0))
+
+            params_gwsignal.update({
+                'eccentricity' : eccentricity * u.dimensionless_unscaled,
+                'longAscNodes' : longitude_ascending_nodes * u.rad,
+                'meanPerAno' : mean_per_ano * u.rad,
+            })
+
+        # TEOBResumSDALI doesn't support lmax_nyquist, so filter it out
+        if "TEOBResum" in self.approximant_str:
+            extra_wf_kwargs_filtered = {k: v for k, v in self.extra_wf_kwargs.items()
+                                       if k != "lmax_nyquist"}
+            params_gwsignal.update(extra_wf_kwargs_filtered)
+        else:
+            params_gwsignal.update(self.extra_wf_kwargs)
 
         if return_target_function:
             # This is a hack to make compatible with LAL version. Target functions for
@@ -1125,14 +1167,22 @@ def generate_FD_waveform(self, parameters_gwsignal: Dict) -> Dict[str, np.ndarra
         frequency_array = self.domain()
         h_plus = np.zeros_like(frequency_array, dtype=complex)
         h_cross = np.zeros_like(frequency_array, dtype=complex)
-        # Ensure that length of wf agrees with length of domain. Enforce by truncating frequencies beyond f_max
+        # Ensure that length of wf agrees with length of domain.
         if len(hp) > len(frequency_array):
             warnings.warn(
                 "GWSignal waveform longer than domain's `frequency_array`"
                 f"({len(hp)} vs {len(frequency_array)}). Truncating gwsignal array."
             )
-            h_plus = hp[: len(h_plus)].value
-            h_cross = hc[: len(h_cross)].value
+            h_plus[:] = hp[: len(h_plus)].value
+            h_cross[:] = hc[: len(h_cross)].value
+        elif len(hp) < len(frequency_array):
+            # For TD models (e.g. TEOBResumS), GenerateFDWaveform chooses the TD
+            # sampling rate based on the waveform's frequency content, not the
+            # requested f_max. High-mass systems merge at lower frequencies, so
+            # the FFT can have fewer bins than the domain. Zero-padding is correct
+            # since there is no signal power above that frequency.
+            h_plus[: len(hp)] = hp.value
+            h_cross[: len(hc)] = hc.value
         else:
             h_plus = hp.value
             h_cross = hc.value
@@ -1212,19 +1262,25 @@ def generate_hplus_hcross_m(
 
         generator = new_interface_get_waveform_generator(self.approximant_str)
         if isinstance(self.domain, UniformFrequencyDomain):
-            # Generate FD modes in for frequencies [-f_max, ..., 0, ..., f_max].
-            if generator.domain == "freq":
+            # Generate FD modes for frequencies [-f_max, ..., 0, ..., f_max].
+            if self.approximant_str in ["SEOBNRv5EHM"]:
+                # SEOBNRv5EHM and TEOBResumSDALI use these conditioning routines
+                # assert LS.SimInspiralImplementedTDApproximants(self.approximant)
                 # Step 1: generate waveform modes in L0 frame in native domain of
-                # approximant (here: FD)
-                hlm_fd, iota = self.generate_FD_modes_LO(parameters)
+                # approximant (here: TD)
+                hlm_td, iota = self.generate_TD_modes_L0(parameters)
 
                 # Step 2: Transform modes to target domain.
-                # Not required here, as approximant domain and target domain are both FD.
+                # This requires tapering of TD modes, and FFT to transform to FD.
+                wfg_utils.taper_td_modes_in_place(hlm_td)
+                hlm_fd = wfg_utils.td_modes_to_fd_modes(hlm_td, self.domain)
 
-            elif (
-                self.approximant_str == "SEOBNRv5PHM"
-                or self.approximant_str == "SEOBNRv5HM"
-            ):
+                # Step 3: Two-sided modes -> polarizations organized by m.
+                pol_m = wfg_utils.get_polarizations_from_fd_modes_m(
+                    hlm_fd, iota, parameters["phase"]
+                )
+
+            elif self.approximant_str in ["SEOBNRv5PHM", "SEOBNRv5HM"]:
                 # Step 1: generate waveform modes in L0 frame in native domain of
                 # approximant (here: TD), applying standard conditioning
                 hlm_td, iota = self.generate_TD_modes_L0_conditioned_extra_time(
@@ -1233,23 +1289,57 @@ def generate_hplus_hcross_m(
 
                 # Step 2: Transform modes to target domain.
                 hlm_fd = wfg_utils.td_modes_to_fd_modes(hlm_td, self.domain)
-            else:
-                # assert LS.SimInspiralImplementedTDApproximants(self.approximant)
+
+                # Step 3: Two-sided modes -> polarizations organized by m.
+                pol_m = wfg_utils.get_polarizations_from_fd_modes_m(
+                    hlm_fd, iota, parameters["phase"]
+                )
+
+            elif self.approximant_str in ["TEOBResumSDALI"]:
+                # Step 1: generate waveform modes in L0 frame in native domain
+                # of approximant (here: TD). Returns gwpy TimeSeries for TEOB.
+                # Then apply sigmoid start taper (TEOB fallback conditioning).
+                hlm_td, iota = self.generate_TD_modes_L0(
+                    parameters, return_lal=False
+                )
+                wfg_utils.taper_td_gwpy_modes_in_place(hlm_td, iota)
+
+                # Step 2: Transform modes to target domain.
+                # FFT without time shift, then apply phase-dependent time
+                # shift at spin_conversion_phase (peak at t≈0 for production).
+                hlm_fd_raw, resized_td_modes = (
+                    wfg_utils.td_modes_to_fd_modes_gwpy(hlm_td, self.domain)
+                )
+                hlm_fd, self._deferred_timeshift_data = (
+                    wfg_utils.apply_time_shift_to_fd_modes(
+                        hlm_fd_raw, resized_td_modes, iota,
+                        parameters["phase"], self.domain,
+                    )
+                )
+
+                # Step 3: One-sided modes -> polarizations organized by m.
+                # TEOB modes from gwpy FFT are one-sided (f >= 0 only), so we
+                # use the non-precessing symmetry h_{l,-m}(t) = (-1)^l h*_{lm}(t)
+                # to directly compute h+ and hx from the one-sided modes.
+                pol_m = wfg_utils.get_polarizations_from_onesided_fd_modes_m(
+                    hlm_fd, iota, parameters["phase"]
+                )
+
+            elif generator.domain == "freq":
                 # Step 1: generate waveform modes in L0 frame in native domain of
-                # approximant (here: TD)
-                hlm_td, iota = self.generate_TD_modes_L0(parameters)
+                # approximant (here: FD)
+                hlm_fd, iota = self.generate_FD_modes_L0(parameters)
 
                 # Step 2: Transform modes to target domain.
-                # This requires tapering of TD modes, and FFT to transform to FD.
-                wfg_utils.taper_td_modes_in_place(hlm_td)
-                hlm_fd = wfg_utils.td_modes_to_fd_modes(hlm_td, self.domain)
+                # Not required here, as approximant domain and target domain are both FD.
 
-            # Step 3: Separate negative and positive frequency parts of the modes,
-            # and add contributions according to their transformation behavior under
-            # phase shifts.
-            pol_m = wfg_utils.get_polarizations_from_fd_modes_m(
-                hlm_fd, iota, parameters["phase"]
-            )
+                # Step 3: Two-sided modes -> polarizations organized by m.
+                pol_m = wfg_utils.get_polarizations_from_fd_modes_m(
+                    hlm_fd, iota, parameters["phase"]
+                )
+
+            else:
+                raise NotImplementedError()
 
         else:
             raise NotImplementedError(
@@ -1261,7 +1351,7 @@ def generate_hplus_hcross_m(
         else:
             return pol_m
 
-    def generate_FD_modes_LO(self, parameters):  # Pending to adapt
+    def generate_FD_modes_L0(self, parameters):  # Pending to adapt
         """
         Generate FD modes in the L0 frame.
 
@@ -1329,7 +1419,7 @@ def generate_FD_modes_LO(self, parameters):  # Pending to adapt
                 f"to test that this works as intended! Ideally, add some unit tests."
             )
 
-    def generate_TD_modes_L0(self, parameters):
+    def generate_TD_modes_L0(self, parameters, return_lal=True):
         """
         Generate TD modes in the L0 frame.
 
@@ -1338,39 +1428,67 @@ def generate_TD_modes_L0(self, parameters):
         parameters: dict
             Dictionary of parameters for the waveform.
             For details see see self.generate_hplus_hcross.
+        return_lal: bool
+            If True (default), return modes as LAL COMPLEX16TimeSeries.
+            If False, return modes as gwpy TimeSeries.
 
         Returns
         -------
         hlm_td: dict
-            Dictionary with (l,m) as keys and the corresponding TD modes in lal format as
+            Dictionary with (l,m) as keys and the corresponding TD modes as
             values.
         iota: float
         """
         # TD approximants that are implemented in L0 frame. Currently tested for:
-        #   52: SEOBNRv4PHM
+        #       SEOBNRV5EHM
+        #       TEOBResumSDALI
 
         parameters_gwsignal = self._convert_parameters(
             {**parameters, "f_ref": self.f_ref}
         )
 
         generator = new_interface_get_waveform_generator(self.approximant_str)
-        hlm_td = gws_wfm.GenerateTDModes(parameters_gwsignal, generator)
-        hlms_lal = {}
+        hlm_td_gwpy = gws_wfm.GenerateTDModes(parameters_gwsignal, generator)
+        iota = parameters_gwsignal["inclination"].value
+
+        # TEOB GenerateTDModes returns dimensionless strain modes;
+        # rescale by nu * M_total / distance * (G/c^2) to get physical units.
+        if "TEOB" in self.approximant_str:
+            m1 = parameters_gwsignal["mass1"]
+            m2 = parameters_gwsignal["mass2"]
+            nu = m1 * m2 / (m1 + m2) ** 2
+            distance_rescaling = (
+                (
+                    nu
+                    * (m1 + m2)
+                    / parameters_gwsignal["distance"]
+                    * ac.G
+                    / ac.c ** 2
+                )
+                .to(u.dimensionless_unscaled)
+                .value
+            )
 
-        for key, value in hlm_td.items():
+        hlm_td = {}
+        for key, value in hlm_td_gwpy.items():
             if type(key) != str:
-                hlm_lal = lal.CreateCOMPLEX16TimeSeries(
-                    "hplus",
-                    value.epoch.value,
-                    0,
-                    value.dt.value,
-                    lal.DimensionlessUnit,
-                    len(value),
-                )
-                hlm_lal.data.data = value.value
-                hlms_lal[key] = hlm_lal
+                if "TEOB" in self.approximant_str:
+                    value = value * distance_rescaling
+                if return_lal:
+                    hlm_lal = lal.CreateCOMPLEX16TimeSeries(
+                        "hplus",
+                        value.epoch.value,
+                        0,
+                        value.dt.value,
+                        lal.DimensionlessUnit,
+                        len(value),
+                    )
+                    hlm_lal.data.data = value.value
+                    hlm_td[key] = hlm_lal
+                else:
+                    hlm_td[key] = value
 
-        return hlms_lal, parameters_gwsignal["inclination"].value
+        return hlm_td, iota
 
     def generate_TD_modes_L0_conditioned_extra_time(self, parameters):
         """
@@ -1551,15 +1669,56 @@ def generate_waveforms_parallel(
     return polarizations
 
 
-def sum_contributions_m(x_m, phase_shift=0.0):
+def sum_contributions_m(x_m, phase_shift=0.0, deferred_timeshift_data=None):
     """
     Sum the contributions over m-components, optionally introducing a phase shift.
+
+    When deferred_timeshift_data is provided and phase_shift != 0, a time shift
+    correction is applied to account for the phase-dependent epoch. The modes in
+    x_m carry a time shift computed at a reference phase; this correction adjusts
+    it to the target phase (phase_ref + phase_shift).
+
+    Parameters
+    ----------
+    x_m: dict
+        Dictionary with integer m keys and dict values containing arrays.
+    phase_shift: float
+        Phase shift to apply via exp(-1j * m * phase_shift).
+    deferred_timeshift_data: dict or None
+        If provided, must contain: resized_td_modes, iota, phase_ref, dt_ref,
+        delta_t, delta_f, frequency_array. Stored on the waveform generator
+        instance as _deferred_timeshift_data after calling
+        generate_hplus_hcross_m for TEOB approximants.
+
+    Returns
+    -------
+    result: dict
+        Dictionary with the same keys as x_m's inner dicts (e.g. "h_plus",
+        "h_cross"), containing the summed arrays.
     """
     keys = next(iter(x_m.values())).keys()
     result = {key: 0.0 for key in keys}
     for key in keys:
         for m, x in x_m.items():
             result[key] += x[key] * np.exp(-1j * m * phase_shift)
+
+    if deferred_timeshift_data is not None and phase_shift != 0.0:
+        target_phase = deferred_timeshift_data["phase_ref"] + phase_shift
+        epoch_target = wfg_utils.compute_epoch_from_resized_td_modes(
+            deferred_timeshift_data["resized_td_modes"],
+            deferred_timeshift_data["iota"],
+            target_phase,
+            deferred_timeshift_data["delta_t"],
+        )
+        dt_target = 1.0 / deferred_timeshift_data["delta_f"] + epoch_target
+        dt_correction = dt_target - deferred_timeshift_data["dt_ref"]
+        correction = np.exp(
+            -1j * 2 * np.pi * dt_correction
+            * deferred_timeshift_data["frequency_array"]
+        )
+        for key in result:
+            result[key] *= correction
+
     return result
 
 
diff --git a/dingo/gw/waveform_generator/wfg_utils.py b/dingo/gw/waveform_generator/wfg_utils.py
index 5a40ae381..a57a03556 100644
--- a/dingo/gw/waveform_generator/wfg_utils.py
+++ b/dingo/gw/waveform_generator/wfg_utils.py
@@ -1,6 +1,8 @@
 import numpy as np
 import lal
 import lalsimulation as LS
+from lalsimulation.gwsignal.core import conditioning_subroutines as cond
+from gwpy.timeseries import TimeSeries
 
 
 def linked_list_modes_to_dict_modes(hlm_ll):
@@ -44,12 +46,12 @@ def get_tapering_window_for_complex_time_series(h, tapering_flag: int = 1):
     )
     h_tapered.data.data = h.data.data.copy().real
     LS.SimInspiralREAL8WaveTaper(h_tapered.data, tapering_flag)
-    eps = 1e-20 * np.max(np.abs(h.data.data))
+    # in case h.data.data is an array of 0's this makes the window return ones
+    eps = 1e-20 * np.max(np.abs(h.data.data)) if np.max(np.abs(h.data.data)) > 0 else 1e-20
     window = (np.abs(h_tapered.data.data) + eps) / (np.abs(h.data.data.real) + eps)
     # FIXME: using eps for numerical stability is not really robust here
     return window
 
-
 def taper_td_modes_in_place(hlm_td, tapering_flag: int = 1):
     """
     Taper the time domain modes in place.
@@ -69,6 +71,90 @@ def taper_td_modes_in_place(hlm_td, tapering_flag: int = 1):
         window = get_tapering_window_for_complex_time_series(h, tapering_flag)
         h.data.data *= window
 
+def get_tapering_window_for_real_array(signal, taper_kind="start"):
+    """Recover the multiplicative gwsignal taper window for a real signal.
+
+    Mirrors `get_tapering_window_for_complex_time_series` (LAL path) for the
+    gwsignal/TEOB path: runs `cond.taper_gwpy_timeseries` on a copy and
+    derives the window as |tapered| / |original| with an eps regularizer.
+    Returning the window (instead of the tapered signal) lets the caller
+    apply it independently to the real and imaginary parts of complex modes.
+
+    Parameters
+    ----------
+    signal: np.ndarray
+        Real-valued 1D array to derive the window from.
+    taper_kind: str
+        Forwarded to `cond.taper_gwpy_timeseries`. Default "start".
+
+    Returns
+    -------
+    window: np.ndarray
+        Real-valued window of the same length as `signal`. All ones if the
+        signal is too short or all-zero (matching gwsignal's no-op behavior).
+    """
+    ts = TimeSeries(signal.copy(), dt=1.0)
+    ts_tapered = cond.taper_gwpy_timeseries(ts, taper_kind)
+    abs_signal = np.abs(signal)
+    max_abs = np.max(abs_signal) if abs_signal.size else 0.0
+    eps = 1e-20 * max_abs if max_abs > 0 else 1e-20
+    return (np.abs(ts_tapered.value) + eps) / (abs_signal + eps)
+
+
+def taper_td_gwpy_modes_in_place(hlm_td, iota, phase=0.0):
+    """
+    Apply sigmoid start taper to gwpy TimeSeries modes in place.
+
+    This replicates the fallback conditioning used by gwsignal for
+    approximants with f_ref_spin=False (e.g. TEOBResumSDALI).
+    See generate_conditioned_td_waveform_from_td_fallback in
+    lalsimulation/gwsignal/core/waveform_conditioning.py:35-64.
+
+    The gwsignal reference applies `taper_gwpy_timeseries` to the combined
+    polarizations hp and hc, which can yield different windows because
+    find_peaks finds different peaks. We therefore derive `W_hp` and `W_hc`
+    separately via `get_tapering_window_for_real_array` and apply `W_hp` to
+    Re(h_lm) and `W_hc` to Im(h_lm) for each mode. For non-precessing
+    systems the modes are purely real in the co-rotating frame, so Re(h_lm)
+    and Im(h_lm) map cleanly to the two polarizations, giving agreement with
+    the reference.
+
+    Parameters
+    ----------
+    hlm_td: dict
+        Dictionary with (l,m) keys and gwpy TimeSeries objects.
+    iota: float
+        Inclination angle in radians.
+    phase: float
+        Reference phase for computing the combined polarizations used to
+        determine the taper window. Defaults to 0.0 since modes should be
+        phase-independent.
+    """
+    max_len = max(len(h) for h in hlm_td.values())
+
+    # Sum modes to get combined hp and hc (unconditioned). The reference
+    # path runs the taper on these summed polarizations.
+    h_complex = np.zeros(max_len, dtype=complex)
+    for (l, m), h in hlm_td.items():
+        data = h.value.copy()
+        if len(data) < max_len:
+            data = np.pad(data, (0, max_len - len(data)))
+        ylm = lal.SpinWeightedSphericalHarmonic(iota, np.pi / 2 - phase, -2, l, m)
+        h_complex += ylm * data
+
+    W_hp = get_tapering_window_for_real_array(h_complex.real)
+    W_hc = get_tapering_window_for_real_array(-h_complex.imag)
+
+    for (l, m) in hlm_td:
+        data = hlm_td[(l, m)].value.copy()
+        if len(data) < max_len:
+            data = np.pad(data, (0, max_len - len(data)))
+        data = W_hp * data.real + 1j * W_hc * data.imag
+        hlm_td[(l, m)] = TimeSeries(
+            data,
+            t0=hlm_td[(l, m)].t0,
+            dt=hlm_td[(l, m)].dt,
+        )
 
 def td_modes_to_fd_modes(hlm_td, domain):
     """
@@ -137,6 +223,163 @@ def td_modes_to_fd_modes(hlm_td, domain):
 
     return hlm_fd
 
+def td_modes_to_fd_modes_gwpy(hlm_td, domain):
+    """
+    Transform dict of gwpy TD modes to dict of one-sided FD modes via FFT.
+    The td modes are expected to be conditioned (tapered/filtered).
+
+    Replicates the FFT pipeline in gwsignal's
+    generate_conditioned_fd_waveform_from_td (waveform_conditioning.py:510-522):
+      resize → fft → normalize by 1/(2*df)
+
+    No time shift is applied to the FD modes. The caller is responsible for
+    computing and applying the time shift (which depends on the target phase
+    via np.argmax of the reconstructed hp). This enables the "deferred time
+    shift" approach where modes are FFT'd once (phase-independently) and the
+    phase-dependent time shift is applied at resummation time.
+
+    Also returns the resized TD mode data (numpy arrays) needed for the
+    deferred time shift computation.
+
+    Since the modes are complex-valued and gwpy's fft uses rfft (real input
+    only), each mode is split into real and imaginary parts which are FFT'd
+    separately and recombined: H_lm(f) = FFT(Re) + i*FFT(Im).
+
+    Parameters
+    ----------
+    hlm_td: dict
+        Dictionary with (l,m) keys and gwpy TimeSeries objects for the
+        corresponding conditioned modes.
+    domain: dingo.gw.domains.UniformFrequencyDomain
+        Target domain after FFT.
+
+    Returns
+    -------
+    hlm_fd: dict
+        Dictionary with (l,m) keys and numpy arrays with the corresponding
+        one-sided FD modes on [0, df, ..., f_max]. No time shift applied.
+    resized_td_modes: dict
+        Dictionary with (l,m) keys and numpy complex arrays of the resized
+        TD mode data, needed for deferred time shift computation.
+    """
+    hlm_fd = {}
+    resized_td_modes = {}
+
+    delta_f = domain.delta_f
+    delta_t = 0.5 / domain.f_max
+    f_nyquist = domain.f_max
+    chirplen = int(2 * f_nyquist / delta_f)
+    frequency_array = domain()  # [0, df, ..., f_max]
+
+    for (l, m), h in hlm_td.items():
+        # Resize to chirplen (waveform_conditioning.py:510)
+        start_id = len(h) - chirplen
+        h_resized = cond.resize_gwpy_timeseries(h, start_id, chirplen)
+        resized_td_modes[(l, m)] = np.asarray(h_resized).copy()
+
+        # Split complex mode into real and imaginary TimeSeries for gwpy fft
+        h_re = TimeSeries(np.asarray(h_resized).real, dt=delta_t)
+        h_im = TimeSeries(np.asarray(h_resized).imag, dt=delta_t)
+
+        # FFT + normalize (waveform_conditioning.py:513-521)
+        hf_re = h_re.fft()
+        hf_re = hf_re / (2 * hf_re.df)
+
+        hf_im = h_im.fft()
+        hf_im = hf_im / (2 * hf_im.df)
+
+        hlm_fd[(l, m)] = (hf_re.value[:len(frequency_array)]
+                          + 1j * hf_im.value[:len(frequency_array)])
+
+    return hlm_fd, resized_td_modes
+
+
+def compute_epoch_from_resized_td_modes(resized_td_modes, iota, phase, delta_t):
+    """Compute the epoch (time of peak) from resized TD modes at a given phase.
+
+    Reconstructs the real h+ polarization from the resized TD mode data using
+    spin-weighted spherical harmonics at the specified phase, then finds the
+    peak via np.argmax. This matches the epoch convention used by gwsignal's
+    resize_gwpy_timeseries (which calls np.argmax on the real hp).
+
+    Parameters
+    ----------
+    resized_td_modes: dict
+        Dictionary with (l,m) keys and numpy complex arrays of resized TD modes.
+    iota: float
+        Inclination angle in radians.
+    phase: float
+        Reference phase for spherical harmonic evaluation.
+    delta_t: float
+        Time step of the TD data.
+
+    Returns
+    -------
+    epoch: float
+        The epoch value (negative of peak index times delta_t).
+    """
+    chirplen = len(next(iter(resized_td_modes.values())))
+    hp = np.zeros(chirplen)
+    for (l, m), h_data in resized_td_modes.items():
+        ylm = lal.SpinWeightedSphericalHarmonic(iota, np.pi / 2 - phase, -2, l, m)
+        hp += (ylm * h_data).real
+    return -np.argmax(hp) * delta_t
+
+
+def apply_time_shift_to_fd_modes(hlm_fd_raw, resized_td_modes, iota, phase, domain):
+    """Apply a phase-dependent time shift to raw (unshifted) FD modes.
+
+    Computes the epoch by reconstructing hp from resized TD modes at the given
+    phase (matching gwsignal's resize_gwpy_timeseries convention), then applies
+    the time shift exp(-i 2pi dt f) to each FD mode.
+
+    Also returns a deferred_timeshift_data dict that can be passed to
+    sum_contributions_m to correct the time shift when a phase_shift is applied.
+
+    Parameters
+    ----------
+    hlm_fd_raw : dict
+        Dictionary with (l,m) keys and numpy arrays of unshifted FD modes.
+    resized_td_modes : dict
+        Dictionary with (l,m) keys and numpy complex arrays of resized TD modes.
+    iota : float
+        Inclination angle in radians.
+    phase : float
+        Reference phase for epoch computation.
+    domain : dingo.gw.domains.UniformFrequencyDomain
+        Frequency domain.
+
+    Returns
+    -------
+    hlm_fd : dict
+        Dictionary with (l,m) keys and time-shifted FD mode arrays.
+    deferred_timeshift_data : dict
+        Data needed by sum_contributions_m for deferred time shift correction.
+    """
+    delta_t = 0.5 / domain.f_max
+    delta_f = domain.delta_f
+    frequency_array = domain()
+
+    epoch = compute_epoch_from_resized_td_modes(
+        resized_td_modes, iota, phase, delta_t
+    )
+    dt = 1.0 / delta_f + epoch
+    time_shift = np.exp(-1j * 2 * np.pi * dt * frequency_array)
+
+    hlm_fd = {lm: hf * time_shift for lm, hf in hlm_fd_raw.items()}
+
+    deferred_timeshift_data = {
+        "resized_td_modes": resized_td_modes,
+        "iota": iota,
+        "phase_ref": phase,
+        "dt_ref": dt,
+        "delta_t": delta_t,
+        "delta_f": delta_f,
+        "frequency_array": frequency_array,
+    }
+
+    return hlm_fd, deferred_timeshift_data
+
 
 def get_polarizations_from_fd_modes_m(hlm_fd, iota, phase):
     pol_m = {}
@@ -174,6 +417,58 @@ def get_polarizations_from_fd_modes_m(hlm_fd, iota, phase):
     return pol_m
 
 
+def get_polarizations_from_onesided_fd_modes_m(hlm_fd, iota, phase):
+    """Combine one-sided (f >= 0) FD modes into polarizations organized by m.
+
+    For one-sided FD modes from non-precessing waveform models, the polarizations
+    are recovered using the symmetry h_{l,-m}(t) = (-1)^l h*_{lm}(t), which
+    relates the negative-frequency content to the positive-frequency modes:
+
+        h+(f) = 0.5 sum_{l,m} [Y_{lm} + (-1)^l conj(Y_{l,-m})] H_{lm}(f)
+        hx(f) = 0.5i sum_{l,m} [Y_{lm} - (-1)^l conj(Y_{l,-m})] H_{lm}(f)
+
+    where Y_{lm} = _{-2}Y_{lm}(iota, pi/2 - phase).
+
+    Each pol_m[m] transforms as exp(-1j * m * phase) under phase shifts.
+
+    Parameters
+    ----------
+    hlm_fd : dict
+        Dictionary with (l, m) keys and one-sided FD mode arrays (numpy arrays
+        or gwpy FrequencySeries) as values, defined for f >= 0 only.
+    iota : float
+        Inclination angle.
+    phase : float
+        Reference phase.
+
+    Returns
+    -------
+    pol_m : dict
+        Dictionary with integer m as keys and
+        {"h_plus": array, "h_cross": array} as values.
+    """
+    pol_m = {}
+    polarizations = ["h_plus", "h_cross"]
+
+    for (l, m), h in hlm_fd.items():
+        if m not in pol_m:
+            pol_m[m] = {k: 0.0 for k in polarizations}
+
+        # Extract numpy array from gwpy FrequencySeries if needed
+        h_data = h.value if hasattr(h, "value") else h
+
+        ylm = lal.SpinWeightedSphericalHarmonic(iota, np.pi / 2 - phase, -2, l, m)
+        ylm_neg_conj = np.conj(
+            lal.SpinWeightedSphericalHarmonic(iota, np.pi / 2 - phase, -2, l, -m)
+        )
+
+        sign = (-1) ** l
+        pol_m[m]["h_plus"] += 0.5 * (ylm + (sign * ylm_neg_conj)) * h_data
+        pol_m[m]["h_cross"] += 0.5j * (ylm - (sign * ylm_neg_conj)) * h_data
+
+    return pol_m
+
+
 def get_starting_frequency_for_SEOBRNRv5_conditioning(parameters):
     """
     Compute starting frequency needed for having 3 extra cycles for tapering the TD modes.
@@ -299,3 +594,5 @@ def taper_td_modes_for_SEOBRNRv5_extra_time(
 
     # return timeseries
     return h_return
+
+
diff --git a/pyproject.toml b/pyproject.toml
index fa510cd77..d774c8415 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -107,6 +107,10 @@ typing = [
     "scipy-stubs>=1.15.3.0",
 ]
 
+teobresums = [
+    "teobresums",
+]
+
 [project.scripts]
 dingo_append_training_stage = "dingo.gw.training:append_stage"
 dingo_build_svd = "dingo.gw.dataset.utils:build_svd_cli"
diff --git a/tests/gw/waveform_generator/test_wfg_m.py b/tests/gw/waveform_generator/test_wfg_m.py
index cc8769939..ce3245430 100644
--- a/tests/gw/waveform_generator/test_wfg_m.py
+++ b/tests/gw/waveform_generator/test_wfg_m.py
@@ -37,14 +37,15 @@ def uniform_fd_domain():
     return domain
 
 
-@pytest.fixture(params=["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM"])
+@pytest.fixture(params=["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM", "TEOBResumSDALI"])
 def approximant(request):
     return request.param
 
 
 @pytest.fixture
 def intrinsic_prior(approximant):
-    if "PHM" in approximant:
+    if approximant in ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM"]:
+        # quasi-circular precessing-spin 
         intrinsic_dict = {
             "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)",
             "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)",
@@ -61,8 +62,21 @@ def intrinsic_prior(approximant):
             "phi_jl": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")',
             "geocent_time": 0.0,
         }
-    else:
-        # Aligned spins
+    elif approximant in ["SEOBNRv5HM"]:
+        # quasi-circular aligned-spin
+        intrinsic_dict = {
+            "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)",
+            "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)",
+            "mass_ratio": "bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0)",
+            "chirp_mass": "bilby.gw.prior.UniformInComponentsChirpMass(minimum=25.0, maximum=100.0)",
+            "luminosity_distance": 1000.0,
+            "theta_jn": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)",
+            "phase": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")',
+            "chi_1": 'bilby.gw.prior.AlignedSpin(name="chi_1", a_prior=Uniform(minimum=0, maximum=0.99))',
+            "chi_2": 'bilby.gw.prior.AlignedSpin(name="chi_2", a_prior=Uniform(minimum=0, maximum=0.99))',
+            "geocent_time": 0.0,
+        }
+    elif approximant in ["SEOBNRv5EHM", "TEOBResumSDALI"]:
         intrinsic_dict = {
             "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)",
             "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)",
@@ -73,15 +87,20 @@ def intrinsic_prior(approximant):
             "phase": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")',
             "chi_1": 'bilby.gw.prior.AlignedSpin(name="chi_1", a_prior=Uniform(minimum=0, maximum=0.99))',
             "chi_2": 'bilby.gw.prior.AlignedSpin(name="chi_2", a_prior=Uniform(minimum=0, maximum=0.99))',
+            "eccentricity": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.3)",
+            "relativistic_anomaly": "bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi)",
             "geocent_time": 0.0,
         }
+    else:
+        raise ValueError(f"Unimplemented approximant {approximant}")
+
     prior = build_prior_with_defaults(intrinsic_dict)
     return prior
 
 
 @pytest.fixture
 def wfg(uniform_fd_domain, approximant):
-    if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM"]:
+    if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM", "TEOBResumSDALI"]:
         wfg_class = NewInterfaceWaveformGenerator
     else:
         wfg_class = WaveformGenerator
@@ -115,7 +134,7 @@ def tolerances(approximant):
         # get should not have a big effect in practice.
         return 2e-2, 1e-5
 
-    elif approximant == "SEOBNRv4PHM":
+    elif approximant in ["SEOBNRv4PHM", "SEOBNRv5EHM", "TEOBResumSDALI"]:
         # The mismatches are typically be of order 1e-5. This is exclusively due to
         # different tapering. The reference polarizations are tapered and FFTed on the
         # level of polarizations, while for generate_hplus_hcross_m, the tapering and FFT
@@ -132,13 +151,21 @@ def tolerances(approximant):
         return 1e-5, 1e-5
 
 
-# Uncomment to test only one approximant.
+approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM"]
+
 try:
-    import pyseobnr
+    import pyseobnr  # noqa: F401
 
-    approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM"]
+    approximant_list += ["SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"]
 except ImportError:
-    approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM"]
+    pass
+
+try:
+    import EOBRun_module  # noqa: F401
+
+    approximant_list += ["TEOBResumSDALI"]
+except ImportError:
+    pass
 
 
 @pytest.mark.parametrize("approximant", approximant_list)
@@ -149,7 +176,13 @@ def test_generate_hplus_hcross_m(intrinsic_prior, wfg, num_evaluations, toleranc
         phase_shift = np.random.uniform(high=2 * np.pi)
 
         pol_m = wfg.generate_hplus_hcross_m(p)
-        pol = sum_contributions_m(pol_m, phase_shift=phase_shift)
+        pol = sum_contributions_m(
+            pol_m,
+            phase_shift=phase_shift,
+            deferred_timeshift_data=getattr(
+                wfg, "_deferred_timeshift_data", None
+            ),
+        )
         pol_ref = wfg.generate_hplus_hcross({**p, "phase": p["phase"] + phase_shift})
 
         mismatches.append(
diff --git a/tests/gw/waveform_generator/test_wfg_mfd.py b/tests/gw/waveform_generator/test_wfg_mfd.py
index 297df73bd..badd9ac59 100644
--- a/tests/gw/waveform_generator/test_wfg_mfd.py
+++ b/tests/gw/waveform_generator/test_wfg_mfd.py
@@ -24,36 +24,70 @@ def mfd():
     return domain
 
 
-@pytest.fixture(params=["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM"])
+@pytest.fixture(params=["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM", "TEOBResumSDALI"])
 def approximant(request):
     return request.param
 
 
 @pytest.fixture
 def intrinsic_prior(approximant):
-    intrinsic_dict = {
-        "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)",
-        "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)",
-        "mass_ratio": "bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0)",
-        "chirp_mass": "bilby.gw.prior.UniformInComponentsChirpMass(minimum=25.0, maximum=100.0)",
-        "luminosity_distance": 1000.0,
-        "theta_jn": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)",
-        "phase": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")',
-        "a_1": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.99)",
-        "a_2": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.99)",
-        "tilt_1": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)",
-        "tilt_2": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)",
-        "phi_12": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")',
-        "phi_jl": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")',
-        "geocent_time": "bilby.core.prior.Uniform(minimum=-0.1, maximum=0.1)",
-    }
+    if approximant in ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM"]:
+        # quasi-circular precessing-spin 
+        intrinsic_dict = {
+            "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)",
+            "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)",
+            "mass_ratio": "bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0)",
+            "chirp_mass": "bilby.gw.prior.UniformInComponentsChirpMass(minimum=25.0, maximum=100.0)",
+            "luminosity_distance": 1000.0,
+            "theta_jn": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)",
+            "phase": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")',
+            "a_1": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.99)",
+            "a_2": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.99)",
+            "tilt_1": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)",
+            "tilt_2": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)",
+            "phi_12": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")',
+            "phi_jl": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")',
+            "geocent_time": 0.0,
+        }
+    elif approximant in ["SEOBNRv5HM"]:
+        # quasi-circular aligned-spin
+        intrinsic_dict = {
+            "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)",
+            "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)",
+            "mass_ratio": "bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0)",
+            "chirp_mass": "bilby.gw.prior.UniformInComponentsChirpMass(minimum=25.0, maximum=100.0)",
+            "luminosity_distance": 1000.0,
+            "theta_jn": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)",
+            "phase": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")',
+            "chi_1": 'bilby.gw.prior.AlignedSpin(name="chi_1", a_prior=Uniform(minimum=0, maximum=0.99))',
+            "chi_2": 'bilby.gw.prior.AlignedSpin(name="chi_2", a_prior=Uniform(minimum=0, maximum=0.99))',
+            "geocent_time": 0.0,
+        }
+    elif approximant in ["SEOBNRv5EHM", "TEOBResumSDALI"]:
+        intrinsic_dict = {
+            "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)",
+            "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)",
+            "mass_ratio": "bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0)",
+            "chirp_mass": "bilby.gw.prior.UniformInComponentsChirpMass(minimum=25.0, maximum=100.0)",
+            "luminosity_distance": 1000.0,
+            "theta_jn": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)",
+            "phase": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")',
+            "chi_1": 'bilby.gw.prior.AlignedSpin(name="chi_1", a_prior=Uniform(minimum=0, maximum=0.99))',
+            "chi_2": 'bilby.gw.prior.AlignedSpin(name="chi_2", a_prior=Uniform(minimum=0, maximum=0.99))',
+            "eccentricity": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.3)",
+            "relativistic_anomaly": "bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi)",
+            "geocent_time": 0.0,
+        }
+    else:
+        raise ValueError(f"Unimplemented approximant {approximant}")
+
     prior = build_prior_with_defaults(intrinsic_dict)
     return prior
 
 
 @pytest.fixture
 def wfg_mfd(mfd, approximant):
-    if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM"]:
+    if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM", "TEOBResumSDALI"]:
         wfg_class = NewInterfaceWaveformGenerator
     else:
         wfg_class = WaveformGenerator
@@ -68,7 +102,7 @@ def wfg_mfd(mfd, approximant):
 
 @pytest.fixture
 def wfg_ufd(mfd, approximant):
-    if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM"]:
+    if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM", "TEOBResumSDALI"]:
         wfg_class = NewInterfaceWaveformGenerator
     else:
         wfg_class = WaveformGenerator
@@ -100,9 +134,11 @@ def tolerances(approximant):
 
 # Uncomment to test only one approximant.
 try:
-    import pyseobnr
+    # import pyseobnr
+    import EOBRun_module
 
-    approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM"]
+    # approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM", "TEOBResumSDALI"]
+    approximant_list = ["TEOBResumSDALI"]
 except ImportError:
     approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM"]