Adaptive Time-Domain Simulation of Optical Cavities with Arbitrary Dynamics

1. Introduction

Time-domain simulation of optical cavities is required whenever the assumptions underlying frequency-domain or adiabatic models break down, such as in the presence of fast transients, short optical pulses, or rapidly varying boundary conditions. In these regimes, the intra-cavity field evolution depends explicitly on the finite round-trip time and the storage of optical energy [1].

This is particularly relevant in precision interferometry and gravitational-wave detectors, where the interaction between delayed intra-cavity fields and moving mirrors leads to non-linear, history-dependent dynamics that cannot be captured by stationary transfer functions [2,3].

Beyond interferometric gravitational-wave detectors, similar requirements arise in pulsed cavity ring-down spectroscopy, fast cavity tuning experiments, and active stabilization of high-finesse resonators, where transient dynamics and non-adiabatic effects determine system performance [4]. In all these contexts, such simulations provide direct access to transient regimes and provide a natural interface for modern data-driven control strategies.

Legacy interferometer simulation frameworks like e2e and SIESTA provide detailed optical modelling capabilities [5,6]. However, their use in real-time and reinforcement-learning applications is limited by their complexity and software design. The simulator presented here is designed to address transient regimes regimes with minimal approximations for a wide range of cavity parameters. It explicitly propagates the intra-cavity field as a sum over round trips, retaining the full temporal memory of the cavity [7,8]. Crucially, both the input optical field and the mirror positions can be modified at every simulation step. This allows arbitrary laser amplitude and phase modulation, as long as their characteristic frequencies remain below the simulator sampling frequency, as well as time-dependent cavity length variations under the same condition.

This level of flexibility is essential for our primary target application: integration into a Gymnasium-compatible reinforcement-learning environment for training real-time artificial intelligence (AI) agents dedicated to lock acquisition of suspended optical cavities currently under development [9,10]. Lock acquisition is inherently a non-stationary control problem, dominated by transients, delays, and non-linear field-mirror coupling [11]. Its stepwise formulation naturally interfaces with discrete-time control algorithms, enabling realistic training conditions and supporting an effective sim-to-real transfer [12].

In addition to the theoretical and numerical developments presented here, the simulator is released as open-source software [? ] under a separate license to ensure transparency and reproducibility.

The remainder of the paper is organized as follows. Sec. 2 introduces the numerical framework and the recursive formulation used to propagate the intracavity field, together with the implementation details of the simulation strategy. In particular we discuss the adaptive selection of the sampling frequency and the different integration regimes.

Sec. 3 presents the validation of the simulator through a comparison with ring down experimental data from the arm cavity of the Virgo interferometer cavity, while Sec. 4 extends the discussion to the broader applicability of the time-domain formulation to optics and photonics in general, as well as to machine learning (ML)-aided control strategies.

Finally, Sec. 5 summarizes the main results and outlines potential applications of the simulator, particularly in the context of real-time control and reinforcement-learning-based lock acquisition strategies.

Additional supporting material is provided in the Appendix, including the full comparison between simulated and experimental data from Virgo auxiliary channel over a 10 $\mathrm{s}$ time window, a reproduction of the pulsed cavity ring-down regimes reported by [4], and the complete benchmark table used to characterize the computational performance of the simulator.

2. Algorithm

The algorithm of the simulator is capable of reproducing non-linear dynamical effects in optical cavities caused by the ring-down phenomenon. The temporal evolution of the electric cavity field is described in [7] by the recursive equation:

$$\begin{array}{c}\hfill E\left(t\right)=t_a{\displaystyle \sum _{n=0}^{N-1}}\left[{\left(r_a{r}_b\right)}^n{e}^{-2ik{S}_n\left(t\right)}{E}_{\mathrm{in}}(t-2nT)\right]+{\left(r_a{r}_b\right)}^N{e}^{-2ik{S}_N\left(t\right)}E(t-2NT),\end{array}$$

(1)

