Second-Order Phase Reduction

• Second-order phase reduction is a method that extends classic phase reduction by incorporating O(ε²) amplitude feedback, capturing higher harmonics and nonpairwise interactions.
• It provides detailed derivations for models like Stuart–Landau networks and delay-coupled systems, revealing triplet and mediated coupling effects absent in first-order approximations.
• This advanced approach overcomes limitations of first-order models by accurately predicting synchrony, chimera states, and bifurcation phenomena in weakly coupled oscillator systems.

Searching arXiv for recent and foundational papers on second-order phase reduction to ground the article in the current literature. Second-order phase reduction is the extension of phase-only modeling from the first-order approximation

$\dot{\varphi}=\omega+\varepsilon Q_1+\cdots$

to

$$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$

for weakly perturbed or weakly coupled limit-cycle oscillators whose long-time dynamics lies on an attracting invariant torus. In this setting, the $O(\varepsilon^2)$ term represents the feedback of $O(\varepsilon)$ amplitude deviations onto phase dynamics; for delay-coupled systems, it also incorporates first-order corrections to the delayed history. Across mean-field complex Ginzburg–Landau dynamics, Stuart–Landau networks, generic two-dimensional oscillators, delay-coupled oscillators, and numerically reconstructed phase models, second-order reduction consistently yields interaction structures and bifurcations that are absent from first-order Kuramoto-, Winfree-, or Kuramoto–Sakaguchi-type descriptions (Rosenblum et al., 2018, León et al., 2019, Bick et al., 2024).

1. First-order reduction and the need for a second order

Classical phase reduction starts from a stable limit cycle and replaces the full state dynamics by an equation for asymptotic phase alone. In the standard weak-coupling setting, one retains only the leading $O(\varepsilon)$ correction. This produces the familiar pairwise first-harmonic coupling structure. For a perturbed oscillator, the literature summarized here repeatedly uses the expansion

$$\dot\varphi=\omega+\varepsilon Q_1(\varphi,\psi)+\varepsilon^2Q_2(\varphi,\psi)+\ldots,$$

with the first-order term reducing to the usual Winfree or Kuramoto–Sakaguchi form when evaluated on the unperturbed limit cycle (Rosenblum et al., 2018, Mau et al., 2023).

The principal reason to go beyond first order is structural, not merely quantitative. In the mean-field complex Ginzburg–Landau equation (MF-CGLE), first-order phase reduction yields the disorder-free Kuramoto–Sakaguchi model

$$\dot{\theta}_j=\Omega+\epsilon\eta R\sin(\Psi-\theta_j+\alpha),$$

which preserves only pairwise sinusoidal mean-field coupling through $Z_1$. In the identical-oscillator case, this model is too degenerate: it yields only full synchrony or incoherence, whereas the original MF-CGLE already exhibits nonuniform incoherent states (NUIS), quasiperiodic partial synchrony (QPS), and routes toward pure collective chaos (León et al., 2019).

The same obstruction appears in other settings. For identical Stuart–Landau oscillator networks, the first-order phase model can be transformed, by going to a corotating frame and rescaling time, into a form independent of the coupling strength $\varepsilon$; it therefore cannot describe any coupling-strength dependence of chimera shape, even though such dependence is present in the full oscillator network (Mau et al., 2023). For delay-coupled oscillators, first order identifies delay with a phase lag, because it evaluates the delayed state on the uncoupled periodic orbit. This equivalence fails at second order, where the delayed coupling samples the first-order corrected state and history rather than the unperturbed cycle (Bick et al., 2024, Bick et al., 31 Oct 2025). In globally coupled two-group Stuart–Landau populations, first-order mean-field theory is also constrained by Watanabe–Strogatz integrability and forbids stable three-cluster $(2+1)$ states that appear in the second-order description (Baibolatov et al., 3 Jun 2026).

A common misconception is therefore that second-order phase reduction merely renormalizes coupling coefficients. The cited work shows that the more important effect is often a change in the admissible functional form of the reduced dynamics: higher harmonics, mediated effective couplings, and genuinely nonpairwise interactions emerge at $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$0 (León et al., 2019, Mau et al., 2023, Genge et al., 2020).

2. Geometric foundations: isochrons, isostables, and the invariant torus

