A stochastic nonlinear model of the dynamics of actively Q-switched lasers

In this paper, we present a novel stochastic and spatially lumped multi-mode model to describe the nonlinear dynamics of actively Q-switched lasers and random perturbations due to amplified spontaneous emission. This model will serve as a basis for the design of (nonlinear) control and estimation strategies and thus a high value is set on its computational efficiency. Therefore, a common traveling-wave model is chosen as a starting point and a number of model-order reduction steps are performed. As a result, a set of nonlinear ordinary differential equations for the dynamic behavior of the laser during a switching cycle is obtained. A semi-analytic solution of these differential equations yields expressions for the population inversion after a switching cycle and for the output energy, which are then used to formulate a nonlinear discrete-time model for the pulse-to-pulse dynamics. Simulation studies including models with different levels of complexity and first experimental results demonstrate the feasibility of the proposed approach.


Introduction
Nanosecond laser pulses are widely used in areas such as materials processing [24], spectroscopy [13], plasma research [28] and medical applications [3], where high pulse energies are required, and also in lidar systems [7]. The most established technology for achieving such pulses is active Q-switching. Since high pulse repetition rates are desirable for many applications [15], current research efforts are focussed on building actively Q-switched lasers which can operate at up to 1 MHz. At such frequencies and usual pump powers, large variations in output pulse energies which, in most cases, are detrimental to the intended application can ultimately damage optical elements inside the cavity. These energy variations are primarily the result of two distinct effects. Firstly, spontaneous emission is a stochastic process with known properties [16] which influences the shape of a pulse as well as its total energy ( [10], Figs. 12 and 13). Secondly, it is already known from regenerative amplifiers [5] that the pulse-to-pulse dynamics of pumping and inversion depletion exhibits instability regions where pulse energies are elements of nontrivial limit cycles: strong pulses which entail substantial population depletion and therefore strongly weaken the subsequent pulse's gain are followed by weak pulses which allow the populations to recover etc. In the case of actively Q-switched lasers, simulations ( [25], Fig.13) and experiments [1] have hinted at similar instabilities which ultimately yield (deterministically) chaotic behavior. Hardware-based techniques have been applied to mitigate stochastic [14] and dynamic effects [6] but limit the flexibility of the system. These limitations could be avoided by employing active feedback control strategies as in [5], which in turn requires accurate modeling of the system's behavior. A review of existing mathematical models of actively Q-switched lasers can be found in [26]. These models always involve some form of rate equation(s) for the lasing population densities and hence the optical gain of the active medium. A second set of equations then describes the evolution of optical energy within the cavity depending on the current gain and losses. On a more detailed level, we will differentiate mathematical models along three mutually independent key characteristics: 1. Spatially distributed or spatially lumped: In the former approach, population densities and optical powers are assumed spatially dependent (especially along the axis of beam propagation), while the latter disregards power transport phenomena and concentrates populations and powers into one point, see, e.g., [29]. Studies have shown [21,26] that spatially distributed modeling achieves superior accuracy especially for long active media.
2. Multi-mode or single-mode: Multi-mode models take into account the spectrally varying and potentially competing gain/loss characteristics of different longitudinal optical cavity modes as pointed out e.g. in [27]. Several techniques have been invented [31] to select a single longitudinal mode and avoid pulse energy variations resulting from mode-hopping in a mode competition. Single-mode models assume a single effective gain and loss and represent radiation by one longitudinal mode, preventing the simulation of spectral selection with the benefit of shorter overall simulation times.
3. Stochastic or deterministic: Contributions of spontaneous emission to the optical power are essential especially for the modeling of self-seeded lasers as they define the initial intracavity power level. If said contributions are represented by stochastic processes (deterministic quantities), we will call a model stochastic (deterministic). To the authors' knowledge, the stochastic nature of spontaneous emission has so far not been systematically considered in the literature on Q-switched lasers.
The current state of the art is the deterministic and spatially distributed multi-mode approach of travelingwave modeling as used recently, e.g., in [12,8,20]. Traveling-wave models feature spatially dependent rate equations as well as power transport equations for every longitudinal mode. The complex behavior resulting from the coupled nonlinear partial differential equations is useful for simulation studies, but renders an analysis of the pulse-to-pulse dynamics difficult. However, its superior accuracy to other existing models also makes traveling-wave modeling a good starting point for this paper. Using significant simplifications, Degnan [4] started off with a deterministic and spatially lumped single-mode model and arrived at a transcendental relation between the population inversions before and after a pulse. This is the first and, to the authors' knowledge, only explicit expression for pulse-to-pulse dynamics so far. Unfortunately, the assumption in this paper that the energy build-up always goes to the threshold (the point where the power gains equal the losses and no more amplification occurs) is no longer valid at higher pulse repetition rates, which makes this work not applicable to our problem. Therefore, we first present a first principles-based model of actively Q-switched lasers which employs a stochastic generalization of traveling-wave modeling where spontaneous emission is represented by a Bose-Einstein distributed noise term. For the inhomogeneous gain, this paper uses the exemplary case of Nd:YAG active media. Taking the general model as a reference, a number of successive simplification steps are then performed to derive a similarly accurate, but spatially lumped continuous-time stochastic multi-mode model and a discrete-time model which explicitly describes the pulse-to-pulse dynamics. Such a discrete-time model is beneficial for three reasons: Firstly, instead of computationally demanding and time-consuming numerical simulations for the complete Q-switching cycle, a single (semi-)analytic expression is given for the population inversion and output pulse energy, respectively. Secondly, the model and its derivation throw light on the question how existing models from the literature are connected, as we will discuss in Section 3.4. Finally, the main motivation for this model is the fact that it serves as a theoretical basis for the design of control and estimation algorithms to actively stabilize the pulse-to-pulse dynamics and to suppress stochastic influences on the intracavity power. Thus, Section 2 first establishes the general framework and the fundamental equations of the model. The step-by-step simplifications are then performed in Section 3. In Section 4, we show simulation results which illustrate the stochastic and dynamic contributions leading to energy fluctuations. Based on these simulation results, agreement with first measurements and theoretical results on pulse statistics and modeling errors introduced by the various simplification steps are assessed in Section 4 as well. The paper concludes with a summary of the presented work and an outlook in Section 5.  2 General model