where, $t_a$ is the amplitude transmissivity of the input mirror, $r_a$ and $r_b$ are the amplitude reflectivities of the input and end mirrors, respectively, $k=2\pi /\lambda$ is the optical wave-number associated with the laser wavelength, and $T=L/c$ is the one-way propagation time inside the cavity of length L. $E_{\mathrm{in}}\left(t\right)$ denotes the input optical field, while $E\left(t\right)$ is the intracavity electric field. The integer n labels the successive cavity round trips and N is the total number of round trips retained in the calculation. Eventually, the quantity $S_{\mathrm{n}}\left(t\right)$ represents the accumulated optical path over the previous n round trips and accounts for the mirror motion and cavity length variations,

$$S_n\left(t\right)={\displaystyle \sum _{p=0}^{n-1}}d(t-2pT),$$

(2)

specifically, $d\left(t\right)=L+(x_b(t-T)-x_a\left(t\right))$ where $x_b\left(t\right)$ and $x_a\left(t\right)$ are the end and input mirror positions. Equation (1) can be used to obtain the optical power $P\left(t\right)={\left|E\left(t\right)\right|}^2$ and the Pound-Drever-Hall (PDH) error signal [13], designed for optical cavity control at resonance. In the approximation where the modulation sidebands are fully reflected by the cavity and act as a stable reference field, the demodulated PDH signal can be expressed as:

$$V_{\mathrm{PDH}}\left(t\right)=-Im\left\{e^{ı\gamma }{E}_{\mathrm{in}}{\left(t\right)}^{*}E\left(t\right)\right\},$$

(3)

where $V_{\mathrm{PDH}}\left(t\right)$ is the PDH error signal, and $\gamma$ is the demodulation phase. $E_{\mathrm{in}}\left(t\right)$ is the complex conjugate of the input optical field, while $E\left(t\right)$ is the intracavity electric field.

These signals are the optical observables of the cavity state and become strongly non-linear under specific conditions, thereby making standard control technique for locking ineffective. In particular, ring down effect occurs when the mirror velocity becomes sufficiently large. Thresholds to this speed, referred to as critical velocities, can be defined [14,15]. In the present treatment, we consider, as a limit to define the non-linear regime, only the threshold obtained by equating the resonance crossing time to the cavity storage time, which represents the average time a photon remains stored in the cavity. This condition corresponds to the critical velocity defined as follows:

$$v_{\mathrm{cr}}\approx {\displaystyle \frac{\lambda \pi c}{4L{\mathcal{F}}^2}},$$

(4)

where, $\lambda$ is the laser wavelength, c is the speed of light, L is the cavity length, and $\mathcal{F}$ is the cavity finesse.

The simulator exploits the structure of Equation (1) to enable the selection of a sampling frequency over a wide range, as well as the independent motion of both the input and end mirrors. The input mirror displacement is modelled as a phase shift of the input laser field, as adopted in the algorithm.

Because the cavity field equation is recursive, each new sample depends only on an internal simulator state rather than on a full re-evaluation of the entire optical history. The cavity memory is therefore represented by a finite buffer of past input fields together with a single stored values of fields and mirrors position. This mechanism supports incremental, state-space-like updates instead of repeated global re-integration. Geometric electric field decay and phase factors are precomputed once during initialization and reused at every step, so the per-step computational load is dominated by multiply-accumulate operations.

The simulator allows the user to select the sampling frequency, but its final value is approximated according to an internal optimization rule designed to preserve consistency with the cavity round-trip time. We will refer to the sampling frequency desired by the user as $f_{\mathrm{calc}}^{\mathrm{desired}}$. The effective sampling frequency $f_{\mathrm{calc}}$ is computed through a rule-based decision tree that compares $f_{\mathrm{calc}}^{\mathrm{desired}}$ with $f_{2\mathrm{T}}=1/2T$, the round trip frequency. The main algorithm pseudocode is shown in Alg. 3. Alg. 1 and Alg. 2 clarify rounding strategies.

Algorithm 1 number_of_subhistories(${\eta }_{2T}$)

Require: ${\eta }_{2T}>0$ 1: if ${\eta }_{2T}<1$ then 2:     return $\mathrm{round}(1/{\eta }_{2T})$ 3: else 4:     return 1 5: end if