The geometric basis of phase reduction is the asymptotic phase defined by isochrons. For the uncoupled Stuart–Landau oscillator

$$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$1

with polar form

$$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$2

the asymptotic phase is given exactly by

$$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$3

On the limit cycle $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$4, one has $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$5, and the nonisochronicity parameter $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$6 measures the tilt of isochrons away from radial lines (León et al., 2019).

A closely related exact construction appears for the Stuart–Landau network formulation with heterogeneity and arbitrary coupling topology. There the phase–isostable transformation is

$$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$7

so that the uncoupled dynamics becomes exactly linear,

$$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$8

with $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$9 the Floquet exponent of the radial mode. In these coordinates, $O(\varepsilon^2)$0 corresponds to the limit cycle, $O(\varepsilon^2)$1 is asymptotic phase, and $O(\varepsilon^2)$2 is an isostable amplitude coordinate (Mau et al., 2023).

For generic smooth two-dimensional oscillators, an analogous phase–amplitude normal form exists locally near the periodic orbit: $O(\varepsilon^2)$3 This formulation makes precise the usual statement that amplitudes are enslaved: on the attracting invariant torus, one writes

$O(\varepsilon^2)$4

or, in perturbative form,

$O(\varepsilon^2)$5

Second-order phase reduction is the first order at which the dependence of the phase equation on these slaved amplitude coordinates affects the reduced phase dynamics (Mau et al., 2023).

When the uncoupled limit cycle is not circular but has phase-dependent amplitude, the geometry changes qualitatively. In the class of oscillators with limit cycle

$O(\varepsilon^2)$6

the transformed system generally loses $O(\varepsilon^2)$7-equivariance for $O(\varepsilon^2)$8. In that case, the reduced phase equations are no longer functions only of phase differences; absolute phases enter already at first order, and the first nontrivial torus correction must be obtained from a transport-type PDE rather than from a simple algebraic slaving relation (Bick et al., 2023).

3. Derivation strategies for second-order phase equations

A recurring analytic pattern is to transform the full dynamics into phase–amplitude or phase–isostable coordinates, expand the invariant torus or slaved amplitudes in powers of the coupling parameter, solve the resulting linear invariance equation for the first amplitude correction, and substitute that correction back into the phase equation. In the MF-CGLE, this is done directly in exact isochron coordinates, leading to the expansion

$O(\varepsilon^2)$9

with $O(\varepsilon)$0 and

$O(\varepsilon)$1

because $O(\varepsilon)$2 for the radial dynamics $O(\varepsilon)$3. The second-order phase equation then follows from the term $O(\varepsilon)$4 in the phase expansion (León et al., 2019).

For arbitrary Stuart–Landau networks with coupling matrix $O(\varepsilon)$5, the reduced phase model is sought in the form

$O(\varepsilon)$6

The first-order amplitude correction $O(\varepsilon)$7 is determined by the linear PDE

$O(\varepsilon)$8

which is solved Fourier mode by Fourier mode using an explicit resolvent operator $O(\varepsilon)$9. The second-order phase term is then

$O(\varepsilon)$0

and, after simplification, decomposes into pairwise and triplet contributions whose coefficients are determined explicitly by the original network architecture (Mau et al., 2023).

For generic coupled two-dimensional oscillators, the same logic yields an abstract second-order formula. Writing

$O(\varepsilon)$1

and assuming $O(\varepsilon)$2 on the invariant torus, one obtains

$O(\varepsilon)$3

Here $O(\varepsilon)$4 is the Green operator solving the linear advection–decay equation for the slaved amplitude. Because $O(\varepsilon)$5, the Fourier-space denominators never vanish, so this construction has no small divisors in the two-dimensional setting (Mau et al., 2023).

Delay-coupled systems require a different, though closely related, mechanism. The delay differential equation is rewritten as an ODE for the current state coupled to a transport equation for the history field $O(\varepsilon)$6, with $O(\varepsilon)$7. One then expands the torus embedding $O(\varepsilon)$8, the history embedding $O(\varepsilon)$9, and the reduced phase flow $$\dot\varphi=\omega+\varepsilon Q_1(\varphi,\psi)+\varepsilon^2Q_2(\varphi,\psi)+\ldots,$$0 order by order. The decisive second-order fact is that the homological equation for the delayed history contains $$\dot\varphi=\omega+\varepsilon Q_1(\varphi,\psi)+\varepsilon^2Q_2(\varphi,\psi)+\ldots,$$1, so second-order phase dynamics depends on the first-order corrected history rather than only on the unperturbed orbit (Bick et al., 2024, Bick et al., 31 Oct 2025).