Rate and transport equations
Contrary to regenerative amplifiers, the cavities of Q-switched lasers are constantly filled with continuouswave optical energy, while the pulse is created via the Q-switching process and is usually of greater duration than the cavity round-trip time. Throughout this work, we will make the usual assumption that no two laser modes of different wavelengths are mutually coherent and thus that no interference phenomena such as mode-locking occur. Therefore, radiation will be represented by powers instead of electric fields. Our starting point will be the (multi-mode) traveling-wave equations, where we use the example of a forward pumped Nd:YAG crystal within a linear (Fabry-Perot) cavity. A sketch of the system is drawn in Fig. 1. The spatial coordinate inside the cavity along which light propagates is defined as z-direction, and the active medium with population densities N i (z, t) resides within the interval [0, L]. Population densities and optical powers will always be assumed to be transversally uniform, which is why there will be no spatial dependences other than on z. At z = z P C , a Pockels cell with a reflection coefficient R(t) performs the Q-switching. All optical components apart from the active medium and the Q-switch determine the boundary conditions at either end of the cavity, which are specified in Section 2.2. For Nd:YAG active media, the most interesting lasing transitions occur at wavelengths λ A = 946 nm and λ B = 1064 nm. Their technological relevance is due to a small quantum defect and hence high theoretical efficiency in the former and a large emission cross-section in the latter [11], which renders the 1064 nm transition the dominant one. Pumping can be achieved at λ p = 809 nm from the ground state to an elevated state, which quickly decays to the upper state for both lasing transitions considered. Thus, our full model of the active medium effectively consists of five population densities as visualised in Fig. 2. The rate equations for the population densities N i (z, t) with Figure 2: The five energy occupation levels that constitute our model and their transitions: pumping (blue) and two lasing transitions (red, green). The black arrows stand for nonradiative thermal transitions.
the mentioned underlying energy level scheme and a total of M longitudinal cavity modes read [26] with relaxation rates γ ij for the populations and absorption and induced emission terms due to the pump power P p (z, t) as well as forward-and backward-running lasing powers P ± (z, t, λ j ). Wavelengths λ j belong to a discrete spectrum of cavity modes. The cross-sectional areas of the pump and the signal beams are given by A p and A s , respectively. The absorption/emission cross-sections (which we set equal for simplicity) are denoted by σ A (λ j ) and σ B (λ j ) for the transition wavelengths λ A and λ B , respectively. They determine the spectral lineshapes of the two transitions (which in general may overlap). N dop is the crystal's total Neodymium dopand density, and N 1 (z, t) can be calculated from the other population densities via the algebraic equation . For setups such as claddingpumped lasers, one also needs to include appropriate overlapping factors Γ, see, e.g., [26]. Because active media usually have very similar refractive indices at comparable wavelengths λ A , λ B , and λ p , the group velocities of radiation in this spectral range are all approximately equal to a common value v.
The power transport equations can then be written as ( [23], Cha. 8) where the first terms on the right-hand sides stand for absorption and induced emission due to interaction with the population densities N i (z, t). The second terms represent effective losses, for example diffraction, with the coefficients α p and α(λ j ), respectively, and α RS is the Rayleigh backscattering coefficient (capture factor included). The final term on the right-hand side of (25b) accounts for spontaneous emission, where 2hc 2 /λ 3 j ∆λ j = 2hν j ∆ν is the contributed power of spontaneously emitted photons of 2 independent polarisations belonging to a mode λ j with a spectral spacing ∆λ j between modes into a solid angle ∆Ω where photons are kept inside the cavity [30]. Moreover, ζ(z, t, λ j ) is random in nature and describes the number of photons of wavelength λ j spontaneously emitted in the region [z, z + dz) at time t per length increment dz. The term will be further discussed in Section 2.3.

Cavity boundary conditions
As can be seen in Fig. 1, optical elements (whose effective round-trip cavity efficiency is η(λ j )) bring along that the pump and beam powers are subject to boundary conditions (B.C.) For R(t) = 1, the Pockels cell introduces a discontinuity at z P C given by with the resulting output power From this and the Q-switching frequency f switch , we get the m-th output pulse energy as

Stochastic seeding via spontaneous emission
First, we define the total scattering cross-section We will use it for the ensemble mean of ζ(z, t, λ j ), which corresponds to the inverse mean free path. It can be shown quantum-mechanically for a highly simplified model [17], where for instance, all longitudinal modes experience identical gains, that amplified spontaneous emission in K independent modes belonging to one bosonic excitation state follows the K-fold degenerate Bose-Einstein distribution, which for large K approaches a Poisson distribution and for low K is significantly broader [18]. Since we simulate every cavity mode individually, we can set K = 1 and therefore propose that the probability of n photons of wavelength λ j being spontaneously emitted within [0, z] at time t follows a single-mode Bose-Einstein distribution with mean Hence, the number of photon counts, is an inhomogeneous jump process following the probability distribution (8).
The dominant mode λ j of a pulse arises from a competition among modes around the maximum of σ A (λ j ) or σ B (λ j ) due to the greatest amount of spontaneously emitted photons in the early phase of the energy build-up. Bragg gratings and other wavelength-selective elements inside the cavity can shape η(λ j ) in the boundary conditions (3a-3c) such that the setup exhibits a few or many significant modes. As a consequence, the experimental output pulses will be approximately Bose-Einstein or Poisson distributed, as was shown in [18].