In particular, Alg. 2 implements a custom quantization rule that maps the continuous variable ${\eta }_{2T}$ to an integer order k used by the simulator. It computes $k_0=\lfloor {\eta }_{2T}\rfloor$, evaluates the switching threshold

$$b={\displaystyle \frac{2k_0(k_0+1)}{2k_0+1}},$$

and returns $k=k_0$ if ${\eta }_{2T}<b$, otherwise $k=k_0+1$. Alg. 1 and Alg. 2 are both used to round the dimensionless round-trip normalized time-step defined as ${\eta }_{2T}={\left(f_{\mathrm{calc}}^{\mathrm{desired}}2T\right)}^{-1}$.

Algorithm 2 round_for_inverse_curve(${\eta }_{2T}$)

Require: ${\eta }_{2T}>0$ 1: $k_0\leftarrow \lfloor {\eta }_{2T}\rfloor$ 2: $b\leftarrow {\displaystyle \frac{2k_0(k_0+1)}{2k_0+1}}$ 3: if ${\eta }_{2T}<b$ then 4:     $k\leftarrow k_0$ 5: else 6:     $k\leftarrow k_0+1$ 7: end if

Within the full algorithm, this function is the core discretization step. Its role is to choose the integer round-trip count N that best approximates the requested sampling condition, thereby controlling the trade-off between physical fidelity and numerical feasibility of the final computed sampling frequency $f_{\mathrm{calc}}=f_{2T}/N$. We will come back to this function further in the description.

Alg. 3 checks if $f_{\mathrm{calc}}^{\mathrm{desired}}>f_{2\mathrm{T}}$, meaning the requested time-step is shorter than a photon round-trip time. In this case, the algorithm interleaves multiple sub-histories, each shifted in phase to increase the effective sampling, and adjusts $f_{\mathrm{calc}}$ to an integer multiple of the round-trip frequency.

For $f_{\mathrm{calc}}^{\mathrm{desired}}<f_{2\mathrm{T}}$, we distinguish between two regimes. The first case occurs when the desired sampling frequency becomes too low compared to the cavity storage time; the corresponding integration window may exceed the effective memory of the cavity. In this case, the number of round-trips included in the summation is capped to a maximum value $N_{max}=5N_{\mathrm{eff}}$, where the effective photon lifetime expressed in round-trips is defined as

$$N_{\mathrm{eff}}={\displaystyle \frac{1}{\left|ln\left(r_a{r}_b\right)\right|}}.$$

The desired sampling frequency is then preserved, while the integration window is allowed to contain a fractional contribution of the last round-trip, ensuring that negligible long-past contributions are not propagated. This strategy guarantees temporal consistency while balancing numerical efficiency against the physical memory of the cavity.

The second case corresponds to intermediate frequencies, where an integer number of round-trips is selected using an inverse-curve rounding rule that finds optimal N, which minimizes the error with respect to the desired sampling frequency. In this case, the effective sampling frequency is set to $f_{\mathrm{calc}}=f_{2T}/N$.

To quantify the impact of the internal frequency-selection strategy, we evaluate the accuracy of the effective sampling frequency with respect to the user-requested value. The result is shown in Figure 1. The plot shows the relative accuracy

$$A_f=1-{\displaystyle \frac{\left|f_{\mathrm{calc}}-f_{\mathrm{calc}}^{\mathrm{desired}}\right|}{f_{\mathrm{calc}}^{\mathrm{desired}}}},$$

highlighting the minimal achievable deviation introduced by the approximation to values compatible with the cavity round-trip time.

The simulator workflow can be summarized in a few simple stages. First, the cavity is defined by specifying its physical parameters: mirror reflectivities, transmissivities, and cavity length, together with the desired sampling frequency $f_{\mathrm{calc}}^{\mathrm{desired}}$. During the initialization phase, the internal frequency-selection rule determines the effective sampling frequency $f_{\mathrm{calc}}$, the corresponding time step $\Theta$, and the number of round-trips included in each integration step. The effective sampling frequency is then returned to the user and must be used as the temporal grid of any external simulation environment. The simulation proceeds iteratively, updating the intracavity field at each step through the recursive field equation.