A numerical alternative is available when explicit phase coordinates are unavailable. The reconstruction method of “Numerical phase reduction beyond the first order approximation” computes the isochron phase of a perturbed state by relaxing it under the autonomous dynamics onto the original limit cycle, differentiates the reconstructed phase time series, fits the effective coupling function $$\dot\varphi=\omega+\varepsilon Q_1(\varphi,\psi)+\varepsilon^2Q_2(\varphi,\psi)+\ldots,$$2 on the torus, and then extracts $$\dot\varphi=\omega+\varepsilon Q_1(\varphi,\psi)+\varepsilon^2Q_2(\varphi,\psi)+\ldots,$$3 and $$\dot\varphi=\omega+\varepsilon Q_1(\varphi,\psi)+\varepsilon^2Q_2(\varphi,\psi)+\ldots,$$4 by polynomial fitting in $$\dot\varphi=\omega+\varepsilon Q_1(\varphi,\psi)+\varepsilon^2Q_2(\varphi,\psi)+\ldots,$$5. This is a numerical effective reduction rather than a closed-form asymptotic derivation (Rosenblum et al., 2018).

4. Interaction structure generated at second order

The main structural result of second-order reduction is that phase dynamics is no longer generically pairwise. In the MF-CGLE, the second-order reduced phase equation is

$$\dot\varphi=\omega+\varepsilon Q_1(\varphi,\psi)+\varepsilon^2Q_2(\varphi,\psi)+\ldots,$$6

with

$$\dot\varphi=\omega+\varepsilon Q_1(\varphi,\psi)+\varepsilon^2Q_2(\varphi,\psi)+\ldots,$$7

The term $$\dot\varphi=\omega+\varepsilon Q_1(\varphi,\psi)+\varepsilon^2Q_2(\varphi,\psi)+\ldots,$$8 only renormalizes first-harmonic pairwise coupling, but the terms

$$\dot\varphi=\omega+\varepsilon Q_1(\varphi,\psi)+\varepsilon^2Q_2(\varphi,\psi)+\ldots,$$9

are genuinely nonpairwise. They can be rewritten as

$$\dot{\theta}_j=\Omega+\epsilon\eta R\sin(\Psi-\theta_j+\alpha),$$0

which makes their three-body character explicit (León et al., 2019).

In arbitrary Stuart–Landau networks, the second-order term $$\dot{\theta}_j=\Omega+\epsilon\eta R\sin(\Psi-\theta_j+\alpha),$$1 comprises two effectively pairwise and two genuine triplet terms. The coefficients of the direct second-order corrections are proportional to

$$\dot{\theta}_j=\Omega+\epsilon\eta R\sin(\Psi-\theta_j+\alpha),$$2

whereas the mediated terms are proportional to

$$\dot{\theta}_j=\Omega+\epsilon\eta R\sin(\Psi-\theta_j+\alpha),$$3

The first class requires oscillator $$\dot{\theta}_j=\Omega+\epsilon\eta R\sin(\Psi-\theta_j+\alpha),$$4 to be structurally connected to both $$\dot{\theta}_j=\Omega+\epsilon\eta R\sin(\Psi-\theta_j+\alpha),$$5 and $$\dot{\theta}_j=\Omega+\epsilon\eta R\sin(\Psi-\theta_j+\alpha),$$6; the second produces mediated or “virtual pairwise” interactions through node $$\dot{\theta}_j=\Omega+\epsilon\eta R\sin(\Psi-\theta_j+\alpha),$$7, including effective coupling between $$\dot{\theta}_j=\Omega+\epsilon\eta R\sin(\Psi-\theta_j+\alpha),$$8 and $$\dot{\theta}_j=\Omega+\epsilon\eta R\sin(\Psi-\theta_j+\alpha),$$9 even if $Z_1$0. The original pairwise graph therefore lifts to a hypernetwork in the reduced phase dynamics (Mau et al., 2023).