Transport equation solution & total rate equations
For control-related tasks, numerically solving the coupled rate (1a-1d) and power transport (25a,25b) equations with the boundary conditions (3a-4a) is not feasible as systematically investigating the pulse-to-pulse dynamics based solely on simulations is difficult. To eliminate spatial dependences, we introduce the total populations and make a number of simplifications which allow to express the dynamical behavior of the laser in terms of these new variables. We start the model simplification by choosing frames t → t ∓ z/v in (25a,25b) which move with a point of the pump or signal beam. Additionally, we make the following two assumptions: (A1.) When the forward-running (the backward-running) laser beam travels through the active medium starting at t k , the population densities and the lasing powers are assumed to be temporally constant, i.e. (2). Moreover, for the Rayleigh backscattering contribution α RS P ∓ (z, t, λ j ), the dependence of P ∓ on z will be neglected, i.e. P − (z, t, λ j ) = P − (L, t, λ j ) and P + (z, t, λ j ) = P + (0, t, λ j ). The transport phenomena can therefore be described completely in terms of the forward-and backward-moving frames t → t ∓ z/v. (A2.) Photons that are spontaneously emitted or backscattered at z ∈ [0, L] join the forward-running (the backward-running) beam at z = L (z = 0). In the supplemental document, we provide the formulae representing this assumption. With initial conditions given by the incident powers P p (0, t k ) = P p (t k ) and P ± ( L∓L 2 , t k , λ j ) for the pump and signal beams before transmission, respectively, and separating each two channels' frequencies by one Figure 3: The overlap between lineshapes of scattering cross-sections σ A (λ) and σ B (λ) is eliminated by lineshape separators g A (λ) and g B (λ).
inverse cavity round-trip time, i.e. ∆ν = 1/t RT , we show in the supplemental document that the power transport equations (2) can be solved analytically provided that the assumptions (A1.) and (A2.) hold. The solutions, i.e. the outgoing powers corresponding to any incoming ones, read where the spontaneously emitted photon count C(t k , λ j ) := C(L, t k , λ j ) is a random number with distribution (8) and mean µ(L, t k , λ j ) according to (9), i.e. the value the counting process has reached between z = 0 and z = L. As a next step, we want to find rate equations for the total populations defined by (11). To this end, we make use of the following assumptions: (A3.) When calculating the depletion of population inversion due to lasing powers or its increase due to the pump power, dependence on z of the loss terms will be neglected, i.e. α(λ j )P + (z, t, λ j ) = α(λ j )P + (0, t, λ j ) and α(λ j )P − (z, t, λ j ) = α(λ j )P − (L, t, λ j ) as well as α p P p (z, t) = α p P p (0, t).
(A4.) For the induced emission contributions to the lower populations N 3 and N 2 in (1a-1d), spectral overlap of the lineshapes σ A (λ) and σ B (λ) can be disregarded. Assumption (A4.) may be implemented using lineshape separator functions g A (λ) and g B (λ) which are either 1 or 0 and switch support at some separation wavelength λ AB . The effect of lineshape separators can be seen in Fig. 3. With assumptions (A1.)-(A4.), we show in the supplemental document that the total rate equations read as