Algorithm 3 sampling_frequency_selection($f_{\mathrm{calc}}^{\mathrm{desired}},f_{2T},N_{\mathrm{eff}}$)

1: $N\leftarrow 1$   2: $n_{\mathrm{sub}}\leftarrow 1$   3: $\mathrm{partial}\_\mathrm{Theta}\leftarrow \mathrm{False}$   4: $\_N_{\mathrm{eff}\_\mathrm{factor}}\leftarrow 5$          ▹ field amplitude has decayed to $\approx 0.67\%$ of initial value   5: ${\eta }_{2T}\leftarrow f_{2T}/f_{\mathrm{calc}}^{\mathrm{desired}}$                  ▹ dimensionless round-trip normalized time-step   6: if $f_{\mathrm{calc}}^{\mathrm{desired}}>f_{2T}$ then   7:     $n_{\mathrm{sub}}\leftarrow \mathrm{number}\_\mathrm{of}\_\mathrm{subhistories}\left({\eta }_{2T}\right)$   8:     $f_{\mathrm{calc}}\leftarrow f_{2T}·n_{\mathrm{sub}}$   9: else 10:     $N_{max}\leftarrow \_N_{\mathrm{eff}\_\mathrm{factor}}·N_{\mathrm{eff}}$ 11:     if $f_{\mathrm{calc}}^{\mathrm{desired}}<{\displaystyle \frac{f_{2T}}{N_{max}}}$ then 12:         $N\leftarrow N_{max}$ 13:         $f_{\mathrm{calc}}\leftarrow f_{\mathrm{calc}}^{\mathrm{desired}}$ 14:         $\mathrm{partial}\_\mathrm{Theta}\leftarrow \mathrm{True}$ 15:     else 16:         $N\leftarrow \mathrm{round}\_\mathrm{for}\_\mathrm{inverse}\_\mathrm{curve}\left({\eta }_{2T}\right)\left[0\right]$ 17:         $f_{\mathrm{calc}}\leftarrow f_{2T}/N$ 18:     end if 19: end if 20: $\Theta \leftarrow 1/f_{\mathrm{calc}}$ 21: return$(f_{\mathrm{calc}},\Theta ,N,\mathrm{partial}\_\mathrm{Theta})$

3. Validation and Benchmarks

The simulator performance is evaluated through comparison with experimental data acquired from the north arm cavity of the Virgo interferometer [16,17] during mode matching measurement using the external excitation method, which consists of applying corrections to excite mechanically the end mirror in the longitudinal degree of freedom and analysing the second order mode during free swinging of the test mass. This mechanical excitation produces non-adiabatic crossings of the cavity resonances, exposing ringing effect for different crossing velocities. These represents optimal conditions for validating the simulation framework.

The selected dataset corresponds to a 10 $\mathrm{s}$ time window of the DC transmission power signal recorded at the cavity auxiliary B7 photodetector. Together with the transmission power, PDH error signal data from the B4 photodetector in reflection has been collected to compare with simulated data obtained from Equation (3).

The experimental data are shown in Figure 2 where it is also depicted the result of its processing. First, the first sample timestamp was set to zero. Then the transmission was detrended by subtracting its median value. Finally, resonance peak identification was performed with an algorithm based on a prominence-based peak detection: only peaks above ${99}^{\mathrm{th}}$ percentile of the detrended signal were retained. Each detected peak corresponds to a resonance condition $2L\left(t_n\right)=N_n\lambda$, with $\lambda =1064\mathrm{n}\mathrm{m}$. Successive resonances therefore correspond to cavity length variations of $\Delta L=\lambda /2$.

An integer resonance index n was assigned to each detected peak to reconstruct the cavity-length variation. The direction of the first crossing was determined from the sign of the local derivative of the PDH signal, after which the index was increased or decreased by one at each subsequent resonance. Turning points were identified based on a threshold for the transmission, which will be higher when the mirror velocity becomes small before changing direction. At these points, the sign of the index increment was reversed. The absolute cavity length for each n was reconstructed and expressed in terms of unity of Free Spectral Range (FSR). The discrete length samples were interpolated using a cubic spline to obtain a continuous estimate $L\left(t\right)$ over the oscillation window. A cubic-spline interpolation provides a smooth reconstruction of the resonance positions, but it is not equivalent to a sinusoidal model. Unlike a sinusoidal fit, the spline does not impose a global physical model on the mirror motion and may therefore introduce non-physical distortions in the derived velocity. For this reason, spline-based reconstruction is useful as a flexible interpolation tool, but can affect the velocity estimation.