The hypergraph viewpoint is developed further for oscillators with phase-dependent amplitude. In the graph-coupled case, second-order interactions are encoded by two directed 3-uniform hypergraphs with tensors

$Z_1$1

Some terms remain pairwise corrections along original edges, some create virtual pairwise edges weighted by two-step paths, and some are genuine triplet interactions. This gives a derived, rather than ad hoc, hypergraph phase model (Bick et al., 2023).

The appearance of indirect and triplet terms is not restricted to mean-field settings. For a chain of three Stuart–Landau oscillators, the explicit second-order phase equations contain second harmonics such as $Z_1$2, mediated effective couplings such as $Z_1$3 even though oscillators $Z_1$4 and $Z_1$5 are not directly linked, and genuine triplet terms such as $Z_1$6. Consequently, phase connectivity and structural connectivity no longer coincide once higher-order corrections are included (Genge et al., 2020).

Delay-coupled systems exhibit a different but equally characteristic second-order structure. For two identical delay-coupled Stuart–Landau oscillators, the second-order phase-difference equation contains both explicit $Z_1$7-dependent coefficients and a new second harmonic $Z_1$8. This is the first order at which one can rigorously say that delay is not equivalent to phase lag: dependence on $Z_1$9 appears not only through $\varepsilon$0, but through coefficients such as $\varepsilon$1 and $\varepsilon$2 (Bick et al., 2024, Bick et al., 31 Oct 2025).

The MF-CGLE results also suggest a broader pattern: truncation at order $\varepsilon$3 yields phase interactions involving up to $\varepsilon$4-body terms and order parameters $\varepsilon$5 only up to $\varepsilon$6. Second order yields three-body terms, third order yields four-body terms, and the many-body hierarchy continues accordingly (León et al., 2019).

5. Dynamical consequences and representative applications

The dynamical payoff of second-order reduction is clearest in problems where the first-order model is degenerate. In the MF-CGLE, second order distinguishes full synchrony, the uniform incoherent state (UIS), and a continuum of NUISs by different stability boundaries. In the thermodynamic limit, a NUIS with prescribed $\varepsilon$7 loses marginal stability at

$\varepsilon$8

while the UIS corresponds to

$\varepsilon$9

and full synchrony has second-order boundary

$(2+1)$0

These curves are asymptotically tangent at $(2+1)$1 to the exact MF-CGLE bifurcation lines. The second-order model also reproduces QPS-like saddle transients near UIS destabilization, with transient times that scale logarithmically with perturbation amplitude, although stable QPS attractors and pure collective chaos are not recovered within the weak-coupling validity regime (León et al., 2019).

In nonlocally coupled oscillator rings, second-order reduction explains coupling-strength-dependent chimera shape. In the full Stuart–Landau ring, both the coherent-domain size $(2+1)$2 and the coherent-domain frequency $(2+1)$3 vary with $(2+1)$4, and chimera states disappear in a fold bifurcation around $(2+1)$5. The first-order identical-oscillator phase model misses this entirely because its dynamics can be made independent of $(2+1)$6. The second-order phase model reproduces the dependence of $(2+1)$7, $(2+1)$8, phase snapshots, and time-averaged frequency profiles up to moderately strong coupling, roughly $(2+1)$9, and also captures the disappearance of chimera states near the fold (Mau et al., 2023).

For oscillators with phase-dependent amplitude, second-order terms determine where deformation of the limit cycle first enters stability calculations. In the family with $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$00, synchrony stability is independent of $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$01 at $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$02, $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$03, and $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$04 order, and depends on deformation only at $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$05. For $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$06, the corresponding critical Floquet multiplier acquires explicit $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$07-dependence only at that order. In the $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$08 splay-state problem, the $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$09 reduction captures a subcritical Neimark–Sacker bifurcation present in the full system, whereas the first-order reduction does not (Bick et al., 2023).

