# Modeling dispersive and active ring resonators with discrete optical components: from injected loops to mode-locked lasers

Nicolas Barre, Marc Brunel

► To cite this version:

Nicolas Barre, Marc Brunel. Modeling dispersive and active ring resonators with discrete optical components: from injected loops to mode-locked lasers. Optics Express, 2025, 33 (7), pp.15855-15877. 10.1364/OE.555083 . hal-05105590

HAL Id: hal-05105590

https://hal.science/hal-05105590v1

Submitted on 10 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Modeling dispersive and active ring resonators with discrete optical components: from injected loops to mode-locked lasers

NICOLAS BARRÉ* AND MARC BRUNEL

Univ Rennes, CNRS, Institut Foton UMR 6082, 35000 Rennes, France

*nicolas.barre@univ-rennes.fr

**Abstract:** This article introduces a general framework for modeling the dynamics of active ring resonators with both continuous and discrete optical elements under the traveling-wave approximation. The framework is built upon a rigorous derivation of the two-time-variable formulation, which serves as the basis for several models, including the Haus master equations for lasers and the Lugiato-Lefever equation for microresonators. Starting from the slowly varying envelope approximation, we demonstrate how to integrate discrete optical components into the model, such as modulators, optical couplers, and saturable absorbers, and address gain dynamics using auxiliary differential equations. By explicitly incorporating these discrete elements rather than relying on mean-field approximations, we enhance simulation accuracy without increasing computational complexity. Additionally, we emphasize the crucial role of selecting an appropriate round-trip time, highlighting its importance in avoiding non-physical solutions. Special attention is given to the subtleties of the model, focusing on its effective applications and limitations. The consistency of the model and its improved accuracy compared to mean-field models are supported by numerical simulations at the end of the article, with additional tutorial-like examples provided in the associated supplementary materials.