The instantaneous mirror velocity was computed from the time derivative of the reconstructed cavity length $v\left(t\right)=\mathrm{d}L\left(t\right)/\mathrm{dt}$.

The computation of the simulated data was performed by initializing the simulation with cavity parameters consistent with arm cavity of the Virgo interferometer [18]: mirror reflectivities $R_a=0.986$, $R_b=0.99999$, input transmission $T_a=0.01377$, and nominal length $L=3000\mathrm{m}$. The input laser field was set to a mean value of $|E_{\mathrm{in}}|=12.7$, obtained by the least squares method fit described below. The calculation frequency was matched to the sampling frequency of the raw data which is 10 $\mathrm{k}$$\mathrm{Hz}$, ensuring temporal resolution well below the resonance crossing time. It is important to specify that in this case the actual calculation frequency approximated by the simulation strategy was $\approx 9.99\mathrm{k}\mathrm{Hz}$, with an approximation accuracy of $99.9$. With this values only one history of the optical field evolution was evaluated and the regime was the aggregated round-trips one, considering the full Equation (1).

The reconstructed velocity was used as input to the simulation, defining a time-dependent cavity length variation. The simulated cavity power was obtained from the intra-cavity electric field through the output mirror transmissivity.

For detailed analysis, one resonance has been chosen, corresponding to $\approx 4.9 \mathsf{\mu }\mathrm{m} \mathrm{s}$, i.e., a factor of $\approx 7.6$ critical velocities. The compared data are showed in Figure 3

In addition to the DC transmission, the simulated and reconstructed PDH error signals exhibit consistent zero-crossing slope and phase structure, confirming that the simulator correctly captures the dynamics induced by high-velocity mirror motion.

In the appendix we present an additional plot for validation of the simulator which consists of a comparison of the whole experimental data time window ($10 \mathrm{s}$) with some observations about integration errors.

In order to find the laser input power, Equation (1) can be approximated by a sum of adiabatic contribution plus a Doppler term:

$$E\left(t\right)\approx D_{0}exp\left(-{\displaystyle \frac{t}{\tau }}-ı{\displaystyle \frac{kv}{2T}}{t}^2\right)+{\displaystyle \frac{t_a{E}_{\mathrm{in}}}{1-r_a{r}_b{e}^{-2ıkvt}}},$$

(5)

where, $D_0$ is the complex amplitude of the transient contribution, $\tau$ is the cavity storage time, v is the mirror velocity, and $E_{\mathrm{in}}$ is the constant incident field. The prefactor $D_0$ is given by:

$$D_0=t_a{E}_{\mathrm{in}}\sqrt{{\displaystyle \frac{ı\pi }{2kvT}}}exp\left({\displaystyle \frac{ıT}{2kv{\tau }^2}}\right).$$

(6)

This approximation holds in the length-scan transient regime for $v\ge v_{cr}$. Its derivation assumes a continuum limit of the discrete round-trip dynamics $n\gg 1$ and timescales $t\gg T$. The immediate resonance region is neglected [7]; therefore, the least-squares fit is performed only on the tail of the ring-down.

In addition to the physical validation against experimental data, the computational performance of the simulator has been systematically assessed through a dedicated benchmark study. The benchmark on almost 400 cases compares the pure Python and Numba-accelerated backends on a test set of 99 cavity configurations spanning different lengths, reflectivities, and calculation frequencies. In all cases, the execution time scales approximately linearly with the number of round trips N, reflecting the dominant cost of the summation over stored field history. The pure Python implementation shows a clear increase in runtime with growing N, reaching several hundred microseconds for high-finesse, short cavities. In contrast, the Numba backend exhibits a much weaker dependence on N, reaching 16 $\mathsf{\mu }$$\mathrm{s}$ for the most demanding cavity tested, with $N=300$.