State reduction via slow-fast dynamics
The model so far considers 5 populations of which only 4 are independent due to the algebraic condition involving N dop . If the relaxation rate γ ik of an excitation state population N i is very large, the state is depopulated nearly instantly to N k . In that case, which is usually a good approximation for low energy gaps ∆E ik in the order of phonon energies, we can apply singular perturbation theory to reduce the order of the model. Assuming that the active medium is at thermal equilibrium (thus introducing a thermalization factor [22] and employing a linear state transform, we show in the supplemental document that the total rate equations (13a-13d) consist of a slow and a fast dynamic subsystem. The limit of infinitely fast relaxation yields a quasi-stationary model for the parameters which results in a slow dynamics that can be expressed solely in terms of the total population In light of (13e), a reasonable choice of the spatially lumped intracavity power can be given by the power directly applying to the active medium, which is given by In the limit (14) with (15), (16) and (6) and defining b := 1/(1 + B 45 ) as well as γ := γ 42 + γ 43 , we show in the supplemental document that one can obtain from (13a-13d) the reduced total rate equation

Continuous-time nonlinear model (CTNM)
Since we are mainly concerned with the evolution of the intracavity power rather than the transport effects inside the cavity, we will average the transmission through the active medium, the cavity losses and the Q-switch influence over half a cavity round-trip for a linear cavity, as was done in [4]. For this, the power transport solution (12b) and the boundary conditions (3b,3c,4a) are utilized and for simplicity, R(t) is applied to both forward-and backward-running beams instead of R(t) to the backward beam. We also assume that amplification and losses also apply to spontaneously emitted or backscattered photons. Then, half a cavity round-trip, i.e. a time increment of ∆t := t RT /2 = L cav /v, affects the powers according to The time derivative of (16) can be approximated by which with (16) and (18) and assuming that spontaneous emission and Rayleigh backscattering are small effects (i.e. ln(1 + x) ≈ x) yields an ordinary differential equation which is derived in the supplemental document and reads as where we have defined the loss rate Again, C(t, λ j ) is a Bose-Einstein noise term with expectation value E(C(t, λ j )) = σ(λ j )N (t) that can be understood in a Langevin sense. More formally, one could write (49b) as a stochastic differential equation involving a renewal-reward process. We show in the supplemental document how the special case of a compound Poisson process with exponentially distributed waiting times and Poisson distributed jumps can be used to model the emission of photons at random times while producing equivalent means and variances after a time increment of t RT /2 as the Bose-Einstein distribution from (8). Neglecting time delays and Rayleigh backscattering during half a round-trip, we show in the supplemental document that the output power, originally given by (4b), can be obtained from the spatially lumped intracavity power by while the output energy is still given by (5). Equ. (49b) and the reduced total rate equation (49a) provide a very compact description of the system dynamics during a switching cycle and constitute what we will call the continuous-time nonlinear model (CTNM) of actively Q-switched lasers. In order to obtain the population and power after a switching cycle (and hence the pulse-to-pulse behavior of the system) from the CTNM, one would have to solve the coupled, nonlinear differential equations (49a,49b), which is not possible without further approximations, as we will discuss in Section 3.5.

Comparison to existing models
It is an interesting observation that our modeling effort so far bridges the gap between spatially distributed traveling-wave modeling (where we started our investigation) and the more traditional point models which are spatially lumped. From the spatially lumped total rate equation (49a), the rate equation used in traditional point models can be derived by isolating the pump and relaxation terms, additionally assuming losslessness and non-pumping-saturation, i.e. α p = 0 and LN dop N (t), whose solution yields equation (5) in [2], and assuming low single-transmission gains, i.e. exp(σ . This indicates that the extra terms in (49a) are effects of the spatially distributed nature of the amplification process of the intracavity power as well as of the absorption of the pump power neglected by traditional point models. These modifications have a significant effect on the model's behavior during a switching cycle and affect the pulse-to-pulse dynamics, especially for low to medium switching frequencies. Furthermore, the power equation used along with a point model rate equation in Degnan's model [4], which has been used extensively for experimental laser optimisation [2], can be retrieved from (49b) by neglecting the last terms accounting for spontaneous emission and backscattering and disregarding multi-mode behavior. To the best of the authors' knowledge, the model proposed in this paper is the first one on Q-switched lasers where the stochastic nature of spontaneous emission is incorporated.

Pulse-to-pulse dynamics
As pointed out in Section 1, the output pulse energies of an actively Q-switched laser are subject to fluctuations due to stochastic and dynamic effects. The former mainly come from amplified spontaneous emission, while the latter are a result of unstable pulse-to-pulse dynamics due to a coupling of subsequent pulses. Since Q-switching (e.g. a Pockels cell represented by R(t)) operates in a cyclic manner, the pulse-to-pulse dynamics can be described by a discrete-time dynamic system where each switching cycle is a discrete event. If t m is the time immediately after the (m − 1)-th pulse has been coupled out, one switching cycle evolves according to t m+1 = t m + 1/f switch . We use the m-th pulse energy E m given by (5) and (21) and the abbreviation N (t m ) := N m . While most generally, N m and P j (t m ) would both be dynamic variables, we can argue that P j (t m ) ≈ 0 for all m and j since all intracavity power has just been coupled out at those times, which reduces the dimension of the dynamical system from 1 + M to 1. With future control applications in mind, typical feedback quantities include the measured pulse energy, i.e. E m is the considered output quantity. As input quantity, we can choose the intracavity power P j (t m + t pump ) =: u j,m after the pumping phase and at the beginning of the m-th energy build-up. This input quantity is either given by spontaneous emission (uncontrolled case) or by a defined level from external seeding or controlled prelasing [14] (controlled case). Thus, defining the pulse-to-pulse and output maps by f and h, respectively, the discrete dynamics are given by These mappings are either found from numerical solutions of the general model discussed in Section 2, numerical solutions of (49a,49b) with (21,5) or an approximate analytic solution. Thus, in order to arrive at explicit expressions for (22), we will separate the switching cycle into its essential sub-processes and approximately solve equations (49a) and (49b) for each sub-process, incorporating additional fitting parameters. The three essential sub-processes of a switching cycle are: 1. pumping and population relaxation for a duration of t pump , 2. energy build-up during the duty cycle t bu = DC/f switch (DC is the relative duty cycle), and 3. coupling-out of the pulse with pulsewidth t coup ≈ t close , where t close is the closing time of the cavity.
The composition of these solutions then constitutes our discrete-time dynamic model of the pulse-to-pulse dynamics. This semi-analytic procedure can be done in different ways, of which the authors explored a few.
We present a detailed derivation in the supplemental document whose results are relatively simple and can be fitted effectively to adjust the dynamics to the observed behavior. What renders the procedure especially difficult is a combination of exponential growth underlying the energy build-up and nonlinearities in both (49a) and (49b). The nonlinearities stemming from spatially distributed amplification and absorption are essential to accurately describe the dynamics, especially for low to medium switching frequencies. Therefore, linearisation or omission of loss terms immediately leads to substantial modeling errors. In the supplemental document, we show that one can find explicit expressions for the pulse-to-pulse dynamics of the form where f 1 represents pumping, f 2 the energy build-up and f 3,N and f 3,E the effect of coupling-out on inversion and pulse energy, respectively. Such explicit expressions provide an efficient tool for analyzing the resulting laser dynamics as well as designing real-time capable control and estimation algorithms.

Simulation results
In the following, we will show simulation results which demonstrate the interplay of stochastic and dynamic effects, agreement with first measurement results and the accuracy of the model simplifications proposed in the previous section.

Pulse statistics and bifurcations
Apart from noise on the intracavity and output powers during a switching cycle, spontaneous emission also causes variations in the output energies even when the system is dynamically stable. To study the interplay of stochastic and dynamic effects in our general model discussed in Section 2 and given by (1)(2)(3)(4)(5), the duty cycle (DC) of the Pockels cell is varied at a fixed switching frequency (f switch = 1 MHz) and pump power (P p = 22 W). By plotting the output energy histograms and deterministic output pulses for comparison, a bifurcation diagram (Fig. 4) is obtained which displays how period-doubling bifurcations and increasing losses due to a longer build-up affect the energy distribution. Since a cavity design with large bandwidth was assumed, we have a large number of contributing modes and hence expect roughly Poissonian output energy statistics, as was pointed out in [18], which is indeed the case, cf. Fig. 5b). On top of the stochastic variations, period-doubling bifurcations at relative duty cycles of around 34% and 42% can be observed. From the third period-doubling onwards, we get deterministic chaos with an intermittent window of a 3limit cycle between 53 and 56 %. The mean output energy diminishes at such build-up times since the amplification phase has no significant gain left and losses dominate instead.

Validity of the CTNM
The continuous-time nonlinear model (CTNM) given by (49a,49b) together with (21,5) simplifies the spatially distributed full model (1)(2)(3)(4)(5) by averaging all processes within the cavity over half a round-trip such that ordinary differential equations remain to describe the whole system dynamics. We will show in Section 4.3 (cf. Fig. 6) that this introduces only small model errors as the full model and the CTNM agree well both in stable and unstable regimes of operation. Fig. 5 demonstrates that the CTNM can reproduce experimental data in terms of dependence on P p and the DC as well as theoretical predictions regarding the photon statistics of output pulses. The measured average pulse energies shown in the top row stem from a double-pass diode-pumped Yb: obtained from stochastic pulse energies E m given by (5), similarly to how it was done in [18] with experimental data. The histograms come from simulations of an Nd:YAG system at 1 MHz switching frequency and 20% DC exhibiting stable pulse-to-pulse dynamics, where the cavity supports around 270 contributing longitudinal modes (bottom left) and a single mode (bottom right). For comparison, the red lines in the bottom plots of Fig. 5 represent a Poisson and a Bose-Einstein distribution and hence the theoretical predictions from [17], respectively. It should also be noted here that the CTNM and the full model exhibit equivalent pulse energy statistics.

Comparison of models
Since we have shown in Section 4.2 that the CTNM (49a,49b,21,5) is able to reproduce the behavior of experimental results and theoretically predicted photon statistics of output pulses, we will now assess whether the CTNM differs significantly from the full model (1)(2)(3)(4)(5) and how further simplifications compromise the model's accuracy by analysing a collection of simulation results shown in Fig. 6. One simplification lies in reducing the system to single-mode operation, neglecting all other cavity modes in (49a,49b) besides the dominating one (i.e. the mode subject to the highest gain). Another simplification of the CTNM is its approximate analytic solution, i.e. the discrete-time model (23) of the pulse-to-pulse dynamics. Fig.  6 contains two exemplary cases for a system at f switch = 1MHz repetition rate and a pump power P p = 22.5 W exhibiting stable (a: 8% DC) and unstable (b: 35% DC) pulse-to-pulse dynamics. They are simulated using the full model (1)(2)(3)(4)(5), the CTNM (49a,49b,21,5), a single-mode version of the CTNM and the discrete-time pulse-to-pulse model (23) discussed in Section 3.5. Here, N s , u s and E s stand for the steady- state population, uncontrolled intracavity power and pulse energy of the pulse-to-pulse dynamics (22). Row 3) of Fig. 6 demonstrates that increasing the DC value alters the shape of the pulse-to-pulse dynamics curve f (N m , u s )/N s . Increasing one of the other two bifurcation parameters f switch or P p moves the steady state (position of the grey arrow) of the pulse-to-pulse dynamics further to the right of the dynamical map f (N m , u s ) (to higher N m ) such that eventually the absolute value of the slope becomes larger than one and the dynamics hence unstable. In Fig. 6, we can see that the single-mode approximation (yellow) introduces major modeling errors as all photons are emitted into the channel with maximum gain. This in turn leads to a premature energy build-up with enormous depletion and power losses. One can further see -most distinctly in the insets of 1a) and 1b) -that the population inversions of the full model and the CTNM differ only by about 1 percent. The discrete-time model matches these dynamics relatively well, but in the instability region of the pulse-to-pulse dynamics, the limit cycles for populations in 1b) as well as pulse energies in 2b) have smaller amplitudes. In terms of pulse energies E m , row 2) again displays a good agreement between the CTNM and the full model. The normalised pulse-to-pulse dynamic maps in rows 3) and 4) of Fig. 6 for the full model, the CTNM and the discrete-time model stay in close proximity. However, it should be emphasized that the discrete-time model fits the dynamics well up to around 1.1 times the steady-state population N s in rows 3) and 4). At larger inversion, too much population is depleted in the discrete-time model, and the output map h(N m , u s ) overshoots significantly. This is due to the fact that losses can only be incorporated via operator splitting, as we showed in the supplemental document. In practice, the overshoot is not problematic at all as such high inversions are never reached if the cavity is being opened and closed periodically, but an analysis on a broad domain of inversion values is still useful for assessing the global dynamical behavior of the models.

Conclusion and Outlook
In this work, a computationally efficient mathematical model for an actively Q-switched laser system including the effects of stochastic spontaneous emission is derived as a basis for the design of advanced control and estimation strategies. A spatially distributed multi-mode traveling-wave model, which is commonly used in the literature, serves as a starting point. Based on this model, the power transport equations were simplified by making physically feasible assumptions and the singular perturbation theory was applied to systematically reduce the number of relevant populations. Then all the processes within the cavity were averaged over half a cavity round-trip. As a result, a set of coupled, nonlinear ordinary differential equations were obtained, which describe the evolution of the inversion and each cavity laser mode at any point in time. Moreover, in order to separately study dynamic and stochastic effects on energy fluctuations, a discrete-time model of the pulse-to-pulse dynamics was derived by solving the differential equations in a semi-analytic way. Simulation results show that the model simplification steps only introduce small errors in the range of sub-percent, except for the discrete-time model at large populations, which does not constitute an issue for the practical application. In particular, the systematic incorporation of multi-mode behavior and spatially distributed amplification, which distinguishes the spatially lumped continuous-time nonlinear model presented in this work from traditional point models known from the literature, provides superior modeling accuracy with respect to the full model at moderate computational costs. Thus, the continuous-time nonlinear model accurately describes the behavior of actively Q-switched lasers during a switching cycle, while the presented discrete-time model provides a simple relation between the initial populations and the corresponding output energies of two different cycles. The simulation results show that stochastic seeding leads to amplified spontaneous emission noise whose pulse statistics correspond to theoretical results from the literature, while the simulated slope efficiencies of output energies for various build-up times agree with first measurement results presented in the paper. Moreover, as it was expected from experimental observations and other papers from the literature, under certain conditions the nonlinear pulse-to-pulse dynamics exhibits period-doubling bifurcations and deterministic chaos. The model presented in this paper aims at laying the foundation to systematically design control and estimation strategies to achieve robustly stable Q-switched lasers at high repetition rates. In particular, the pulse-to-pulse dynamics can be stabilised or otherwise modified via active feedback methods and prelasing strategies can be employed to suppress stochastic seeding due to sponta- &710VLQJOHPRGH GLVFUHWHWLPHPRGHO

A stochastic nonlinear model of the dynamics of actively Q-switched lasers: supplemental document
A Approximate solution of the power transport equations As mentioned in the paper, we will present here how the assumptions (A1.) and (A2.) lead to an approximate solution of the power transport equations (2) in the paper. In this section, we will omit wavelengths λ j in all function arguments for conciseness. Starting with (2), we introduce the moving frame by t → t ∓ := t ∓ z/v. As a consequence, the partial derivatives in z transform as by the chain rule (for the pump power analogously). The power transport equations then take the form We will now focus on forward-running signal beams (pump and backward-running beams work analogously).
By assumption (A1.), (25b) can be solved analytically, and we get We see that the exponential within the integral cancels a part of the outer exponential. This means that spontaneously emitted or backscattered photons experience a gain which depends on where they are emitted. By virtue of assumption (A2.), photons joining the P + beam do so at the end of the active medium (z = L), and hence, Thus, the exponentials exactly cancel, and writing out time arguments again, we can directly evaluate all integrals since we have defined the total populations N tot i (t) as integrals over N i (z, t) and the photon count C(t) as integral over ζ(z, t). As written in the paper, the wavelength spacing is chosen correspondingly to each two frequencies being separated by an inverse cavity round-trip time. Hence, we obtain which, if we insert it into (26) together with (27), yields the solutions given by (12a,12b) in the paper.

B Derivation of the total rate equations
We will now derive the total rate equations given by (13a-13e) in the paper by using assumptions (A1.)-(A4.). For demonstrational purposes, we will show the derivation of , but the concept applies to all equations equivalently and the computations are analogous. First, we write the rate equation (1c) as For the second row of (29), we implement assumption (A4.) by employing the lineshape separator function g B (λ). For the rate equation (1d), this would be g A (λ), and for (1b), it would be g A (λ) + g B (λ) ≡ 1. This means that in (29), we have If we now look at the power transport equation in a moving frame, i.e. (25b), with (28) already applied, we observe that Inserting (30) into (29), which yields and integrating both sides over z, the integrals of the time derivative and relaxation terms can be performed immediately due to the definition of N tot i (t) in (11) of the paper. For the depletion term Γ(z, t, λ j ), the integral of (31) reads where we used the definition of the stochastic term provided in (10) of the paper. For P + (L, t, λ j ) and P − (0, t, λ j ) in (32), we insert the power transport solution given by (12b) in the paper and derived in Section A. The spontaneous emission terms then cancel, as do the α RS terms due to assumption (A1.). The integral of the term in (32) which involves α(λ j ) can be directly computed because of assumption (A3.), and we obtain L 0 α(λ j ) P + (z, t, λ j ) + P − (z, t, λ j ) dz = α(λ j )L P + (0, t, λ j ) + P − (L, t, λ j ) , Finally, the resulting differential equation yields (13c) if we make the substitution

C State reduction via slow-fast dynamics
In the paper, we stated that the reduced total rate equation (17) can be derived as a quasi-stationary model from (13) by applying singular perturbation theory. More concretely, the limit of infinitely fast relaxation within the upper and lower state manifolds which is expressed by (14) leads to a slow and a fast subsystem of (13). We will now pursue this path in detail, derive (17) and assess the approximation order of the reduced model via Tikhonov's theorem (Theorem 11.4 in [9]).
In this section, we introduce the Boltzmann distribution thermalization factors B ij := exp −∆Eij k B T and use the abbreviations The total rate equations (13) then read as As a next step, we employ a linear state transformation where we have y 4 + y 1 = LN dop . Inversion of the transformation yields By differentiation of (34) and insertion of (33), one obtains After inserting (35) into (36), we set up the slow-fast dynamics by assuming thermal equilibrium, i.e. If we now define the singular perturbation parameter := 1/γ 54 , we get where we have not transformed the arguments of H and G j for brevity. Multiplying (37a,37c,37d) by and defining the slow state variable x := y 4 as well as the fast state variables the system (37) takes the singular perturbation standard form with the dynamical maps −(B 12 a 31 y 3 + (1 + B 12 )a 21 y 2 ) The quasi-stationary model is derived from (38) via the limit → 0 and reads as dx r dt = f 1 (t, x r , z r , 0) , with z r = g(t, x r ) as the solution of the algebraic equations f 2 (t, x r , z r , 0) = 0. This yields z r = g(t, x r ) = 0 and hence the quasi-stationary model reads as According to Tikhonov's theorem (see Theorem 11.4 in [9]), the so-called boundary layer model in the fast time scale τ must be uniformly exponentially stable in t and x s , which is the case because Applying Tikhonov's theorem further proves that there exist finite values of up to which the reduced system (40) is of approximation order O( ) with respect to (38). For more details on the singular perturbation theory, the reader is referred to, e.g., [9]. Finally, the slow dynamics (40) is transformed back to N i via (35) and by employing y 5 = y 3 = y 2 = 0. This yields which, keeping in mind that y 4 = x and defining N := N 4 , endows us with the reduced total rate equation One final approximation lies in setting B 12 = B 13 = 0 (large energy gaps ∆E 12 and ∆E 13 such that upwards relaxation is impossible), which entails B = 1 and after which we arrive at (17) in the paper.

D Continuous-time nonlinear model
We will now show how (20) can be obtained from (18) and (19) in the paper. First, we write (18) for forward-and backward-running powers separately, i.e.
Taking the sum of (42a) and (42b) and using the definition (16) of P j (t), we obtain and thus, by virtue of (43). The argument of the logarithm in (44) consits of three factors which can be separated via ln(abc) = ln(a) + ln(b) + ln(c). Multiplication of (44) by Pj t RT /2 yields the first two terms of (20). The third and fourth terms come from the approximation that spontaneous emission and backscattering are small effects, i.e.
by a Taylor approximation, which completes the derivation of (20).
In (21) of the paper, we give an expression for the output power as a function of the intracavity power under the approximation of disregarding time delays inside the cavity as well as Rayleigh backscattering. These approximations must be made to avoid an infinite regression in time and obtain a closed-form expression. The given formula can be derived by firstly stating that forward-and backward-running powers must equal the intracavity power and secondly applying the boundary conditions (3,4) as well as the power transport solutions (12) given in the paper. In total, this means that we first solve the system of equations that is Note that by (4a), R(t) appears here instead of its square root, which was applied to both forwardand backward-running beams, e.g. in (42a) and (42b), as an approximation to derive (43). Solving for P − (L, t, λ j ), we get the output power from the backward-running power by first applying (4a), i.e.
and then (4b), which together reads as This yields the result (21) given in the paper.

E Compound Poisson process for spontaneous emission
In order to model spontaneous emission occurring in the active medium, it is intuitively clear that some stochastic jump process will be necessary. In Section 2.3 of the paper, we state that it has been derived from a highly simplified model of amplifier chains that the number of emissions into a single mode (M = 1) is Bose-Einstein distributed with a mean number µ of emitted photons. The variance of this distribution is µ 2 + µ.
A jump process which can be designed to fit these two moments is the compound Poisson process given by where N t is a Poisson process with exponentially distributed waiting times. Y n are identically distributed random numbers independent of N t . The choice Y n ≡ 1 yields a Poisson process.
In view of the continuous-time nonlinear model from Section 3.3, we propose average waiting times of t RT /2 where an average of σ(λ j )N (t) (see Section 2.3) photons is emitted per longitudinal mode j. This information can be linked to (46) using Wald's formula such that [19] which suggests the rate parameter r = 2/t RT and E(Y 1 ) = σ(λ j )N (t). As mentioned above, the Bose-Einstein distribution fixes the variance of emitted photons after t RT /2 to be (σ(λ j )N (t)) 2 + σ(λ j )N (t). For (46), applying the Blackwell-Girshick equation yields [19] Var(X t ) = rt(E(Y 1 ) 2 + Var(Y 1 )) .
As a consequence, consistency with the Bose-Einstein distribution and E(Y 1 ) = σ(λ j )N (t) requires Therefore, Y n must be random numbers with both mean and variance equal to σ(λ j )N (t). This is satisfied by the Poisson distribution. We conclude that the stochastic process we use to model spontaneous emission is a compound Poisson process X t with rate parameter 2/t RT and Y n drawn from a Poisson distribution with intensity parameter σ(λ j )N (t).
The noise term C(t, λ j ) is finally obtained from the (rigorously non-existing) derivative with the factor t RT 2 for agreement with (20).

F Derivation of the discrete-time model
We stated in the paper that the discrete-time pulse-to-pulse dynamics can be approximately expressed as where f 1 represents the pumping of inversion for a duration t pump , f 2 the inversion and average intracavity power after the energy build-up for t bu and f 3,N and f 3,E the effect of coupling-out on inversion and pulse energy for t coup . We will now show how such expressions can be derived from the continuous-time nonlinear model given by (17,20) of the paper, which we repeat here as Since (49) cannot be solved analytically, we will proceed by solving simplified versions of (49) for the three sub-processes identified in Section 3.5 and then composing the respective solutions. At different stages of a switching cycle, different terms in (49) dominate, and hence, non-dominant ones will be neglected as a first simplification. Secondly, the spectrally distributed intracavity power P j (t) will be replaced by a spectral average P (t). To further simplify the functional structure of (49), nonlinear depletion terms in (49a) will be used only to first or second order. For the effect of coupling-out on the population, a shape of P (t) will be assumed in accordance with simulations to avoid solving (49b) with time-varying R(t). Finally, losses to P (t) will be included by an operator-splitting approach.
Pumping. For the first sub-process, the dominant terms in (49a) are the ones including P p and γ as the intracavity power P j is negligibly small for all j. Further assuming that 1. α p = 0 (no pump losses within the active medium), 2. LN dop N (t) (only a small fraction of all dopand atoms is in the excited state) and 3. P p = const. (the pump power is not varied during a switching cycle), we have the differential equation which allows the simple solution Energy build-up. The second sub-process is the energy build-up at high cavity quality R(t) = R max . N (t) and P j (t) start at f 1 (N m ) and u j,m := P (t m + t pump , λ j ), respectively. The dominant terms in (49) are now the depletion of N (t) due to P j (t) in (49a) and the amplification of P j (t) as well as losses induced by τ (t, λ j ) in (49b), which results in As mentioned above, the following steps will be taken to approximately solve (51): 1. P j (t) and 1/τ (t, λ j ) are replaced by spectrally averaged quantities P (t) and 1/τ := M j=1 (2α(λ j )L − ln(η(λ j )R max )) /(M t RT ) which describe the effective evolution of total energy within the cavity.
2. The exponential in (51a) is Taylor expanded to second order such that after cancellations, 3. To find an approximate solution for P (t), one needs to omit the quadratic terms in (52) which represent second-order effects as usually, σ(λ j )N (t) < 1 and α(λ j )L 1.
4. Losses are disregarded in (51b) and taken into account for P (t) by including 1/τ through an operatorsplitting approach.
5. The nonlinear depletion terms' effect on N (t) is approximated by inserting the expression for P (t) from the previous step into the simplified total rate equation from step 2.
6. The resulting expressions for N (t) and P (t) are finally augmented by fitting parameters s 1 and s 2 that mitigate modeling errors due to the simplifications performed above.
After approximation steps 1 and 2 and with the auxiliary quantities the sums over j in (51a) can be evaluated and one obtains If steps 3 and 4 (omission of nonlinear depletion and losses) are performed on (51b), both sides can be multiplied by λ j (1 − α(λ j )L) and summed over j. This results in For N (t) and P (t) starting the build-up at N (t m + t pump and u m := M j=1 u j,m /M with q 2 N 2 (t) ≈ q 0 ≈ 0 by step 3, (53,54) can be solved analytically, with the solution and the auxiliary parameters a N := 2hcA s a P := bΛt RT .
Using an operator-splitting approach to finish the fourth step, we refine the solution P lin (t) and re-incorporate the losses 1/τ from (51b) via Composing the solution of this differential equation (exponential decay) with (55b), we get As a next step (step 5), we can take the nonlinearity of (53) due to q 2 and q 0 into account by plugging P lin,OS (t) from (57) into (53) instead of P (t) and solving for N (t). Without theτ term, the approximate build-up solution for N (t) reads as Note that the solution exists for finiteτ as well. It is just a longer expression involving hypergeometric functions.
Simulations show that (58) fits depletion curves from more complex models relatively well. In a sixth step, the fit can be greatly improved, though, by incorporating a fitting parameter s 2 which accounts for nonlinearities beyond quadratic order in N (t) in (51a). If s 2 is placed into (58) in N (t m + t pump + t bu ) = 1 2q 2 tanh ln   a N N (t m + t pump ) + a P u m ) exp 2bq1 a N a P (a N N (t m + t pump ) + a P u m ) s 2 t bu a N N (t m + t pump ) + a P u m   × × q 2 1 − 4q 0 q 2 2q 1 + atanh q 1 + 2q 2 N (t m + t pump ) as an additional factor to t bu , the quasi-exponential fall in population can be shifted in its time of occurrence to fit some desired behavior. Shifting the time when N (t) starts to be significantly depleted via the same s 2 improves the accuracy of (57) as well. However, a second fitting parameter s 1 can be introduced to adjust how strongly the intracavity power rises, which is required due to disregarded nonlinear depletion terms from (51a) and the crude inclusion of losses via (56). The approximate build-up solution for P (t) is thus given by P (t m + t pump + t bu ) = u m   a N N (t m + t pump ) + a P u m a N N (t m + t pump ) exp − 2bq1 a N a P (a N N (t m + t pump ) + a P u m )) s 2 t bu + a P u m   s1 .
In total, the energy build-up contributes to populations N (t) and spectrally averaged intracavity powers P (t) in (48a,48b) according to with N (t m + t pump + t bu ) and P (t m + t pump + t bu ) from (59) and (60), respectively.
Coupling-out. Solving (49) with time-varying R(t) and hence time-varying τ (λ j , t) is not feasible. Therefore, for the third sub-process, 1. steps 1 and 2 from the energy build-up are repeated, resulting in (53).
2. It can be seen from simulations that for R(t) falling from R max to R min , P (t) falls from some maximal value to approximately zero within a similar timeframe as well, while population is being further depleted. Thus, the shape of P (t) during coupling-out is approximated by a line falling from P (t m + t pump + t bu ) to 0 within t coup .
These assumption entail that for the intracavity power, we can use the simple linear ansatz which we again insert into (53) and solve for N (t). Adding a third fitting parameter s 3 as a factor of P (t), the population after coupling-out (here without theτ term) is thus given by The purpose of s 3 is to account for the fact that P (t) is often not yet at its maximum at the moment the cavity is opened and instead continues to further rise for a short time.
Output pulse energy. We are also interested in finding a simple representation of the output energy as given by (48b). The most relevant contribution to this comes from output coupling. Hence, we again use the spectrally averaged intracavity power ansatz from (61) with (60) as well as N (t m + t pump + t bu ) from (58) and insert them into the deterministic version of the output power given by (21) from the paper. Finally, we convert this to the approximate m-th output energy by taking the integral over time while R(t) is falling, i.e. t m + t pump + t bu ≤ t < t m + t pump + t bu + t close , and when R(t) has fallen to its minimum value R min , i.e. between t close and t coup . The formula for this reads as  with R(t) = R max − (R max − R min ) t − (t m + t pump + t bu ) t close for t m + t pump + t bu ≤ t < t m + t pump + t bu + t close R(t) = R min for t m + t pump + t bu + t close ≤ t < t m + t pump + t bu + t coup = t m+1 P (t) = P (t m + t pump + t bu ) (1 − (t − (t m + t pump + t bu )) /t coup ) for t m + t pump + t bu ≤ t < t m+1 .
The integral can be carried out analytically, which symbolically endows us with f 3,E and thus completes the derivation of (48a,48b).