Published by Optica Publishing Group under the terms of the [Creative Commons Attribution 4.0 License](https://creativecommons.org/licenses/by/4.0/). Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

## 1. Introduction

Optical ring resonators are fundamental components in a broad range of photonic systems, from mode-locked lasers to Kerr microresonators [1] and frequency-shifting loops (FSLs) [2,3]. Modeling their dynamics accurately is crucial for understanding phenomena like pulse formation, soliton dynamics, and modulation effects. Traditionally, this modeling has been done using mean-field approximations, most notably the Haus master equation for lasers [4,5] and the Lugiato-Lefever equation (LLE) [6,7] for Kerr microresonators. These equations are based on scalar and slowly-varying envelope approximations (SVEA), where the dynamics is first-order in the longitudinal spatial coordinate (z-axis) and higher-order in time. By discretizing the field over the time variable and assuming periodic boundary conditions in time, these equations are particularly well-suited to efficient numerical methods like the split-step Fourier method [8]. However, while this approach works well for many systems, it relies on assumptions, such as a well-defined and fixed repetition rate $T_R$, that may not always hold in real-world systems.

In practice, this assumption of time-periodic boundary conditions and the pre-determined knowledge of the repetition rate $T_R$ can limit the applicability of these models. The round-trip time is not always intuitive to define, especially in resonators containing optical elements with their own independent dynamics. Moreover, treating the entire system as a continuous medium, as in the mean-field approach, neglects the finer details that arise when explicitly modeling a sequence of discrete or continuous optical elements, such as gain, modulators, saturable absorbers, or optical couplers, which do not necessarily commute. An alternative approach
is to discretize the field along the propagation direction, where the physical geometry of the resonator as a closed loop naturally supports space-periodic boundary conditions. This is the basis of the finite-difference time-domain (FDTD) family of techniques [9,10]. However, these methods require more complex numerical schemes than the split-step Fourier method, such as Lax-Wendroff, leapfrog, or upwind schemes [11,12]. To ensure stability, these schemes must obey the Courant-Friedrichs-Lewy (CFL) condition, requiring the time step to be small enough so that the wave propagates no more than one spatial grid cell per time step. This constraint often results in extremely small time steps, making simulations impractical for large systems, such as meter-long fiber loops. While finite-difference time-domain (FDTD) methods provide a more rigorous approach and are indispensable when high spatio-temporal resolution is required, their substantial computational cost remains a significant drawback. As a result, FDTD methods are typically reserved for cases where alternative approaches fail to achieve the necessary accuracy or resolution.

We represent these two approaches in Fig. 1, which illustrates two perspectives of a ring resonator. Fig. 1(a) depicts the resonator as a closed loop of length $l$, where the optical field is discretized along the propagation direction $z$, evolving over time $t$. In simple cases, the field primarily undergoes a spatial translation at the group velocity $v_g$, moving through the loop and returning to its starting point after traversing the full length of the resonator, as depicted in Fig. 1(b). This provides a clear physical interpretation of the system dynamics. In contrast, Fig. 1(c) shows the unfolded resonator that is implicitly considered in Haus or LLE models, where the system is represented as an infinite sequence of identical optical systems of length $l$, with the field discretized in time and time-periodic boundary conditions assumed. This unfolded view is typically analyzed in a moving reference frame at the group velocity $v_g$, where the field appears quasi-stationary, as shown in Fig. 1(d). While this perspective simplifies the mathematics, it inherently assumes that the round-trip time and group velocity are constant, which may not always hold, especially in systems with complex or nonlinear dynamics.

Given the widespread use of the Haus master equation and LLE models, along with their ease of implementation, there is a natural tendency to prefer these methods despite their limitations, and to describe the system according to the view in Fig. 1(c). In this article, we derive a rigorous two-time-variable formulation (Section 2) to model active ring resonators that incorporate discrete optical components, such as modulators, couplers, and saturable absorbers. We demonstrate how each optical element can be modeled to enhance the basic ring model, either by adding continuous terms to the propagation equation or by introducing discrete operations on the field at specific locations, depending on the nature of the optical element (Section 3). Additionally, we show how to couple joint differential equations to the main field propagation equation when the inserted optical elements are active and follow their own dynamics (Section 4). We also include a discussion (Section 5) explaining how our model differs from traditional approaches based on the Ikeda map [13], particularly in the presence of modulation, where the latter fails to preserve invariance under reference frame transformations due to the approximations it requires. While the main article focuses on theoretical modeling and a specific numerical case study dedicated to an injected phase-modulated loop (Section 6), the accompanying Supplement 1 provides key examples demonstrating the versatility of this formulation in modeling a wide range of systems. These include active [14,15] and passive mode-locked lasers [16–20] (Sections S1 and S2), as well as Kerr microresonators [21–24] (Section S3). A formal derivation of the LLE based on the Ikeda map is also included (Section S4) to explicitly outline all the underlying hypotheses and approximations, serving as a comparison to the approach we follow in this article. The numerical examples presented in both the main article and the Supplement 1 allow us to discuss the limitations of the assumptions regarding time-periodic boundary conditions, ensuring that readers are aware of potential pitfalls, such as misjudging the round-trip time, which may lead to non-physical drift or spurious solutions. By emphasizing both the limitations and effective
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Two perspectives of a ring resonator with optical coupler (C), gain (G), modulator (M), optical filter (H) and saturable absorber (A). (a) Closed loop representation, where the optical field is discretized along the propagation direction z with space-periodic boundary conditions. (b) Schematic of pulse propagation at the group velocity $v_g$ as described by the model (a). (c) Unfolded resonator model, depicted as an infinite sequence of identical systems, with the field discretized in time and time-periodic boundary conditions assumed. (d) Illustration of a stationary pulse in the co-moving frame at speed $v_g$, with time-periodic boundary conditions.</figcaption></figure>

applications of existing mean-field approaches, we provide a versatile framework for resonator modeling that is computationally efficient while capturing a broad range of physical effects. A key limitation of this work is the traveling-wave assumption, which models only unidirectional field propagation and neglects the effects of spatial hole burning.

## 2. Derivation of the two-time variables formulation

In order to model the evolution of a field $E(z, t)$ under the SVEA, two different sign conventions are often used in the literature to describe the fast variations of the phase:

$$E(z, t) = u(z, t)e^{\sigma i(\omega_0 t - k_0 z)} + \text{c.c.}, \quad (1)$$

with $\sigma = \pm 1$. If one chooses $\sigma = 1$, the complex spectrum $\tilde{E}(z, \omega)$ of the field is obtained by a Fourier transform of $E(z, t)$ with respect to the time variable $t$, whereas if $\sigma = -1$ the complex spectrum must be defined by the inverse Fourier transform of $E(z, t)$. In the following, we do some generic derivations, keeping the sign parameter $\sigma$, and $\tilde{E}(z, \omega)$ refers to either the direct or the inverse Fourier transform of $E(z, t)$, depending on the value of $\sigma$. Conversely, we will represent the inverse transform of the tilde operator with an overline such that $\tilde{E}(z, \omega)(t) = E(z, t)$.
Starting from Maxwell equations and assuming an homogeneous linear susceptibility $\chi$, we easily derive the well-known Helmholtz equation for the Fourier components of $E(z, t)$:

$$ \frac{\partial^2 \tilde{E}(z, \omega)}{\partial z^2} + \frac{\omega^2}{c^2}(1 + \chi(\omega))\tilde{E}(z, \omega) = 0. \qquad (2) $$

Note that this equation does not depend on the sign convention $\sigma$.

According to Eq. (1) we have $\tilde{E}(z, \omega) = \tilde{u}(z, \omega - \omega_0)e^{-\sigma i k_0 z} + \tilde{u}^*(z, \omega + \omega_0)e^{\sigma i k_0 z}$, where the second term possessing higher frequency components can be neglected [25]. By injecting this expression in Eq. (2), and by neglecting the second-order derivative $\partial^2/\partial z^2$, we obtain the SVEA equation:

$$ \frac{\partial \tilde{u}(z, \omega - \omega_0)}{\partial z} = -\sigma \frac{i}{2} \left( \frac{1 + \chi(\omega)}{k_0} \frac{\omega^2}{c^2} - k_0 \right) \tilde{u}(z, \omega - \omega_0). \qquad (3) $$

By defining $\beta(\omega)^2 = (1 + \chi(\omega))\omega^2/c^2 \approx \beta(\omega_0)^2 + 2\beta(\omega_0)\beta(\omega)$ and choosing $k_0 = \beta(\omega_0)$, we get the most common formulation:

$$ \frac{\partial \tilde{u}(z, \omega - \omega_0)}{\partial z} = -\sigma i(\beta(\omega) - \beta(\omega_0))\tilde{u}(z, \omega - \omega_0). \qquad (4) $$

In principle, the choice of the optical carrier pulsation $\omega_0$ can be arbitrary, but still within the range of the frequencies constituting the pulses that are going to be generated. In case of an injected resonator with a monochromatic wave $E_{in}e^{\sigma i \omega_{in}}$, it makes sense to choose $\omega_0 = \omega_{in}$, and in case of a laser containing a gain modeled by a two-level atomic transition of energy $E_{at} = \hbar\omega_{at}$, it makes sense to choose $\omega_0 = \omega_{at}$, although these specific choices are not mandatory and must not influence the computed solution.

We can now make a Taylor expansion of $\beta(\omega)$, leading to

$$ \beta(\omega) - \beta(\omega_0) = (\omega - \omega_0) \frac{\partial \beta}{\partial \omega}(\omega_0) + \beta_R(\omega), \qquad (5) $$

where $\beta_R(\omega)$ contains only higher-order terms ($\ge 2$) in $\omega - \omega_0$. When coming back to the temporal domain, these higher-order terms represent higher-order time derivatives.

By definition, $\frac{\partial \beta}{\partial \omega}(\omega_0)$ represents the inverse of the group velocity $v_g(\omega_0)$ of a wave packet centered around the pulsation $\omega_0$. In the case of a ring resonator of length $\ell$, the group velocity $v_g$ can be used to define the cavity proper time $T_{\Omega}(\omega_0) = \ell/v_g(\omega_0)$, which corresponds to the round-trip time of a pulse centered around the pulsation $\omega_0$ in the passive cavity. As mentioned previously, $\omega_0$ can be chosen arbitrarily, and $T_{\Omega}(\omega_0)$ can differ from the repetition rate of the solution that will be established if additional active elements are inserted into the loop, such as an amplitude or phase modulator, external injection, or a gain medium.

Let us introduce the detuning time $\tau_d(\omega_0) = T_R - T_{\Omega}(\omega_0)$, where $T_R$ is a free parameter, independent of the repetition rate for the moment, to preserve the generality of the following demonstration. Eq. (5) becomes:

$$ \beta(\omega) - \beta(\omega_0) = (\omega - \omega_0) \left( \frac{T_R}{\ell} - \frac{\tau_d(\omega_0)}{\ell} \right) + \beta_R(\omega). \qquad (6) $$

We now substitute the expression (6) into the Eq. (4), leading to

$$ \left( \frac{\partial}{\partial z} + \sigma i(\omega - \omega_0) \frac{T_R}{\ell} \right) \tilde{u}(z, \omega - \omega_0) = \sigma i \left( (\omega - \omega_0) \frac{\tau_d(\omega_0)}{\ell} - \beta_R(\omega) \right) \tilde{u}(z, \omega - \omega_0), \qquad (7) $$

which is equivalent to the following differential equation in the temporal domain:

$$ \left( \frac{\partial}{\partial z} + \frac{T_R}{\ell} \frac{\partial}{\partial t} \right) \left( u(z, t) e^{\sigma i \omega_0 t} \right) = \frac{\tau_d(\omega_0)}{\ell} \frac{\partial}{\partial t} \left( u(z, t) e^{\sigma i \omega_0 t} \right) - \sigma i \left( \overline{\beta_R(\omega)}(t') * u(z, t') e^{\sigma i \omega_0 t'} \right) (t). \qquad (8) $$
Note that the exponential term $e^{\sigma i \omega_0 t}$ can be factored out of the convolution product since

$$ \left( \overline{\beta_R(\omega)(t')} * u(z, t') e^{\sigma i \omega_0 t'} \right) (t) = e^{\sigma i \omega_0 t'} \int \overline{\beta_R(\omega)(t')} e^{-\sigma i \omega_0 t'} u(z, t-t') dt' \\ e^{\sigma i \omega_0 t'} \left( \overline{\beta_R(\omega + \omega_0)(t')} * u(z, t') \right) (t). \qquad (9) $$

By introducing the convolution property (9) into the propagation equation (8) and by deriving the product $u(z, t)e^{\sigma i \omega_0 t}$ with respect to time on both sides, we obtain:

$$ \left( \frac{\partial}{\partial z} + \frac{T_R}{\ell} \frac{\partial}{\partial t} \right) u(z, t) = - \sigma i \frac{\omega_0 T_\Omega(\omega_0)}{\ell} u(z, t) + \frac{\tau_d(\omega_0)}{\ell} \frac{\partial u(z, t)}{\partial t} \\ - \sigma i \left( \overline{\beta_R(\omega + \omega_0)(t')} * u(z, t') \right) (t). \qquad (10) $$

In order to simplify the differential operator on the left hand side of Eq. (10), we introduce the change of variables

$$ T = \frac{T_R}{\ell} z, \quad \tau = t - \frac{T_R}{\ell} z, \qquad (11) $$

such that

$$ z = \frac{\ell}{T_R} T, \quad t = T + \tau, \qquad (12) $$

thereby giving $T_R$ the physical interpretation of the round-trip time of the moving reference frame traveling at a constant velocity $v_R = \ell/T_R$, in which the propagation equation will be now expressed. According to the chain differentiation rule, we have:

$$ \frac{\partial}{\partial T} = \frac{\partial z}{\partial T} \frac{\partial}{\partial z} + \frac{\partial t}{\partial T} \frac{\partial}{\partial t} = \frac{\ell}{T_R} \frac{\partial}{\partial z} + \frac{\partial}{\partial t}, \qquad (13) $$

$$ \frac{\partial}{\partial t} = \frac{\partial T}{\partial t} \frac{\partial}{\partial T} + \frac{\partial \tau}{\partial t} \frac{\partial}{\partial \tau} = \frac{\partial}{\partial \tau}, \qquad (14) $$

which allows to rewrite Eq. (10) under the form:

$$ \frac{T_R}{\ell} \frac{\partial}{\partial T} u \left( \frac{\ell}{T_R} T, T + \tau \right) = - \sigma i \frac{\omega_0 T_\Omega(\omega_0)}{\ell} u \left( \frac{\ell}{T_R} T, T + \tau \right) + \frac{\tau_d(\omega_0)}{\ell} \frac{\partial}{\partial \tau} u \left( \frac{\ell}{T_R} T, T + \tau \right) \\ - \sigma i \left( \overline{\beta_R(\omega + \omega_0)(t')} * u \left( \frac{\ell}{T_R} T, t' \right) \right) (T + \tau). \qquad (15) $$

This motivates the definition of a two-time variable field function $u_{v_R}(T, \tau) = u(\ell T/T_R, T + \tau)$, defined in the frame moving at the velocity $v_R$ in the z direction, which obeys the following evolution equation:

$$ T_R \frac{\partial u_{v_R}(T, \tau)}{\partial T} = - \sigma i \omega_0 T_\Omega(\omega_0) u_{v_R}(T, \tau) + \tau_d(\omega_0) \frac{\partial}{\partial \tau} u_{v_R}(T, \tau) \\ - \sigma i \ell \left( \overline{\beta_R(\omega + \omega_0)(\tau')} * u_{v_R}(T, \tau') \right) (\tau), \qquad (16) $$

where the convolution is performed with respect to the second time variable $\tau$.

Under the SVEA, Eq. (16) is exact, but is not directly usable without further approximations. Indeed, integrating this equation over the slow time variable $T$ requires knowledge of the field $u_{v_R}(T, \tau)$ at all times $\tau$, which is not feasible. To address this, we restrict the domain of $\tau$ to a finite interval $\tau \in [0, \Delta\tau]$ and seek appropriate boundary conditions that minimally influence the solution. If the system exhibits a quasi-periodic behavior, it is natural to set $\Delta\tau$ equal to the
system's repetition rate, as all relevant information is contained within this interval. Furthermore, if the desired solution varies only slightly from one round-trip to the next, assuming periodic boundary conditions is a reasonable choice.

It is crucial to distinguish between assuming periodic boundary conditions and assuming the periodicity of the solution. Before reaching a steady-state, where the solution may become stationary and periodic, periodic boundary conditions do not restrict the evolving field from changing shape or evolving across iterations. However, periodic boundary conditions are not transparent and affect the computation of higher-order derivatives at the edges of the fast time window. Nevertheless, this approach is consistent with Haus's original assumption that the field should not undergo significant variations over a single round-trip. Importantly, while these boundary conditions enforce continuity of the field and its derivatives at the edges of the interval, they do not enforce a stationary or periodic solution. Instead, they allow the system to evolve freely, with the solution naturally converging to a steady-state or periodic behavior over time.

In the literature, a shortcut to such a formulation is typically obtained from the Ikeda map [13], which is a discrete model where the variable $T$ can only take discrete values, $T = n\Delta\tau$ with $n \in \mathbb{N}$, thereby implicitly assuming a moving reference frame with a round-trip time $T_R$ equal to the repetition rate $\Delta\tau$. Then, $T$ is treated as a continuous variable, which gives the impression that an approximation is required. In contrast, our derivation requires no approximation beyond the SVEA to obtain Eq. (16). Moreover, this derivation provides a clear interpretation of each time variable in traditional two-time variable models, which is not as evident in the Ikeda map formulation, where values of $T$ that are not integer multiples of $T_R$ lack a well-defined physical meaning. To better illustrate the conceptual difference between our derivation of Eq. (15) and the additional assumptions required in the Ikeda map formulation, we refer the reader to Section S4 of the Supplement 1.

According to our derivation, the first time variable $T$ is rigorously related to the position $z$ in the loop with a one-to-one mapping $z = v_R T$, depending on the travel speed $v_R = \ell/T_R$ of the reference frame. Obviously, one needs to consider $z$ modulo $\ell$ in order to recover the spatial position of the field $u_{v_R}(T, \tau)$ inside a ring of length $\ell$. $T$ is traditionally called the slow time variable as it usually requires a low sampling in order to integrate Eq. (16) with a good accuracy. Conversely, $\tau$ is called the fast time variable because it is directly related to the true time variable $t = T + \tau$, with both times differing only by the offset $T$.

Above, we stated that the Ikeda map formulation inherently identifies the periodicity of the solution $\Delta\tau$ with the round-trip time $T_R$ of the moving frame. However, in our derivation, $T_R$ can be chosen independently of the solution's period $\Delta\tau$. For simplicity, we also assume throughout this work that the round-trip time of the reference frame is set to $T_R = \Delta\tau$. For completeness, we include an example in Section 6.2 involving a phase-modulated loop, demonstrating that selecting $T_R = T_\Omega$ or $T_R = \Delta\tau$ yields consistent results. Since we did not apply any approximation to simplify the derivative operators, Eq. (16) remains invariant under a change of reference frame. We will explain in the discussion that this is not the case of the master equation derived from the discrete Ikeda map in [5].

Moreover, it is important to note that we do not know the expected periodicity $T_R = \Delta\tau$ a priori. In the following, we will show that in the presence of an active modulation with period $T_M$ inside the loop we can assume $T_R = T_M$. However, since $T_\Omega(\omega_0) \neq T_M$ in general, this leads to a non-zero term $\tau_d(\omega_0)$. In the absence of active modulation, we have no precise knowledge of $T_R$, other than the fact that it should not deviate significantly from $T_\Omega(\omega_0)$. In this case, the initial choice of $T_R$ must be verified a posteriori by checking whether the steady-state solution is self-consistent, as it should be, or whether it drifts after a round-trip.

Finally, it is worth noting that the demonstration we just made in free space by neglecting transverse effects is also valid for waveguides, by using the effective refractive index $n_{\text{eff}}$ in the
definition of $\beta$, given by:

$$ \beta(\omega) = \frac{\omega}{c} n_{\text{eff}}(\omega) = \frac{2\pi}{\lambda} n_{\text{eff}}(\omega). \qquad (17) $$

The main difference in this case is that a mode solver is now required in order to compute the frequency dependence of the effective refractive index $n_{\text{eff}}(\omega)$ associated with a given transverse mode.

### 3. Modeling discrete spatial elements

In this section, we focus on extending the system with optical elements placed at discrete positions inside the loop, such as a light modulator, an optical filter or an optical coupler. We describe how to update the main Eq. (16) with a set of right-hand contributions.

#### 3.1. Modulator

Let us consider the presence of a light modulator (phase or amplitude) at a discrete position $z_s$. Its action on the field $u(z, t)$ can be modeled by the discrete transformation

$$ u(z_s + \frac{\epsilon}{2}, t) = u(z_s - \frac{\epsilon}{2}, t) e^{iS_M(t)}, \qquad (18) $$

where $S_M(t)$ is the modulation signal and $\epsilon$ is a parameter that tends towards zero. For a phase modulation $S_M(t)$ is real, while for an ideal amplitude modulation it is purely imaginary. The continuous differential equation describing accurately the discrete relation (18) reads:

$$ \frac{\partial u(z, t)}{\partial z} = iS_M(t)u(z, t)\delta(z - z_s), \qquad (19) $$

where $\delta$ is the Dirac function. In order to show that Eq. (19) allows to recover the discrete relation (18), we can use a separation of variables in (19) and integrate between $z_s - \epsilon/2$ and $z_s + \epsilon/2$:

$$ \int_{z_s - \epsilon/2}^{z_s + \epsilon/2} \frac{\partial u(z, t)}{\partial z} (u(z, t))^{-1} dz = \left[ \ln(u(z, t)) \right]_{z_s - \epsilon/2}^{z_s + \epsilon/2} = iS_M(t) \int_{z_s - \epsilon/2}^{z_s + \epsilon/2} \delta(z - z_s) dz. \qquad (20) $$

Then, we use the fact that $\forall \epsilon$, $\int_{z_s - \epsilon/2}^{z_s + \epsilon/2} \delta(z - z_s) dz = 1$.

We can now apply the change of variables (12) that we applied earlier, and consider that the modulation signal is applied periodically in $z$ every $z_s + k\ell$, in order to find the equation for $u_{vR}(T, \tau)$:

$$ T_R \frac{\partial u_{vR}(T, \tau)}{\partial T} = i\ell u_{vR}(T, \tau) \sum_{k=0}^{\infty} S_M(T + \tau) \delta \left( \ell \frac{T}{T_R} - z_s - k\ell \right), \qquad (21) $$

where $k \in \mathbb{N}$ is an index counting the number of round-trips made along the ring cavity at the speed $v_R$ of the reference frame.

By introducing the time offset $t_s = T_R z_s / \ell$, we can rewrite Eq. (21) under the more convenient form:

$$ T_R \frac{\partial u_{vR}(T, \tau)}{\partial T} = iT_R u_{vR}(T, \tau) \sum_{k=0}^{\infty} S_M(t_s + kT_R + \tau) \delta(T - t_s - kT_R). \qquad (22) $$

To go further, we can assume that the presence in the loop of a modulator generating a periodic signal of period $T_M$ forces the periodicity of the system to be $\Delta\tau = T_M$, so that the definition domain of $\tau \in [0, T_M]$ can be set with confidence. Then, choosing $T_R = T_M$ for the reference frame round-trip time simplifies the modulation function in Eq. (22) to $S_M(t_s + kT_R + \tau) = S_M(t_s + \tau)$, making it independent of the discrete index $k$. However, in Section 6.2, we demonstrate that the
choice of $T_R = T_M$ is not mandatory. For instance, choosing $T_R = T_Ω$ would also yield consistent results, although it requires keeping track of the round-trip index $k$ and translating the modulation function by an additional $T_R$ at each iteration. We refer the reader to Section 6.3, where it is shown that choosing a different domain for the definition of $τ$ in the presence of modulation leads to spurious artifacts in the computed solution.

To facilitate a comparison between the exact discrete model and a mean-field approximation, such as those commonly favored in the physics community, we average each Dirac delta function appearing in (22) over a period $T_R$. Although this approach is not strictly rigorous from a mathematical perspective, it provides the most consistent mean-field formulation in comparison to existing literature, particularly in [5]:

$$T_R \frac{\partial u_{v_R}(T, \tau)}{\partial T} = iS_M(t_s + kT_R + \tau)u_{v_R}(T, \tau), \quad \text{for each roundtrip } k, \qquad (23)$$

which, in the particular case $T_R = T_M$, simplifies to:

$$T_M \frac{\partial u_{v_M}(T, \tau)}{\partial T} = iS_M(t_s + \tau)u_{v_M}(T, \tau), \qquad (24)$$

where $v_M = \ell/T_M$ denotes the speed of the reference frame synchronized with the modulation period $T_M$.

## 3.2. Optical filter

Let us consider an optical filter with transfer function $H(\omega + \omega_0)$ placed at a discrete position $z_f$. Its action on the Fourier transform of the slowly-varying envelope $\tilde{u}(z, \omega)$ is given by:

$$\tilde{u}\left(z_f + \frac{\epsilon}{2}, \omega\right) = \tilde{u}\left(z_f - \frac{\epsilon}{2}, \omega\right) H(\omega + \omega_0), \qquad (25)$$

where $\epsilon$ is again a parameter as small as possible, and $\omega$ represents the relative pulsation around the carrier $\omega_0$.

Similarly to what we have done previously, we can build a continuous differential equation for $\tilde{u}(z, \omega)$ that describes the discrete relation (25):

$$\frac{\partial \tilde{u}(z, \omega)}{\partial z} = \ln(H(\omega + \omega_0))\tilde{u}(z, \omega)\delta(z - z_f). \qquad (26)$$

In the reference frame moving at the speed $v_R$ around the ring, taking into account the contributions at every round-trip, we get:

$$T_R \frac{\partial \tilde{u}_{v_R}(T, \omega)}{\partial T} = T_R \ln(H(\omega + \omega_0)) \tilde{u}_{v_R}(T, \omega) \sum_{k=0}^{\infty} \delta(T - t_f - kT_R), \qquad (27)$$

where $t_f = T_R z_f / \ell$.

We can also make a mean-field approximation by averaging every delta function over a period $T_R$, leading to the simpler form:

$$T_R \frac{\partial \tilde{u}_{v_R}(T, \omega)}{\partial T} = \ln(H(\omega + \omega_0)) \tilde{u}_{v_R}(T, \omega). \qquad (28)$$

We note that equations (27) and (28) are expressed in the Fourier domain, while the main Eq. (16) is expressed in the time domain. Numerically, this is not a problem since it is common to solve them together with a split-step Fourier scheme. However, it is always possible to do a Taylor expansion of the logarithm of the transfer function and go back into the time domain with
higher order time derivative operators. It is also interesting to note that the filter transfer function plays the same role as $\beta_R(\omega + \omega_0)$ in Eq. (16). Adding a filter can thus be seen as engineering the gain and dispersion properties of the propagation medium. With this in mind, it may be useful to subtract the first-order term in $\omega$ from the imaginary part of the logarithm of the filter transfer function $\ln(H(\omega + \omega_0))$, and include its contribution in $T_\Omega(\omega_0)$ instead.

Regarding the modeling of intracavity filters, we must make a final remark regarding the convention used for the Fourier transform, with $\sigma = \pm 1$. When $\sigma = 1$ and the complex spectrum is defined as the Fourier transform of the field, it is known in signal processing that the real and imaginary parts of a causal filter transfer function $H(\omega)$ must satisfy Kramers-Kronig relations. However, when $\sigma = -1$, the filter transfer function should be conjugated compared to the other convention, and its real and imaginary parts must follow a modified Kramers-Kronig relation, with a sign-swap. Such subtleties are not much discussed in the literature and can easily lead to sign errors in numerical simulations. It is also important to note that most papers and books dealing with optical pulse propagation use the convention $\sigma = -1$, which dictates the use of ifft instead of fft to calculate numerical optical spectra.

### 3.3. Optical coupler

We now want to model an optical coupler inserted at a position $z_c$, with a transfer matrix $T_c$ defined by:

$$T_c = \begin{pmatrix} t_{11} & t_{12} \\ t_{21} & t_{22} \end{pmatrix}. \qquad (29)$$

The first input is an injection field of the form $u_{\text{in}}(t)e^{\sigma i\omega_0 t}$ and the first output is a field of the form $u_{\text{out}}(t)e^{\sigma i\omega_0 t}$. The second input is an end point of the loop and the second output is an entry point of the loop. The different inputs and outputs are related by the discrete relations:

$$u_{\text{out}}(t) = t_{11}u_{\text{in}}(t) + t_{12}u\left(z_c - \frac{\epsilon}{2}, t\right), \qquad (30)$$

$$u\left(z_c + \frac{\epsilon}{2}, t\right) = t_{21}u_{\text{in}}(t) + t_{22}u\left(z_c - \frac{\epsilon}{2}, t\right), \qquad (31)$$

where $\epsilon$ is again a small parameter that tends towards zero.

Constructing a continuous differential equation that describes accurately the relation (31) is more subtle than in the previous cases. To illustrate this, we can show that if the injection of $u_{\text{in}}(t)$ and the attenuation of $u(z, t)$ by a factor $t_{22}$ take place at two slightly different positions, namely $z_c$ and $z'_c$, the continuous equivalent equation reads:

$$\frac{\partial u(z, t)}{\partial z} = t_{21}u_{\text{in}}(t)\delta(z - z_c) + \ln t_{22} u(z, t)\delta(z - z'_c), \qquad (32)$$

whereas if they occur at the same position $z_c$, a proper integration leads to:

$$\frac{\partial u(z, t)}{\partial z} = \left( t_{21} \frac{\ln t_{22}}{t_{22} - 1} u_{\text{in}}(t) + \ln t_{22} u(z, t) \right) \delta(z - z_c), \qquad (33)$$

where a correction factor appears in front of the injection signal $u_{\text{in}}(t)$.

This highlights the inconsistency that arises when attempting to construct a mean-field model with an injected field. Although an approximate formulation, the LLE mean-field equation with a constant drive term $t_{21}u_{\text{in}}$ serves as the foundation for many simulation codes [26,27] within the microresonator community. In the absence of an optical modulator in the loop, we verified that both formulations, as given by equations (32) and (33), closely approximate the exact solution
obtained by treating the injection as a discrete event. However, when modulation is present, the more accurate mean-field formulation is expressed as:

$$T_R \frac{\partial u_{v_R}(T, \tau)}{\partial T} = t_{21} \frac{\ln t_{22} + iS_M(t_s + \tau) - \sigma i\omega_0 T_{\Omega}(\omega_0)}{t_{22}e^{iS_M(t_s+\tau)-\sigma i\omega_0 T_{\Omega}(\omega_0)} - 1} u_{\text{in}} + \ln t_{22} u_{v_R}(T, \tau), \quad (34)$$

which derivation is provided in Section S4.B of the Supplement 1. In the presence of modulation, the correction factor preceding the drive field $u_{\text{in}}$ becomes significant and leads to incorrect solutions if it is not accounted for, as illustrated later in the example of Section 6.4. Moreover, if the drive field $u_{\text{in}}(t)$ varies over time, no exact analytical expression exists for the correction factor, particularly when a nonlinear Kerr effect is included.

In summary, the use of a mean-field formulation should be strictly avoided when dealing with a drive field. While such formulations are highly useful for deriving analytical solutions and performing stability analyses, they offer no computational advantage from a purely numerical perspective and can lead to inaccuracies. Instead, the injection field should be treated as a discrete event by halting the integration of the main field evolution Eq. (16) at discrete times $T_k$ corresponding to the discrete position $z_c$ of the optical coupler in the loop, and simply adding the entire drive field at this position.

Numerically, we implement this by interrupting the continuous integration described previously at discrete times $T_k$ and by applying the following transfer function:

$$u^+(T_k, \tau) = t_{22}u^-(T_k, \tau) + t_{21}u_{\text{in}}, \quad (35)$$

where $t_{21}$ and $t_{22}$ are coupling coefficients of the transfer matrix introduced in Section 3.3, and $u^-$ and $u^+$ refer to the field distribution $u$ just before and after the injection at $T = T_k$, respectively.

## 4. Modeling saturation dynamics

In this section, we describe how to extend the system by incorporating dynamical elements governed by their own set of differential equations into the resonator. This involves coupling the field propagation equation with additional time-dependent evolution equations. To illustrate this approach, we focus on modeling a saturated medium, such as gain, distributed across a portion or the entire length of the cavity. Furthermore, we demonstrate how this method can be applied to discrete elements, such as a saturable absorber placed inside the cavity. It is important to note that the approach followed here is not limited to modeling saturation dynamics, and can be generalized to other differential equations coupled with the system.

The amplification process originates from the material polarization defined as $P(z, t) = p(z, t)e^{\sigma i(\omega_0 t - k_0 z)} + \text{c.c.}$, which contributes a new source term to Eq. (10):

$$\left(\frac{\partial}{\partial z} + \frac{T_R}{\ell} \frac{\partial}{\partial t}\right) u(z, t) = -\sigma \frac{i}{2} \frac{\omega_0}{\epsilon_0 n(\omega_0) c} p(z, t). \quad (36)$$

In the simplest two-level atom model, the evolution of the polarization $p(z, t)$ and the density of the population inversion $\Delta n(z, t)$ follow the Maxwell-Bloch equations [28]:

$$\frac{\partial p(z, t)}{\partial t} = -(\Omega_G + \sigma i \delta_G) p(z, t) + \sigma i \frac{d^2}{\hbar} \Delta n(z, t) u(z, t), \quad (37)$$

$$\frac{\partial \Delta n(z, t)}{\partial t} = \frac{1}{T_G}(\Delta n_0(z) - \Delta n(z, t)) - \frac{2}{\hbar} \Im(\sigma u^*(z, t) p(z, t)), \quad (38)$$

where $\Im$ denotes the imaginary part, $\Omega_G$ is the gain bandwidth of the Lorentzian atomic line, $\delta_G = \omega_0 - \omega_{at}$ is the detuning of the carrier pulsation with respect to the atomic transition pulsation, $T_G$ is the upper-state transition lifetime, and $d$ is the elementary dipolar moment of the
material. $\Delta n_0(z)$ is a source term for the density of the population inversion that depends on the pumping rate. Although the two-level atom model is simplified compared to multilevel doped crystals, it is often a good approximation, for a given pump rate, to reduce the multilevel model to a two-level system with an effective pump rate and an effective lifetime $T_G$ for the upper state of the transition [29,30].

When performing the variable change (12), according to the chain rule relation (14), we have:

$$ \frac{\partial}{\partial t} p(z, t) = \frac{\partial}{\partial \tau} p\left(\frac{\ell}{T_R} T, T + \tau\right), \qquad (39) $$

and

$$ \frac{\partial}{\partial t} \Delta n(z, t) = \frac{\partial}{\partial \tau} \Delta n \left( \frac{\ell}{T_R} T, T + \tau \right). \qquad (40) $$

In order to simplify the notations and to rewrite the equations (37) and (38) taking into account this variable change, we define two functions in the new reference frame moving at the speed $v_R$:

$$ s_{v_R}(T, \tau) = -\sigma i \frac{\omega_0}{\epsilon_0 n(\omega_0) c} p\left(\frac{\ell}{T_R} T, T + \tau\right), \qquad (41) $$

$$ g_{v_R}(T, \tau) = \frac{1}{\Omega_G \epsilon_0 n(\omega_0) c} \frac{d^2}{\hbar} \Delta n \left( \frac{\ell}{T_R} T, T + \tau \right). \qquad (42) $$

This gives us the new set of equations:

$$ T_R \frac{\partial u_{v_R}(T, \tau)}{\partial T} = \frac{\ell}{2} s_{v_R}(T, \tau), \qquad (43) $$

$$ \frac{\partial s_{v_R}(T, \tau)}{\partial \tau} = -(\Omega_G + \sigma i \delta_G) s_{v_R}(T, \tau) + \Omega_G g_{v_R}(T, \tau) u_{v_R}(T, \tau), \qquad (44) $$

$$ T_G \frac{\partial g_{v_R}(T, \tau)}{\partial \tau} = g_0(T) - g_{v_R}(T, \tau) - \frac{2d^2 T_G}{\hbar^2 \Omega_G} \mathcal{R}(u_{v_R}^*(T, \tau) s_{v_R}(T, \tau)), \qquad (45) $$

where $\mathcal{R}$ denotes the real part.

### 4.1. Slow gain approximation

When modeling a gain medium, a common approximation is to assume that the population inversion, described by the function $g_{v_R}(T, \tau)$, does not vary significantly over the round-trip time $T_R$. This approximation holds when $T_G \gg T_R$, and involves replacing the term $g_{v_R}(T, \tau)$ by its mean value over a round-trip $g_m(T) = \langle g_{v_R}(T, \tau) \rangle_{T_R}$ in Eq. (44). Doing so, we can establish a Fourier relation between the field $u_{v_R}(T, \tau)$ and its source term $s_{v_R}(T, \tau)$:

$$ \tilde{s}_{v_R}(T, \omega) = g_m(T) \mathcal{L}(\omega + \delta_G) \tilde{u}_{v_R}(T, \omega), \qquad (46) $$

where $\mathcal{L}(\omega)$ is a complex Lorentzian profile defined by:

$$ \mathcal{L}(\omega) = \frac{\Omega_G}{\Omega_G + \sigma i \omega}. \qquad (47) $$

It is now possible to get rid of Eq. (44) by injecting the relation (46) in the equations (43) and (45). After averaging Eq. (45) over $\tau \in [0, T_R]$ and assuming $\frac{g_{v_R}(T, T_R) - g_{v_R}(T, 0)}{T_R} \approx \frac{\partial g_m(T)}{\partial T}$, we
obtain the new set of coupled equations:

$$T_R \frac{\partial \tilde{u}_{vR}(T, \omega)}{\partial T} = \frac{\ell}{2} g_m(T) \mathcal{L}(\omega + \delta_G) \tilde{u}_{vR}(T, \omega), \quad (48)$$

$$T_G \frac{\partial g_m(T)}{\partial T} = g_0(T) - g_m(T) \left( 1 + \frac{1}{T_R} \int \mathfrak{R}(\mathcal{L}(\omega + \delta_G)) \frac{|\tilde{u}_{vR}(T, \omega)|^2}{I_{\text{sat}}} d\omega \right), \quad (49)$$

where $I_{\text{sat}}$ is the saturation intensity defined as

$$I_{\text{sat}} = \frac{\hbar^2 \Omega_G}{2d^2 T_G}, \quad (50)$$

with dimensions of $V^2 \cdot m^{-2}$ in SI units. This is consistent with the squared modulus of an electric field but differs from the conventional units of intensity, $W \cdot m^{-2}$. For numerical applications, since Eqs. (48), (49) no longer depend on fundamental physical constants, the saturation intensity can be expressed in arbitrary units of $|u|^2$, which saves from multiplying the squared fields everywhere by the constant factor $\epsilon_0 c / 2$. Moreover, since $u$ does not depend on the transverse dimensions in most simple models, it is often practical to interpret $|u|^2$ directly as electric power in watts (W) rather than as intensity in $W \cdot m^{-2}$. In this case, $I_{\text{sat}}$ should be interpreted as saturation power, also in watts, and must be determined experimentally.

Eq. (48) expresses the fact that the Lorentzian gain profile acts as a bandpass filter on the field $u_{vR}(T, \tau)$. This can be solved in combination with the main field evolution Eq. (16) by a split-step Fourier method, while simultaneously addressing the dispersion term $\beta(\omega)$.

In many situations, the mean-gain evolution equation (49) can be further simplified by neglecting the Lorentzian profile in the integral, which, by application of Parseval's theorem, results in:

$$T_G \frac{\partial g_m(T)}{\partial T} = g_0(T) - g_m(T) \left( 1 + \frac{I_m(T)}{I_{\text{sat}}} \right), \quad (51)$$

where $I_m(T) = \langle |u_{vR}(T, \tau)|^2 \rangle_{T_R}$ is the average intensity over a round-trip.

Concerning the numerical simulation of Eq. (49), the numerical scheme used to perform one time-step evolution of the gain $g_m$ is as follows:

$$g_m(T + dT) = g_{\text{sat}}(T) + (g_m(T) - g_{\text{sat}}(T)) \exp \left\{ - \left( 1 + \frac{I_{\text{avg}}(T)}{I_{\text{sat}}} \right) \frac{dT}{T_G} \right\}, \quad (52)$$

with

$$g_{\text{sat}}(T) = \frac{g_0(T)}{1 + \frac{I_{\text{avg}}(T)}{I_{\text{sat}}}}, \quad (53)$$

$$I_{\text{avg}}(T) = \frac{1}{T_R} \int \mathfrak{R}(\mathcal{L}(\omega + \delta_G)) |\tilde{u}(T, \omega)|^2 d\omega, \quad (54)$$

where the average intensity $I_{\text{avg}}$ can be replaced by $I_m = \langle |u|^2 \rangle_{T_R}$ if the simplified equation (51) is used instead of (49).

## 4.2. Fast dipole approximation: adiabatic elimination

Another approximation involves considering that the atomic response to the field variations is so fast that the dipoles reach their stationary state instantaneously. This leads to the adiabatic elimination of the polarization by setting $\frac{\partial s_{vR}(T, \tau)}{\partial \tau} = 0$ in Eq. (44). This assumption holds if the inverse of $\Omega_G$ is much smaller than both the round-trip time $T_R$ and the population inversion
relaxation time $T_G$. In this case, Eq. (44) yields a simple relation in the time domain between the field $u_{vR}(T, \tau)$ and its source term $s_{vR}(T, \tau)$:

$$s_{vR}(T, \tau) = g_{vR}(T, \tau) \mathcal{L}(\delta_G) u_{vR}(T, \tau). \qquad (55)$$

By substituting the relation (55) into the equations (43) and (45), we obtain the new set of coupled evolution equations:

$$T_R \frac{\partial u_{vR}(T, \tau)}{\partial T} = \frac{\ell}{2} g_{vR}(T, \tau) \mathcal{L}(\delta_G) u_{vR}(T, \tau), \qquad (56)$$

$$T_G \frac{\partial g_{vR}(T, \tau)}{\partial \tau} = g_0(T) - g_{vR}(T, \tau) \left( 1 + \mathfrak{R}(\mathcal{L}(\delta_G)) \frac{|u_{vR}(T, \tau)|^2}{I_{\text{sat}}} \right), \qquad (57)$$

where $I_{\text{sat}}$ is again given by Eq. (50). Under the fast dipole approximation, we observe that the gain no longer acts as a bandpass filter and is instead saturated by the instantaneous intensity rather than the average intensity over a round-trip. This is not a good approximation for slow gain media, but it is well-suited to describe saturable absorbers instead.

In contrast to the slow gain approximation case, Eq. (57) still depends on the fast time variable $\tau$ and necessitates a proper integration scheme in order to be solved alongside the field evolution equation (56). The initial conditions for the set of equations (56), (57) are $u_{vR}(0, \tau)$ for $\tau \in [0, T_R]$ and $g_{vR}(T, 0)$ for $T \in [0, T_R]$. At a given time $T < T_R$, in order to do an integration step of Eq. (56), we need to know $g_{vR}(T, \tau)$ for all $\tau \in [0, T_R]$. For this, we need to integrate Eq. (57) over an interval of length $T_R$, starting from $\tau = 0$ with $g_{vR}(T, 0)$ as initial condition. If $T \ge T_R$, there is no more initial condition predefined for $g_{vR}(T, 0)$. In this case, we have to refer to the value computed at the last round-trip at $T - T_R$, corresponding to the same relative position $z$ in the ring. In order to do this correctly, we should note that by the definition (42) we have $g_{vR}(T, \tau) = g_{vR}(T - T_R, \tau + T_R) = g_{vR}(T - T_R, \tau)$. Therefore, this time the already calculated value $g_{vR}(T - T_R, 0)$ should be used as the initial condition to integrate (57) over a new interval of length $T_R$.

If neither the slow gain approximation nor the fast dipole approximation hold, the same integration scheme can be applied to the full set of equations (43–45).

### 4.3. Saturable absorber modeling

In order to model the dynamics of a saturable absorber, we can rely on the fast dipole approximation and the equations (56) and (57) that we just derived, with the assumption that the detuning $\delta_g = 0$. The only difference here is that the saturable absorber is considered as a discrete element positioned at $z = 0$, leading to a discrete formulation of Eq. (56), as follows:

$$u^+(T_k, \tau) = u^-(T_k, \tau) e^{-\frac{1}{2}a(T_k, \tau)}, \qquad (58)$$

$$T_A \frac{\partial a(T_k, \tau)}{\partial \tau} = a_0 - a(T_k, \tau) \left( 1 + \frac{|u^-(T_k, \tau)|^2}{I_{\text{Asat}}} \right), \qquad (59)$$

where $T_A$ is the relaxation time, $a_0$ is the unsaturated absorption and $I_{\text{Asat}}$ is the saturation intensity of the absorber. Similarly to the discrete injection case, the system of equations (58) and (59) needs to be solved only at discrete times $T_k = kT_R$.
In the particular case of a fast saturable absorber, we can perform an adiabatic elimination of Eq. (59) and obtain a single field evolution equation:

$$u^{+}(T_k, \tau) = u^{-}(T_k, \tau) \exp\left\{-\frac{1}{2}a_{\text{sat}}(T_k, \tau)\right\}, \qquad (60)$$

with

$$a_{\text{sat}}(T_k, \tau) = \frac{a_0}{1 + \frac{|u^{-}(T_k, \tau)|^2}{I_{\text{Asat}}}}. \qquad (61)$$

Otherwise, we need to solve Eq. (59) at every $T_k$ for $\tau \in [0, T_R]$. To do so, we can employ the same numerical scheme as for Eq. (52):

$$a(T_k, \tau + d\tau) = a_{\text{sat}}(T_k, \tau) + (a(T_k, \tau) - a_{\text{sat}}(T_k, \tau)) \exp\left\{-\left(1 + \frac{|u^{-}(T_k, \tau)|^2}{I_{\text{Asat}}}\right) \frac{d\tau}{T_A}\right\}, \qquad (62)$$

where the initial value $a(T_k, 0)$ at each $T_k$ should be set to the last value $a(T_{k-1}, T_R)$ from the previous round-trip at $T_{k-1}$.

## 4.4. Inhomogeneous gain

Under the slow gain approximation, it is straightforward to generalize Eqs. (48), (49) to account for inhomogeneous gain. Here, we consider an inhomogeneous gain characterized by a continuous set of atomic transitions $\omega_{at}$, but it can be easily adapted to account for a discrete set of dipole frequencies.

To do so, we introduce a new parameter $\omega' = \omega_{at} - \omega_0 = -\delta_G$. Eq. (49) simply needs to be rewritten to explicitly show the dependence of $g_m(T)$ on $\omega'$, which was previously hidden in the fixed parameter $\delta_G$. Then, Eq. (48) can be extended by integrating over $\omega'$, considering all the atomic contributions. This leads to the new set of equations:

$$T_R \frac{\partial \tilde{u}_{v_R}(T, \omega)}{\partial T} = \frac{\ell}{2} \left( \int g_m(T, \omega') \mathcal{L}(\omega - \omega') d\omega' \right) \tilde{u}_{v_R}(T, \omega), \qquad (63)$$

$$T_G \frac{\partial g_m(T, \omega')}{\partial T} = g_0(T, \omega') - g_m(T, \omega') \left(1 + \frac{1}{T_R} \int \mathfrak{R}(\mathcal{L}(\omega - \omega')) \frac{|\tilde{u}_{v_R}(T, \omega)|^2}{I_{\text{sat}}} d\omega'\right). \qquad (64)$$

It is important to note that the function $g_0(T, \omega')$ cannot be directly measured experimentally. However, at very low power, Eq. (63) allows for the measurement of the gain spectrum, which is a convolution of $g_0(T, \omega')$ with the complex Lorentzian profile. The gain profile $g_0(T, \omega')$ can then be reconstructed by a deconvolution method, assuming we know the Lorentzian linewidth $\Omega_G$, which is assumed to be the same for all the considered atomic transitions.

## 4.5. Noise modeling

Finally, accurate noise modeling is crucial for realistic simulations of laser systems, particularly to capture the effects of amplified spontaneous emission (ASE). ASE noise arises from spontaneous photon emissions within the gain medium, introducing randomness in the laser output. Properly accounting for ASE ensures the initiation of laser oscillation, especially when the intra-cavity field is initially zero.

ASE noise is effectively characterized by a complex Gaussian distribution, suitable for modeling thermal noise processes. The real $u_r$ and imaginary $u_i$ components of the complex noise signal $u_{\text{noise}} = u_r + i u_i$ are independent Gaussian random variables with zero mean and variance $\sigma_{\text{noise}}^2$. The relationship between the noise power $P_{\text{noise}}$ and the variance is given by $P_{\text{noise}} = 2\sigma_{\text{noise}}^2$.
To effectively introduce the Gaussian noise field $u_{\text{noise}}$ into the simulation, care must be taken to maintain consistency. A reliable approach is to add a newly generated random field $u_{\text{noise}}(kT_R, \tau)$ to the main field $u(kT_R, \tau)$ after each round-trip time $T_R$. This ensures that the noise power injected within a round-trip does not depend on the discretization of the slow time variable $T$.

Instead, one might consider adding a new stochastic source term in Eq. (16), simplified as:

$$T_R \frac{\partial u}{\partial T} = u_{\text{noise}}(T, \tau), \qquad (65)$$

where $u_{\text{noise}}(T, \tau)$ represents a newly generated random noise distribution at each slow time $T$.

To integrate Eq. (65) such that the total total noise power remains independent of the number of steps $N$, a suitable integration scheme like the Euler-Maruyama method [31] should be employed. This leads to the following evolution rule:

$$u(T + dT, \tau) = u(T, \tau) + \sqrt{\frac{dT}{T_R}} u_{\text{noise}}(T, \tau), \qquad (66)$$

which yields a physically correct result.

## 5. Discussion

In the presence of modulation, to model a resonator incorporating dispersion, gain, and Kerr nonlinearity, we construct a composite field propagation equation expressed in the reference frame moving at the velocity $v_M = \ell/T_M$. This equation is derived from the set of equations (16), (24), (28) and (48), with the inclusion of a nonlinear Kerr effect characterized by a nonlinear coefficient $\gamma$ [8,25]:

$$T_M \frac{\partial u}{\partial T} = \left( \frac{\ell}{2} \alpha - \sigma i \delta_0 + i S_M - \sigma i \ell \gamma |u|^2 \right) u + \frac{\left( \sigma i (\omega (T_M - T_\Omega) - \ell \beta_R) + \ln H_R + \frac{\ell}{2} g_m \mathcal{L}_R \right) \tilde{u}}{2 g_m \mathcal{L}_R}, \qquad (67)$$

where we have omitted the function parameters for the sake of conciseness. We have introduced a linear loss parameter $\alpha<0$, assuming that the losses are independent of frequency, which allows us to consider $\beta_R(\omega + \omega_0)$ as a purely real function. The functions $H_R(\omega)$ and $\mathcal{L}_R(\omega)$ differ from the original filter transfer function $H(\omega)$ and the complex Lorentzian profile $\mathcal{L}(\omega)$ in that we have subtracted the constant phase and group velocity contributions at $\omega = \omega_0$. Instead, these contributions should be incorporated into the time constant $T_\Omega(\omega_0) = \ell/v_g(\omega_0)$, as they directly contribute to the group velocity of the cold cavity. More explicitly, we have:

$$H_R(\omega) = H(\omega) \exp\{-i \arg(H(\omega_0))\} \exp\left\{\left(-i(\omega - \omega_0)\frac{\partial \arg(H)}{\partial \omega}(\omega_0)\right)\right\}, \qquad (68)$$

$$\mathcal{L}_R(\omega) = \mathcal{L}(\omega) - i\mathfrak{F}(\mathcal{L}(\omega_0)) - i(\omega - \omega_0)\frac{\partial \mathfrak{F}(\mathcal{L})}{\partial \omega}(\omega_0). \qquad (69)$$

It is worth noting that Perego et al. [5] recently derived a master equation for laser mode-locking very similar to Eq. (67), starting from the Ikeda map formulation and then applying a change of variables. In doing so, they had to simplify the equation and neglect certain derivatives to obtain their result. In contrast, we did not perform any additional analytical approximations when deriving Eq. (67). Without detuning, both our approach and the Ikeda map formulation surprisingly lead to the same master equation, suggesting that the additional approximations in the latter case did not have significant consequences. However, in the presence of modulation
and detuning, our result differs slightly from the master equation obtained in [5]. Indeed, in the reference frame moving at the velocity $v_M$, the time prefactor of the field derivative on the left-hand side of Eq. (67) is $T_M$ in our case, whereas the time prefactor predicted by Perego et al. [5] in their master equation corresponds to $T_Ω$ in our notation (we use $T_Ω$ where they use $T_R$). This difference, which looks minor at first glance, comes from the simplifications applied to the derivative operators after a change of variables was introduced on top of the discrete Ikeda map formulation. However, such simplifications result in an equation that is no longer invariant under a change of reference frame, which can lead to noticeable differences for large detuning. In this case, an appropriate experiment could help distinguish between the two formulations, which we reserve for future investigation. Despite this discrepancy, the critical detuning term proportional to $τ_d = T_M - T_Ω$ on the right-hand side remains the same in both formulations.

To solve Eq. (67), we divide the round-trip time $T_R = T_M$ into $N$ steps $dT$ such that $T_M = NdT$, and employ a classical split-step Fourier method at each time step. Most of the time, when the gain is not too strong and the nonlinear coefficient $γ$ is zero or sufficiently small, we have observed that $N = 1$ is sufficient. Furthermore, Eq. (67) needs to be solved together with the gain evolution equation (49).

To incorporate discrete elements, such as injection through an optical coupler or a saturable absorber, the numerical integration must be paused after each round-trip. At this point, the discrete field evolution operations (35) and (60) are applied to the field $u_{v_M}(T, τ)$. Two exceptions were made for the modulation function and the optical filter, which are also discrete elements but were incorporated into the continuous differential equation (67) using their mean-field formulations. This choice was made for the sake of elegance in the formulation, and no reduction in accuracy was observed in the simulations that we performed due to this treatment.

Finally, we offer recommendations regarding the choice of the periodicity of the solution $Δτ$, as selecting an incorrect value can have significant implications, including the establishment of non-physical solutions.

In the presence of modulation, $Δτ$ should be set to the modulation period $T_M$ for consistency. This represents a typical case where a stationary solution appears to propagate at a speed different from the group velocity $v_g$ over a round-trip, as $T_M$ can differ from $T_Ω$ in general. Examples of such scenarios are provided in Section S2 of the Supplement 1 for active mode-locking and in the next Section 6.1 for the simple phase-modulated loop case. Moreover, we demonstrate in the numerical example of Section 6.3 that an incorrect choice of the period $Δτ$ in the presence of modulation results in spurious artifacts.

In the absence of modulation, the logical choice for the period is $Δτ = T_Ω$. A slightly incorrect value of $Δτ$ can be identified by observing a drifting pulse after each round-trip. This drift can be corrected by adjusting $Δτ$ so that the detuning term $τ_d(ω_0)$ in Eq. (16) precisely cancels the observed drift. When setting a new value for $Δτ$, the algorithm should be restarted from the beginning, with the expectation that it will follow the same convergence path as before and eventually reach a self-consistent solution.

If achieving perfect self-consistency proves difficult, one may have to settle for a solution that drifts slightly after each round-trip but is still considered sufficiently accurate. Alternatively, one could abandon the two-time-variable model and instead use a finite difference time domain (FDTD) method, which more accurately handles system geometry and boundary conditions. However, this approach should be viewed as a last resort, as FDTD simulations tend to be significantly more resource-intensive in terms of both memory consumption and computational time.

## 6. Numerical simulation of an injected phase-modulated loop

Frequency-Shifting Loops (FSLs) are optical resonators that enable high repetition rate pulse generation through precise frequency shifting. Typically, an FSL consists of a continuous-wave
laser injected into a loop containing a modulator, an amplifier, and a filter. While early FSL systems used acousto-optic frequency-shifters [2,3], recent implementations have incorporated electro-optic modulators [32,33], offering greater bandwidth and tunability.

FSLs are employed in applications such as real-time Fourier transformation, high-speed optical communication, and pulse generation. The system exploits the temporal Talbot effect to convert a continuous-wave laser into Fourier-transform-limited pulses, enabling flexible pulse repetition rates tied to the cavity round-trip frequency. These features make FSLs valuable for advanced optical signal processing.

In this section, we will simplify the example as much as possible to focus on the subtleties of the model, particularly regarding modulation detuning with respect to the cavity proper time $T_Ω(ω_0)$, the effect of injection on a pulse, the choice of the reference frame for tracking the traveling pulse, the selection of the fast time window width, and the breakdown of the traditional mean-field model under large modulation frequency detuning. To do so, we will neglect dispersion and nonlinear effects and exclude gain from the simulation, thereby considering a simple phase-modulated loop.

We consider a ring resonator with a length of $ℓ = 10$ m injected by a continuous-wave at the pulsation $ω_0$, with an intensity $I_{in} = 1$ (expressed in units of $|u|^2$), and phase-modulated by a sinusoidal signal $S_M(t) = Γ_m \sin(2πf_m t)$ at the frequency $f_m = 1/T_M$, where $T_M$ is the modulation period, and the modulation depth $Γ_m = 5$. We set the coupling coefficient of the external field to $t_{21} = \sqrt{0.05}$ and the transmission coefficient of the loop after a round-trip to $t_{22} = \sqrt{0.95}$. We assume that all losses are accounted for by the transmission coefficient $t_{22}$, and we neglect the linear losses of the loop. The cavity proper time is set to $T_Ω(ω_0) = 50$ ns. A schematic of this simple FSL configuration is shown in Fig. 2.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Illustration of an injected phase-modulated loop.</figcaption></figure>

## 6.1. Influence of modulation period tuning

In the presence of modulation, we assume that the repetition rate $Δτ$ matches the modulation period $T_M$. In this first example, we assume for now that we are working in a reference frame with $T_R = Δτ$. When the modulation period $T_M$ also matches the cavity proper time $T_Ω$, a simple expression for the pulse amplitude can be derived:

$$|u(\tau)| = \left| \frac{t_{21} e^{iS_M(\tau)}}{1 - t_{22} e^{iS_M(\tau)}} \right| \quad (70)$$

We run a first simulation with $T_M = T_Ω$ over 500 round trips, which is sufficient to reach the steady-state solution. In Fig. 3(a), we observe that the solution in the fast time window consists of two pulses, both coinciding with the zeros of the modulation function. The first pulse is associated with a positive slope of the modulation function, while the second pulse is associated with a negative slope. In Fig. 3(b), we verified that the central pulse indeed follows the theoretical analytical expression given in Eq. (70).

In a second simulation, we introduced a mismatch between the modulation time $T_M$ and the cavity proper time $T_Ω$, setting $T_M = T_Ω/1.005$. After 500 round trips, we observe in Fig. 4(a) that
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** (a) Real amplitude profile $\Re(u(\tau))$ of the field $u$ in the fast time window of width $\Delta\tau = T_R = 50$ ns, shown alongside the modulation function. Each pulse coincides with a zero crossing of the modulation function, one associated with a positive slope and the other with a negative slope. (b) Zoomed view of the central pulse, fitted with the theoretical expression of Eq. (70).</figcaption></figure>

both pulses are now chirped. Note that such chirps have already been observed experimentally in injected active phase-modulated loops [34,35]. The role of the injection in achieving a self-consistent solution after a round trip can be investigated further by comparing the pulses just before and just after injection. Fig. 4(b) shows that, apart from a small amplification of the pulse, which is expected when adding energy to the system, a time shift is also observed. A thorough numerical investigation shows that this time delay is precisely $T_M - T_\Omega$. The discrete injection induces a time shift to precisely compensate for the mismatch between $T_M$ and $T_\Omega$, allowing the solution to be self-consistent after a round trip of period $T_M$, even though the pulse is traveling along the ring at a speed corresponding to the group velocity $v_g = \ell/T_\Omega$. We refer the reader to Ref. [36] for a more detailed study of this phenomenon.

When we demonstrated the two-time-variable model and derived Eq. (16), we used a reference frame moving at the speed $v_R = \ell/T_R$, where $T_R$ is the round-trip time of the reference frame. In contrast, $\Delta\tau$, the periodicity of the fast-time window, can also be interpreted as the repetition rate of the system. It represents the time it takes for a pulse to become self-consistent again, during which it may change shape and even drift slightly as it traverses the various optical elements in the loop. Here, by choosing $T_R = \Delta\tau$, we eliminate the need to distinguish between the round-trip time of the pulse and the round-trip time of the moving reference frame. However, we find it useful to study below a case where this distinction matters, as it highlights the invariance of our formulation under a change of reference frame, in contrast to the derivations obtained from the Ikeda map.

## 6.2. Equivalence of reference frame change

For didactic purposes, since our example is simple enough, we will now explore what changes occur in the model if we choose a different reference frame and verify that both approaches lead to the same physical results. We will now consider a second reference frame moving at the group velocity $v_g = \ell/T_\Omega$. The first consequence is that the detuning term $\tau_d(\omega_0)$ present in Eq. (16) cancels out, and all occurrences of the time $T_R$ are replaced by $T_\Omega$ and $v_R$ by $v_g$. Concerning Eq. (22) modeling the modulation, we can no longer eliminate the reference frame round-trip index $k$ because the modulation function is not periodic in $T_\Omega$ when $T_\Omega \neq T_R$ (we still assume $T_R = T_M$). It is then necessary in the simulation to keep track of the round-trip index $k$ in order to
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** (a) Chirped pulses observed for a modulation period detuning $T_M = T_\Omega/1.005$. (b) Effect of the injection on the pulse intensity profile. Apart from a small amplification, both pulses also undergo a time shift.</figcaption></figure>

translate the initial modulation function by an amount $kT_\Omega$ after each round trip of the reference frame, in accordance with Eq. 23 as explained in Section 3.1. We anticipate that this discrete translation of the modulation signal will play an equivalent role to the term $\tau_d(\omega_0)\partial/\partial\tau$ present in the first approach.

We ran two sets of simulations with a detuning $T_M = T_\Omega/1.02$, one using our initial model with a reference frame moving at the speed $v_R = \ell/T_R$, and another using the alternative model with a reference frame moving at speed $v_g = \ell/T_\Omega$. The resulting amplitude profiles in the fast time window are represented in Figures 5(a) and (b), corresponding to 500 and 501 round-trips, respectively.

A notable fact is that for 500 round-trips of their respective reference frames, the stationary solutions superimpose perfectly. The reason is that the total elapsed time $\Delta T$ is a multiple of $T_R$, the true round-trip time of a pulse, in both cases. Indeed, in the first case, we have $\Delta T^{(a)} = 500 T_R$, while in the second case $\Delta T^{(b)} = 500 T_\Omega = 510 T_R$. However, when running the simulation for a time corresponding to 501 round-trips of their respective reference frames, both stationary solutions are shifted one with respect to the other. The reason is that this time, $\Delta T^{(b)} = 501 T_\Omega = 511.02 T_R$, which is not an integer multiple of $T_R$. This confirms that both approaches are physically consistent, as long as the boundary conditions imposed by the fast time window width $\Delta\tau = T_R = T_M$ are correct.

It is also worth noting that when the reference frame moves at speed $v_g$, the pulse drifts in the fast time window after each round-trip. However, this drift should not be interpreted as an inconsistency in the choice of $T_R$ in this case. This example was introduced for discussion purposes, and in practice, it is preferable to select a reference frame whose round-trip time $T_R$ matches the repetition rate $\Delta\tau$, ensuring that the pulse remains self-consistent and exhibits no drift after each round-trip.
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** The detuning is $T_M = T_\Omega/1.02$. (a) Stationary real amplitude profiles obtained after 500 round-trips of the respective reference frames, using both approaches. The dotted line corresponds to a reference frame moving at speed $v_R$, and the solid line corresponds to a reference frame moving at speed $v_g$. (b) Stationary real amplitude profiles obtained after 501 round-trips of the respective reference frames, moving at speeds $v_R$ and $v_g$, respectively.</figcaption></figure>

### 6.3. Consequence of incorrect boundary conditions

In the presence of modulation, the repetition rate $\Delta\tau$ of the system is forced to match the modulation period $T_M$. We can wonder what happens then if we set incorrect boundary conditions through a wrong fast time window width $\Delta\tau \neq T_M$. To illustrate this, we conducted another set of simulations with a detuning $T_M = T_\Omega/1.01$ and a fast time window width set to $\Delta\tau = 4 T_\Omega$. The results are shown in Fig. 6, comparing the two reference frame speeds considered earlier. We observe that incorrect boundary conditions result in non-physical artifacts in the computed solutions. The simulation forces the signal to fit within the user-defined window, even when it does not correspond to an integer multiple of the modulation period.

### 6.4. Breakdown of the mean-field injection model in the presence of modulation

In Section 3.3, while modeling an optical coupler, we stated that the mean-field formulation, which consists of adding a source term $u_{\text{in}}$ in the main propagation Eq. (16), is incorrect and requires the prefactor correction introduced in Eq. (34). Without modulation, this correction is negligible, as demonstrated by Eq. (S14) in Section S4.B of the Supplement 1, which explains the success of the LLE in predicting solutions for Kerr microresonators. However, we will see that in the presence of modulation, if one truly wishes to adhere to a mean-field formulation, this correction becomes essential. In order to emphasize the effect of the correction, we choose a larger detuning $T_M = T_\Omega/1.5$ compared to the previous examples.

The comparison between the predictions of the discrete injection model, the traditional mean-field model, and the corrected mean-field formulation is displayed in Fig. 7. We observe that the corrected mean-field formulation matches perfectly with the discrete injection model, whereas the traditional mean-field model can be completely off, exhibiting minima where maxima should occur. The correction introduced in Eq. (34) is derived from a simple linear model and
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Stationary solutions obtained using an incorrect fast time window width of $\Delta\tau = 4 T_\Omega$, for two reference frames moving at speeds $v_R$ and $v_g$, respectively. Only the most problematic period is shown.</figcaption></figure>

may also break down in the presence of strong nonlinear effects, either in the time or Fourier domain. For this reason, we recommend always avoiding mean-field formulations for accurate simulations. Instead, mean-field formulations should be reserved for analytical studies and a higher-level understanding of the system.

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Stationary solutions obtained using a detuning $T_M = T_\Omega/1.5$, for the different injection models: the discrete model, the mean-field model, and the corrected mean-field model.</figcaption></figure>

In conclusion, we have demonstrated that even the simplest passive resonator, which includes a modulator and an external drive field, must be modeled carefully to produce accurate simulation results. This example also illustrated the consistency of our derivation in Section 2, by highlighting the distinctions between the round-trip time of the reference frame, the repetition rate of a self-consistent pulse, and the cavity proper time at a given pulsation $\omega_0$. Additionally, it provided insight into how a pulse that travels at a speed corresponding to the group velocity within the loop can be temporally shifted by a discrete external injection in order to align with the modulation periodicity when the latter is detuned from the cavity proper time.

## 7. Conclusion

In this article, we rigorously derived the two-time-variable formulation which serves as a foundational framework for a range of dynamical models, including Haus master equations for
lasers and the Lugiato-Lefever equation for microresonators. Instead of accepting this formulation as an ansatz based on an Ikeda map, we started from the slow-varying envelope approximation and clarified the necessity of introducing time-periodic boundary conditions as the most effective approximation. We further expanded this formulation by incorporating auxiliary differential equations to model gain dynamics and by showing how to explicitly integrate discrete elements within the resonator, such as modulators, optical couplers, and saturable absorbers. In particular, our analysis demonstrates that explicitly treating an optical coupler as a discrete element, rather than using a mean-field model, improves accuracy without adding computational complexity.

We specifically highlighted that converting original spatial boundary conditions into time boundary conditions has significant implications for the solutions obtained and must be approached with care. This requires that the concept of round-trip time is clearly defined before tackling the problem. Moreover, in the transition regime, the solution should not exhibit substantial variations from one round-trip to the next, as initially proposed by Haus. This ensures that the assumption of time periodic boundary conditions does not give rise to unphysical solutions.

A primary limitation of the model presented here is its reliance on the traveling wave approximation, which considers only a single propagation direction and neglects spatial hole burning in the saturated gain medium. In future work, this framework could be extended to include bidirectional wave propagation by introducing corrections to the saturated intensity [37,38]. Such an extension would provide a more comprehensive description of systems where counter-propagating waves play a significant role.

Another promising direction for future research involves investigating systems with more complex dynamics, such as the coexistence of multiple pulse solutions traveling at different velocities or involving distinct polarizations [39]. Extending the model to account for these scenarios would enhance its applicability to real-world systems and provide deeper insights into their behavior. These perspectives offer a pathway to broaden the scope and accuracy of the foundational framework presented in this study.

Finally, through a didactic exploration of simple yet fundamental optical systems, including those detailed in the Supplement 1, we have gained a clearer understanding of the potential to generalize this method for modeling more complex resonators. This includes the ability to incorporate various loop sections with different materials and discrete optical elements placed at diverse locations. We anticipate that this method can be extended to arbitrary optical designs by employing a framework of cascaded differential equations with connected boundary conditions, like a block diagram, which can be solved iteratively.

**Funding.** Agence Nationale de la Recherche (ANR-21-CE42-0025-02 COCOA).

**Disclosures.** The authors declare no conflicts of interest.

**Data availability.** Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

**Supplemental document.** See Supplement 1 for supporting content.

**References**

1. A. Pasquazi, M. Peccianti, L. Razzari, et al., “Micro-combs: A novel generation of optical sources,” Phys. Rep. **729**, 1–81 (2018).

2. H. Guillet de Chatellus, E. Lacot, W. Glastre, et al., “Theory of talbot lasers,” Phys. Rev. A **88**(3), 033828 (2013).

3. H. Guillet de Chatellus, O. Jacquin, O. Hugon, et al., “Generation of ultrahigh and tunable repetition rates in CW injection-seeded frequency-shifted feedback lasers,” Opt. Express **21**(13), 15065–15074 (2013).

4. H. A. Haus, “Mode-locking of lasers,” IEEE J. Sel. Top. Quantum Electron. **6**(6), 1173–1185 (2000).

5. A. M. Perego, B. Garbin, F. Gustave, et al., “Coherent master equation for laser modelocking,” Nat. Commun. **11**(1), 311 (2020).

6. L. A. Lugiato and R. Lefever, “Spatial dissipative structures in passive optical systems,” Phys. Rev. Lett. **58**(21), 2209–2211 (1987).

7. L. Lugiato, F. Prati, M. Gorodetsky, et al., “From the Lugiato-Lefever equation to microresonator-based soliton Kerr frequency combs,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences **376**(2135), 20180113 (2018).
8. G. P. Agrawal, *Applications of Nonlinear Fiber Optics* (Academic press, 2020).

9. F. Teixeira, C. Sarris, Y. Zhang, et al., "Finite-difference time-domain methods," *Nature Reviews Methods Primers* **3**(1), 75 (2023).

10. A. F. Oskooi, D. Roundy, M. Ibanescu, et al., "Meep: A flexible free-software package for electromagnetic simulations by the FDTD method," *Comput. Phys. Commun.* **181**(3), 687-702 (2010).

11. R. J. LeVeque, *Numerical Methods for Conservation Laws*, vol. 214 (Springer, 1992).

12. R. J. LeVeque, *Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-state and Time-dependent Problems* (SIAM, 2007).

13. K. Ikeda, "Multiple-valued stationary state and its instability of the transmitted light by a ring cavity system," *Opt. Commun.* **30**(2), 257-261 (1979).

14. D. Kuizenga and A. Siegman, "FM and AM mode locking of the homogeneous laser - Part I: Theory," *IEEE J. Quantum Electron.* **6**(11), 694–708 (1970).

15. D. Kuizenga and A. Siegman, "FM and AM mode locking of the homogeneous laser - Part II: Experimental results in a Nd:YAG laser with internal FM modulation," *IEEE J. Quantum Electron.* **6**(11), 709-715 (1970).

16. H. Haus, "Theory of mode locking with a slow saturable absorber," *IEEE J. Quantum Electron.* **11**(9), 736-746 (1975).

17. H. A. Haus, "Theory of mode locking with a fast saturable absorber," *J. Appl. Phys.* **46**(7), 3049-3058 (1975).

18. U. Keller, *Ultrafast Lasers* (Springer, 2021).

19. F. X. Kärtner and U. Keller, "Stabilization of solitonlike pulses with a slow saturable absorber," *Opt. Lett.* **20**(1), 16-18 (1995).

20. A. Matsko, A. Savchenkov, and L. Maleki, "On excitation of breather solitons in an optical microresonator," *Opt. Lett.* **37**(23), 4856-4858 (2012).

21. P. Del'Haye, A. Schliesser, O. Arcizet, et al., "Optical frequency comb generation from a monolithic microresonator," *Nature* **450**(7173), 1214-1217 (2007).

22. T. Herr, V. Brasch, J. D. Jost, et al., "Temporal solitons in optical microresonators," *Nat. Photonics* **8**(2), 145-152 (2014).

23. T. Herr, M. L. Gorodetsky, and T. J. Kippenberg, "Dissipative Kerr solitons in optical microresonators," *Nonlinear Optical Cavity Dynamics: From Microresonators to Fiber Lasers* (Wiley, 2016), pp. 129-162.

24. T. Hansson and S. Wabnitz, "Dynamics of microresonator frequency comb generation: models and stability," *Nanophotonics* **5**(2), 231-243 (2016).

25. R. W. Boyd, *Nonlinear Optics, Third Edition* (Academic Press, Inc., USA, 2008), 3rd ed.

26. G. Moille, Q. Li, X. Lu, et al., "pyLLE: a fast and user friendly Lugiato-Lefever equation solver," *J. Res. Natl. Inst. Stand. Technol.* **124**, 124012 (2019).

27. H. Guo, M. Karpov, E. Lucas, et al., "Universal dynamics and deterministic switching of dissipative Kerr solitons in optical microresonators," *Nat. Phys.* **13**(1), 94–102 (2017).

28. I. I. Khanin, *Fundamentals of Laser Dynamics* (Cambridge Int Science Publishing, 2006).

29. E. Lacot, F. Stoeckel, and M. Chenevier, "Dynamics of an erbium-doped fiber laser," *Phys. Rev. A* **49**(5), 3997-4008 (1994).

30. U. Keller, "Ultrafast all-solid-state laser technology," *Appl. Phys. B* **58**(5), 347-363 (1994).

31. G. L. Lippi, J. Mørk, and G. P. Puccioni, "Numerical solutions to the laser rate equations with noise: technical issues, implementation and pitfalls," in *Nanophotonics VII*, vol. 10672 (SPIE, 2018), pp. 82-95.

32. H. Yang, M. Vallet, H. Zhang, et al., "Pulse doublets generated by a frequency-shifting loop containing an electro-optic amplitude modulator," *Opt. Express* **27**(13), 18766-18775 (2019).

33. H. Yang, M. Brunel, M. Vallet, et al., "Analysis of frequency-shifting loops in integer and fractional talbot conditions: electro-optic versus acousto-optic modulation," *JOSA B* **37**(11), 3162-3169 (2020).

34. H. Yang, L. Wang, C. Zhao, et al., "Sinusoidal frequency-modulated waveforms generated by a phase-modulated frequency-shifting loop," *J. Lightwave Technol.* **39**(10), 3112-3120 (2021).

35. M. Brunel, L. Frein, G. Loas, et al., "Nonlinear frequency chirps from a stabilized injected phase-modulated fiber laser loop," in *EPJ Web of Conferences*, vol. 287 (EDP Sciences, 2023), p. 07015.

36. H. G. de Chatellus, G. Loas, L. A. de Borggraeff, et al., "Slow and fast light effects induced by interference with a control light field," *Opt. Express* **32**, 35228-35235 (2024).

37. C. Sun, N. Mangan, M. Dong, et al., "Stable numerical schemes for nonlinear dispersive equations with counter-propagation and gain dynamics," *J. Opt. Soc. Am. B* **36**(12), 3263-3274 (2019).

38. N. Opačak and B. Schwarz, "Theory of frequency-modulated combs in lasers with spatial hole burning, dispersion, and Kerr nonlinearity," *Phys. Rev. Lett.* **123**(24), 243902 (2019).

39. M. Saha, S. Roy, and S. K. Varshney, "Polarization dynamics of a vector cavity soliton in a birefringent fiber resonator," *Phys. Rev. A* **101**(3), 033826 (2020).