As a result, the speed-up factor increases systematically with N, ranging from 4.67 for $N\approx 1$ up to $17.61$ for $N=300$, demonstrating that just-in-time (JIT) compilation is necessary for efficient simulations in the high-finesse regime. Detailed results are shown in Table A1.

4. Relevance to Optics and Photonics

Beyond Fabry-Perot (FP) cavities, the time-domain simulator developed in this work may be of relevance to a broader class and physical systems characterized by energy storage, ringing effects, and finite propagation time, all of which exhibit memory-driven dynamics analogous to those of optical cavities. Within this framework, the same delay-based formulation naturally extends to integrated optical resonators, such as ring resonators [19], which can be treated as a unidirectional cavities where the intra-cavity field is recursively re-injected after each round-trip with attenuation and phase accumulation. Ring resonators constitute key building blocks of integrated photonic, including frequency-comb micro-resonators, Kerr non-linear devices, optical buffers, all-optical switching architectures, and quantum photonic circuits.

As an additional validation of the generality of our time-domain approach, we reproduce the characteristic ring-down dynamics obtained in [4], namely the transition between pulsed and smooth exponential decay regimes as a function of pulse duration and detuning. Figure A2 in the Appendix A presents simulated ring-down signals for different normalized pulse durations and detunings, showing qualitative agreement with the temporal feature described in [4], including the emergence of transient peaks and intensity modulation arising from the partial temporal overlap of successive photon round trips.

Another promising application of the simulator lies in the optimization of locking schemes for extremely high-finesse cavities ($\mathcal{F}\sim 10000$). In this regime, resolving the temporal evolution of the intra-cavity field is crucial to understand transient dynamics and control performance. Machine learning techniques, particularly reinforcement learning (RL), can exploit this information to identify optimal control strategies, avoiding direct interaction with experimental hardware, which can be both time-consuming and potentially harmful [9,10].

Finally, such simulations can support the design and optimization of gravitational-wave interferometers such as Virgo and Einstein Telescope, improving their design and further upgrades [20,21,22]. In particular, in a situation where phenomena like transient dynamics, mirror motion, radiation pressure effects,play a key role in maintaining stable operation.

5. Conclusions and Applications

In this work we presented the logic and the validation of a fast time domain optical simulator, capable of reproducing the dynamical evolution of the intra-cavity field in optical resonators. The simulator allows the user to flexibly define cavity parameters, boundary conditions, and mirror motion.

Computational efficiency is achieved through a simple yet effective simulation strategy, in which the user-defined sampling frequency is approximated to the nearest value compatible with internal optimization constraints. This adjustment ensures consistency between the discrete time step and the cavity round-trip time, preserving causal propagation while avoiding unnecessary evaluations of redundant round-trip histories. Efficiency is further improved thanks to the use of JIT compilation.

The resulting framework is designed to provide a computationally efficient and physically consistent simulation tool, spanning applications from large-scale interferometric systems to integrated photonic devices. Computational efficiency is particularly important for predictive modelling, but it becomes even more critical when considering real-time control development and design of RL-based lock acquisition protocols. Lock acquisition in suspended optical cavities is intrinsically a non-stationary control problem, in which the cavity state is observed through signals that deviate from the linear regime assumed by classical controllers.

The feature of a step-by-step evolution of the intracavity field and its associated observables (DC transmission, PDH error signal), under arbitrary time-dependent mirror motion and input field modulation, naturally matches the interaction paradigm required by modern RL algorithms, where an agent observes the system state at discrete time steps, applies an action (e.g., actuation on mirror position or laser frequency), and receives a reward based on the effectiveness of its actions in reaching the desired goal.

The development of framework empowered with this simulator enables large-scale training campaigns, in which thousands of interaction steps are required to learn robust policies. Moreover, addressing one of the biggest challenge in RL applications, the so called sim-to-real transfer, requires high-fidelity simulation of the physical system. Such fidelity is mandatory to train RL agents under realistic dynamics, thereby ensuring a safe and more reliable deployment on the real optical setups.