For driven generic two-dimensional oscillators, higher-order phase reduction improves finite-coupling synchronization predictions even when the problem is not one of collective states. In the harmonically forced van der Pol oscillator, the $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$10 phase model predicts the Arnold tongue with high accuracy, and numerical comparisons show that the second-order approximation remains almost perfectly aligned with the full model even at $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$11, where the first-order approximation has already shifted noticeably. The same framework also predicts higher-order synchronization regions such as $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$12 and $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$13 tongues (Mau et al., 2023).

In delay-coupled Stuart–Landau oscillators, the second-order phase-difference equation

$$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$14

predicts curved delay-dependent stability boundaries and bistability of in-phase and anti-phase synchronization. This bistable region is absent from the first-order phase-lag approximation and was later proved within an arbitrary-order delay-reduction framework for two delay-coupled Stuart–Landau oscillators (Bick et al., 2024, Bick et al., 31 Oct 2025).

A closely related phenomenon appears in the two-group globally coupled Stuart–Landau model analyzed by Mau, Omel'chenko, and Rosenblum. There the second-order Adler equation

$$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$15

produces bistability of in-phase and anti-phase locking near $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$16, where the first-order term vanishes. The same second-order mean-field structure also permits stable $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$17 three-cluster states, which are forbidden in the first-order mean-field description by Watanabe–Strogatz theory (Baibolatov et al., 3 Jun 2026).

6. Validity, limitations, and computational practice

Second-order phase reduction remains a weak-coupling theory. The standard assumptions are a stable hyperbolic limit cycle for each uncoupled oscillator, weak coupling relative to transverse stability, smoothness sufficient for Taylor or Fourier expansion, and persistence of an attracting invariant torus on which amplitudes are smooth functions of phases. In the Stuart–Landau network formulation, this smallness requirement is expressed relative to the Floquet exponents $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$18. In delay systems, the expansion is in $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$19, not in $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$20; delays may be of order one relative to the period, provided the invariant torus persists (Mau et al., 2023, Bick et al., 2024).

Model-specific validity estimates can be quite restrictive. For the MF-CGLE, a conservative estimate proposed in the paper is

$$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$21

Numerically, for $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$22, the second-order approximation remains accurate up to about $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$23, with accuracy improving when $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$24 is smaller. The same study also emphasizes that some phenomena predicted by the reduced model, such as stable clustering at moderate coupling, appear only outside the validity range, while stable QPS attractors and collective chaos are not recovered at second or third order within the weak-coupling regime (León et al., 2019).

The phase-only description can also fail for geometric reasons. In the numerical reconstruction framework, validity requires that the perturbed motion remain on a smooth attracting torus and that the phase based on unperturbed isochrons remain a proper monotone coordinate. The method applies in the quasiperiodic regime, not in synchronization or locking, because the coupling function must be reconstructed from a densely sampled torus. In the driven Rössler example, breakdown occurs when the torus shifts so far relative to the original isochrons that the “good phase condition” $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$25 is violated (Rosenblum et al., 2018).

The numerical approach of Pikovsky and Rosenblum provides a practical route when closed-form higher-order reduction is unavailable. It reconstructs the full effective phase equation

$$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$26

fits $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$27 either on a $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$28 grid or via a finite double Fourier series, and then extracts

$$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$29

from data at multiple coupling strengths. In the Rayleigh oscillator, the full coupling-function model yields a fit error around $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$30 over the explored range, the first-order approximation is accurate only for roughly $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$31, and the third-order approximation reduces the truncation error below $$\dot{\varphi}=\omega+\varepsilon Q_1+\varepsilon^2 Q_2+\cdots,$$32 across the full tested range (Rosenblum et al., 2018).

A final limitation is algebraic complexity. Even for analytically tractable models such as Stuart–Landau oscillators, formulas become cumbersome already at second order and grow rapidly thereafter. This is why the recent literature has concentrated on special classes—exactly reducible Stuart–Landau oscillators, generic two-dimensional oscillators with one transverse mode, or explicit delay constructions—rather than on a single universal higher-order formalism. The cumulative evidence nevertheless suggests a common conclusion: whenever the observed weak-coupling dynamics depends on coupling strength, on delay duration, or on collective states forbidden by purely pairwise first-harmonic phase models, second-order phase reduction is the minimal phase-only description that can recover the missing structure (León et al., 2019, Mau et al., 2023, Bick et al., 2024).