Nearly Hamiltonian dynamics of laser systems

The Arecchi-Bonifacio (or Maxwell-Bloch) model is the benchmark for the description of active optical media. However, in the presence of a fast relaxation of the atomic polarization, its implementation is a challenging task even in the simple ring-laser configuration, due to the presence of multiple time scales. In this Article we show that the dynamics is nearly Hamiltonian over time scales much longer than those of the cavity losses. More precisely, we prove that it can be represented as a pseudo spatio-temporal pattern generated by a nonlinear wave equation equipped with a Toda potential. The existence of two constants of motion (identified as pseudo energies), thereby, elucidates the reason why it is so hard to simplify the original model: the adiabatic elimination of the polarization must be accurate enough to describe the dynamics correctly over unexpectedly long time scales. Finally, since the nonlinear wave equation with Toda potential can be simulated on much longer times than the previous models, this opens up the route to the numerical (and theoretical) investigation of realistic setups.

The Arecchi-Bonifacio (or Maxwell-Bloch) model is the benchmark for the description of active optical media.However, in the presence of a fast relaxation of the atomic polarization, its implementation is a challenging task even in the simple ring-laser configuration, due to the presence of multiple time scales.In this Article we show that the dynamics is nearly Hamiltonian over time scales much longer than those of the cavity losses.More precisely, we prove that it can be represented as a pseudo spatio-temporal pattern generated by a nonlinear wave equation equipped with a Toda potential.The existence of two constants of motion (identified as pseudo energies), thereby, elucidates the reason why it is so hard to simplify the original model: the adiabatic elimination of the polarization must be accurate enough to describe the dynamics correctly over unexpectedly long time scales.Finally, since the nonlinear wave equation with Toda potential can be simulated on much longer times than the previous models, this opens up the route to the numerical (and theoretical) investigation of realistic setups.
The computation time can, however, be exceedingly long in the presence of a comparably fast polarization dynamics: a rather common feature encountered in lightamplifying devices [25,26] and optical networks [27][28][29][30].It is, therefore of paramount importance to overcome the difficulty posed by the presence of multiple time scales.Unfortunately, in spite of its extremely fast relaxation time, the atomic-polarization variable cannot be plainly adiabatically eliminated, since the resulting model exhibits unphysical divergencies.In mode-locking setups, the problem is not perceived, since they are typically based on the Haus Master equation [31].This is a phenomenological equation which does not include the atomic variables, consequently bypassing the unpleasant instabilities.Moreover, the Haus equation does not only lack a first-principle justification; it is also inaccurate in certain important cases: as noted in [32], the RNGH in-stability is not reproduced in the case of a plain ring laser, even when an equation for the population is added a posteriori.In fact, a more rigorous strategy requires a refined adiabatic elimination, as for instance proposed in [32,33].In this Article, we follow a similar strategy with a couple of differences: we start from the simpler and yet accurate "spaceless" delayed representation [34] with no periodic intra-cavity modulations.Time-delayed models had been already proposed for both the unidirectional propagation in a ring laser [35] and mode-locked lasers [36].However, they do not rigorously account for the atomic variables, as done in Ref. [34], which instead includes the polarization dynamics.By implementing a clean perturbative approach, the spaceless delayed model is progressively simplified, finally leading to a nonlinear wave equation with Toda potential (NWT).This result generalizes the observation made in the '80s that a single-mode laser characterized by a large polarization decay-rate can be viewed as a slowly relaxing Toda oscillator [37,38].The extension is nontrivial since the Hamiltonian here involves an infinity of modes and, furthermore, the "space" variable is actually a mixture of different time-scales.This latter property reveals also a substantial difference with the early work by Haken to describe multi-mode lasers [39], where a "free" energy was derived in the standard space-time representation.Last but not least, it is remarkable that the NWT model is to a large extent valid independently of the population decay rate and of the cavity losses which enter only to define the effective length and, indirectly, the corresponding NWT energy.
We start by briefly recalling the AB equations, since they are the fundamental reference: this is the model which takes into account the atomic-variable dynamics as from the theory of two-level atoms.For the sake of simplicity, we restrict ourselves to the resonant case and arXiv:2303.01061v1[physics.optics] 2 Mar 2023 thus assume that the atomic polarization P and the electric field F are real, as well as the population inversion D. In a comoving frame, the model can be written as, [34] where y is the scaled spatial variable (y The parameter γ is the ratio γ /γ ⊥ , where γ and γ ⊥ are the population and polarization decay rates, respectively, The time t is expressed in units of γ −1 ⊥ , while a is the pump parameter, controlling the amount of energy injected into the laser.The boundary condition is R being the mirror reflectivity and T the round-trip time across the cavity.
In Ref. [34], we have shown that this set of equations is generally well approximated by a spaceless delayed model, This is a slightly simplified version of the equations derived in [34]; we have, however, verified that they reproduce the regime generated by the original AB model (see below the discussion of numerical simulations).Here, P ( t) = P(y, t) , D = D(y, t) (the angular brackets denote a spatial average), and finally In a large fraction of laser devices, γ 1; these are the so-called class-B lasers.From now on, we consider this class of systems.It is well known that although the polarization relaxes quickly, its adiabatic elimination does not lead to a meaningful model [33].It is convenient to rescale the time, defining t = Γ t, where Γ = √ γ (accordingly T = ΓT ), and to introduce the rescaled pump-tolosses parameter Next, by following an observation made in Ref. [33] about the amplitude of the oscillations of P and D around their equilibrium value, we introduce a new set of variables, namely f = F/ √ I, u = P (I + 1) − F /(Γ √ I), and g = D(I + 1) − 1 /Γ.As a result, model (2) can be rewritten as We now adiabatically eliminate u by setting u = 0, Although this relationship is more accurate than the one obtained by setting the time derivative of P equal to 0, (we would have missed the last term in Eq. ( 5)), it is still not sufficiently accurate to reproduce all the relevant details.However, Eq. ( 5) is an important step forwards as it allows understanding the origin of the singularity of the limit Γ → 0, thereby identifying the backbone of the dynamical evolution.By plugging Eq. ( 5) into Eq.( 4), we obtain where we have introduced the smallness parameter and leave the correction terms in Eq. ( 7) unspecified.
Later we show that they are not relevant for the identification of the leading order of the dynamics.
We now transform the dependence on the single variable t into a spatio-temporal representation.This can be done by formally using the multiscale approach, which introduces slow and fast timescales εt and t, as described in [40][41][42].In such a case, the fast timescale plays the role of space and the slow timescale plays the role of time.However, we proceed here with a more phenomenological approach as in [43] (leading to the same result) by introducing (ζ, ξ), i.e., f (t) = f (ζ, ξ) and g(t) = ḡ(ζ, ξ), where ζ = t mod T is a pseudo-spatial variable, while ξ = t/T is a pseudo-temporal variable.This transformation can be intuitively interpreted as "wrapping" the time axis around a cylinder of circumference T in such a way that the longitudinal coordinate (the new time) increases by one unit after a full rotation -see more details in [41][42][43][44].It is easily seen that periodic boundary conditions hold in the space-time representation: f (ζ = 0, ξ) = f (ζ = T, ξ) (analogously for ḡ).Although this transformation was originally introduced and mainly used in the context of long delays, it is justified for arbitrary delay, like here.
This transformation is very useful since, from Eq. ( 6) If, moreover, we rescale the time axis, introducing, θ = εξ we can write, where the subscript denotes a derivative with respect to the corresponding variable.The change of variables from t to (ζ, θ) implies also that By substituting this into Eq.( 6) we obtain fθ + fζ = ḡ f , where we have neglected quadratic and higher order terms in ε.As terms of order O(ε) are absent in the first equation, we are entitled to neglect them in Eq. ( 7) too.As a result we obtain Eqs. (11,12) constitute the minimal description of the laser dynamics in the small Γ limit.The model is completed by recalling that the two fields satisfy periodic boundary conditions on the domain of size ζ L = T .A further, non-trivial simplification can be achieved by introducing q = ḡ/ √ I, s = ln f , and rescaling the variables As a result, the explicit dependence on I is removed, and it appears only in the definition of the "spatial" length L = T √ I.This equation can be given a simple physical interpretation by introducing the new time-like variable y = (τ + σ)/2 and the new space-like variable x = (3τ − σ)/2.In fact, the model can be written as a single second-order PDE (the subscript denotes a derivative with respect to the variable s).This is a nonlinear wave equation equipped with a Toda potential V (s) = e 2s − 2s − 1 [45].Notice that the validity of the Hamiltonian model depends only on the smallness of ε, irrespective of the presence of finite cavity losses (R < 1).Equation ( 14) generates a Hamiltonian dynamics.Being the NWT model translationally invariant in space, we expect energy and impulse to be both conserved.In the more physical (σ, τ ) description, the simplest independent conservation laws can be formulated in terms of the two densities With the help Eq.( 13) and recalling the definition of V (s), one can prove that the integral of h K is constant, where we have exploited the periodic boundary conditions.Analogously, By recalling that σ corresponds to the original "fast" time scale, h K can be read as a "kinetic" energy (KE) density.On the other hand, h P can be interpreted as a density of "potential" energy (PE), which accounts for the field fluctuations weighted according to the Toda potential.The existence of these two conservation laws is the primary reason why the adiabatic elimination of the polarization is such a delicate problem: the leading terms controlling the dynamics of the KE and PE must be captured correctly, and they require going one order beyond in the perturbation expansion.
In the limit Γ 1, the NWT model is the backbone of an amplitude equation for the laser dynamics.In a sense, this work complements the analysis of Ref. [46] performed for arbitrary Γ, but restricted to the region close to the RNGH threshold, where dynamical spatial structures emerge, and explains the singularity of the normal-form therein derived as a manifestaton of vanishing losses.Being the NWT model Hamiltonian, it is not structurally stable: arbitrarily small Γ generically brings in dissipation and amplification, via the herein neglected higher-order terms.Such terms affect the temporal evolution of the energy densities KE and PE, which are no longer constant.Nevertheless, over time-scales smaller than τ l ≈ 1/( √ Iε) perturbations should be negligible and the energies thereby conserved.
We now start validating he spaceless model ( 2) by comparing it with the original AB equations 1.In Fig. 1 the AB pattern generated after discarding a long transient is plotted (left) alongside with the outcome of the corresponding spaceless model (right).The pump value a = 7 is much above the RNGH threshold (linear stability shows that the instability arises above a ≈ 0.45 [17]) and long enough delay (τ = 124.1) to yield a multi-mode dynamics as also testified by the irregular spatial structure.It is transparent that the delayed model provides a fairly accurate representation of the laser dynamics and can, thereby, be considered as a reliable starting point.
Next, we compare the spaceless model with the Hamiltonian one (13).The delayed equations (2) have been simulated for the parameters reported in Fig. 1, with Spatio-temporal representation of the logarithm of the field amplitude, in grayscale from black (0.2) to white (2.3) for the left column, and from black (0.1) to white (3.0) for the right column.The two upper panels are obtained by integrating model ( 2); the lower ones by integrating (14).Left and right columns correspond to two different initial conditions.Parameters are as in Fig. 1, except for γ = 10 −4 /64 (Γ = 0.00125) and a total of τ = 117.9effective time units.γ = 10 −4 /64 (Γ = 0.00125) [47].We have selected a couple of instantaneous configurations in the stationary regime and used them as initial conditions to generate spatio-temporal patterns.The outcome is presented in the two upper panels of Fig. 2. We have then simulated the NWT model, rewriting it as, This equation can be seen as a first order PDE for the field q σ (σ, τ ), where q(σ, τ ) is determined via a spatial integration, under the condition that the spatial average of q is zero, as required by the boundary conditions.The NWT equation has been integrated by using a fourth-order finite-difference algorithm, with a time step δ = 0.005 and N = 8000 points to discretize the space.The resulting patterns, obtained by starting from the same initial conditions, are presented in the two lower panels of Fig. 2. The excellent agreement confirms the Hamiltonian-like nature of the underlying dynamics.
Since direct numerical simulations display a highfrequency instability, we have performed a spatial smoothing of amplitude 10 −4 every 100 time steps, which does not affect the conservation laws over the explored time scales.Whether this instability is the consequence of an ill posed Cauchy problem or not, it is immaterial, since higher-order perturbative terms are neverthe- less present, and the overall stability should be judged after including such corrections.This will be the task of future work.
Next, we have tested the behavior of the delayed model over yet longer time scales.The results obtained for different Γ values are reported Fig. 3.In the upper panel, the total energy density e tot = h K + h P is plotted vs. time (the other parameters are as in Fig. 2).There, we see that e tot is approximately constant and independent of Γ.The temporal fuctuations are a consequence of the neglected terms: the few jumps are quite likely induced by pattern selection mechanisms.
Interestingly, the smallness of the self-selected energy density (≈ 0.3/0.4) is suggestive of weak chaos, since in the linear limit with vanishing energy density the model is obviously integrable.This is consistent with the regularity of the patterns displayed in Figs.1,2.
In the lower panel of Fig. 3, we plot the ratio between KE and PE.It is nearly constant and close to 0.35.We have verified that the same holds true in a relatively broad range of a-pump values from 1 to 8. In equilibrium statistical mechanics, the ratio is 1.A refinement of our theory, including higher-order terms is necessary to explain the deviation.
In conclusion, we have shown that the ring-laser dynamics is well reproduced by the NWT model over the fast "spatial" scale t f ≈ (γ γ ⊥ ) −1/2 , and the longer Hamiltonian scale τ ≈ 1, which corresponds to t H ≈ T p R/[Γ a(1 − R)], where T p is the round trip time (all in physical units).Moreover, from the size of the neglected terms, we can predict the yet longer time scale t l ≈ t H τ l ≈ T p /[Γ 2 a(1 − R)], which is also the scale over which the two energies are expected to fluctuate.A fully quantitative analysis requires going one order beyond in the perturbative analysis: we leave this task to future work.The separation of time scales is better appreci-ated by considering a common Erbium laser configuration [48].Referring to effective two-level-system parameters (see, e.g., [20]) γ ⊥ = 2π × 10 12 s −1 , γ = 5 × 10 3 s −1 and T p = 100 ns, together with a = 7 and R = 0.95 here considered, we obtain t f ≈ 5.6 ns, t H ≈ 5.7 ms and t l = 3.6 s, with a span of almost 9 orders of magnitude.Notice that t H is much longer than the inverse of the field decay rate (≈ 10 µs), showing that the Hamiltonian character occurs over time scales when the field dynamics has been fully affected by its dissipative losses.
Remarkably, the Hamiltonian conservation laws (16,17) can be interpreted as the signature of a marginal stability of the system at this level of expansion.As a consequence, we have provided a rigorous ground to the common habit to move at least to the second order in the adiabatic elimination.
Finally, a further important advantage of the NWT equation: the presence of the much compressed time scale (the variable τ ) allows integrating the dynamical equations by using a time step much longer (by 3-4 orders of magnitude) than in the best current model (the PDE derived in Ref. [32]).
It is instructive to compare our results with those in Ref. [49], where light propagation is investigated in a Kerr medium.In that paper too, the authors start from a delay-algebraic system, to then derive a Hamiltonian amplitude-equation, which, in their case, has the form of a nonlinear Schrödinger equation.The similarity of the conclusions suggests that delay-algebraic models may possess general, still unearthed, properties.In fact, in spite of the analogies, the two models are substantially different: unlike in Ref. [49], our discussion of the propagation takes into account the role of polarization and does not require the long delay limit.
Our results pave the way to a series of routes that should be explored.Distributed losses arising from light propagation are often negligible; to what extent, their inclusion can change the overall scenario?Is the nearly Hamiltonian representation still valid in the presence of detuning as it happens in the single mode case [50]?Finally, what happens in semiconductor lasers, where the pump parameter a needs to be multiplied by a complex phase factor [51]? Minor qualitative changes are expected (see e.g.[32]), but the conclusion needs to be validated.

FIG. 1 .
FIG. 1. Spatio-temporal representation of the logarithm of the field amplitude in grayscale, from black (0.2) to white (2.3), as obtained from the integration of AB model (1) (left) and Eqs.(2) (right).Parameters are: a = 7, R = 0.95, roundtrip time T equivalent to L = 10 √ 139 ≈ 117.9, and γ = 10 −4 (Γ = 10 −2 ). 2 × 10 4 cells (delays) are shown, corresponding to τ = 124.1 time unit.The initial condition for the delayed model has been fixed by setting P and D equal to their value at the end of the cavity in the AB equations.