The Proper Orthogonal Decomposition for Dimensionality Reduction in Mode-Locked Lasers and Optical Systems

The onset of multipulsing, a ubiquitous phenomenon in laser cavities, imposes a fundamental limit on the maximum
energy delivered per pulse. Managing the nonlinear penalties in the cavity becomes crucial for increasing the energy and suppressing the multipulsing instability. A proper orthogonal decomposition (POD) allows for the reduction of governing equations of a mode-locked laser onto a low-dimensional space. The resulting reduced system is able to capture correctly the experimentally observed pulse transitions. Analysis of these models is used to explain the sequence of bifurcations that are responsible for the multipulsing instability in the master mode-locking and the waveguide array mode-locking models. As a result, the POD reduction allows for a simple and efficient way to characterize and optimize the cavity parameters for achieving maximal energy output.


Introduction
Ultrashort lasers have a large variety of applications, ranging from small scale problems such as ocular surgeries and biological imaging to large-scale problems such as optical communication systems and nuclear fusion.In the context of telecommunications and broadband sources, the laser is required to robustly produce pulses in the range of 200 femtoseconds to 50 picoseconds.The generation of such short pulses is often referred to as mode-locking [1,2].One of the most widely used mode-locked lasers developed to date is a ring cavity laser with a combination of waveplates and a passive polarizer.The combination of such components acts as an effective saturable absorber, providing an intensity discriminating mechanism to shape the pulse [3][4][5][6][7].It was first demonstrated experimentally in the early 90s that ultrashort pulse generation and robust laser operation could be achieved in such a laser cavity [3,4].Since then, a number of theoretical models have been proposed to characterize the mode-locked pulse evolution and its stability, including the present work which demonstrates that the key phenomenon of multipulsing can be completely characterized by a lowdimensional model with modes created by a proper orthogonal decomposition (POD).
Mode-locking can be achieved through a variety of mechanisms, for example, nonlinear polarization rotation [1] or waveguide arrays [2].In the context of mode-locking with nonlinear polarization rotation, the master mode-locking equation proposed by Haus [3,[8][9][10][11], which is the complex Ginzburg-Landau equation with a bandwidth-limited gain, was the first model used to describe the pulse formation and mode-locking dynamics in the ring cavity laser.The Ginzburg-Landau equation was originally developed in the context of particle physics as a model of superconductivity and has since been widely used as a prototypical model for nonlinear wave propagation and pattern-forming systems.In the context of mode-locking with waveguide arrays, the nonlinear mode coupling inherent in semiconductor waveguide arrays is exploited to produce the saturable absorber needed for mode-locking [12][13][14].The governing equation in this context is the waveguide array mode-locking model (WGAML).
In both of these systems, a large number of stationary solutions are supported by this equation [14][15][16].These stationary solutions can be categorized as single-pulse, doublepulse, and in general n-pulse solution depending on the strength of the cavity gain.The phenomenon for which a single mode-locked pulse solution splits into two or more pulses is known as multipulsing [1,2,17] and has been observed in a wide variety of experimental and theoretical configurations [12,[18][19][20][21][22][23][24].Specifically, when the pulse energy is increased, the spectrum of the pulse experiences a broadening effect.Once the spectrum exceeds the bandwidth limit of the gain, the pulse becomes energetically unfavorable, and a doublepulse solution is generated which has a smaller bandwidth.The transition from a single-pulse to a double-pulse solution can be via various mechanisms: abrupt jump from single pulse to double pulse, creation of periodic structures and chaotic/translational behavior.It is the aim of this manuscript to characterize the translational behavior in the nonlinear polarization rotation and WGAML systems and provide a framework for the study of the multipulsing phenomena when parameters are varied.
A number of analytical and numerical tools have been utilized over the past fifteen years to study the mode-locking stability and multipulsing transition.One of the earliest theoretical approaches was to consider the energy rate equations derived by Namiki et al. that describe the averaged energy evolution in the laser cavity [25].Other theoretical approaches involved the linear stability analysis of stationary pulses [12,[26][27][28].These analytical methods were effective in characterizing the stability of stationary pulses but were unable to describe the complete pulse transition that occurs during multipulsing.
Computationally, there is no efficient (algorithmic) way to find bifurcations by direct simulations, since solutions that possess a small basin of attraction are difficult to find.Alternative, local approach is to follow the single-pulse solution using continuation methods; however, to the limitation of the analytic tools, the transition region is difficult to characterize computationally even with the fastest and the most accurate algorithms available [29].Such difficulties has led to the consideration of reduction techniques that simplify the full governing equation and allow for a deeper insight of the multipulsing instability.One obvious method for a lowdimensional reconstruction of the dynamics is to use an appropriate eigenfunction expansion basis for characterizing the system.Such eigenfunction expansions represent a natural basis for the dynamics as have been demonstrated in various photonic applications [30][31][32].For nonlinear systems; however, it is not always clear what eigenfunction basis should be used nor is it clear that the optimal basis set can be found by such a technique.Nevertheless, it provides a powerful and effective method for recasting the dynamics in a lower-dimensional framework.Among other reduction techniques, the variational method was first used by Anderson et al. in the context of nonlinear evolution equations with nonhamiltonian dissipative perturbations [33,34].Since then, it has been widely used to study various aspects in nonlinear optics such as soliton resonance and pulse interactions [35][36][37].In the context of mode-locking, the method (or the associated method of moments) was successful in characterizing the single-pulse dynamics of the governing systems [37][38][39].
The limitation of the variational method is that the accuracy of the results depends largely on the approach used to describe the true dynamics.When multipulsing occurs, the form of the solution in the laser system changes significantly, and the prediction by the variational reduction is ambiguous.Unlike the variational approach [38,39], the lowdimensional model considered here is constructed using the method of POD [40].The proper orthogonal decomposition (POD), also known as principal component analysis (PCA) or the Karhunen-Loéve (KL) expansion, is a singular-value decomposition-(SVD-) based technique often used to generate a low-rank, orthogonal basis that optimally (in an L 2 -sense [41]) represents a set of data.Although the POD technique can be applied to generic data sets, it is often used on data obtained from systems with a physical, biological, and/or engineering application [40,42,43].The resulting set of the POD modes is often paired with the standard Galerkin projection technique in order to generate a set of ODEs for the mode amplitudes [40].This combination has proven to be very powerful and is used in a number of fields including nonlinear optics [44,45], turbulence [40], fluid flow [46][47][48][49], convective heating [50], and even neuroscience [51].As a result, the POD method provides a generic, analytic tool that can be used efficiently, unlike the direct methods.The resulting POD system is able to capture the whole bifurcation sequence from the single mode-locked pulse to the doublepulse solution observed in numerical simulations of both systems of mode-locking.The key idea behind the method is that there exists an "ideal" (optimal) basis in which a given system can be written (diagonalized) so that in this basis, all redundancies have been removed, and the largest variances of particular measurements are ordered.In the language being developed here, this means that the system has been written in terms of its principle components or in a proper orthogonal decomposition.As already mentioned, such techniques are commonly used as a tool in exploratory data analysis and for making predictive models.
This paper is outlined as follows: Section 2 describes the procedure for obtaining the POD modes from a given set of data.A simple example of the POD technique is applied to N-soliton dynamics.This illustrates the key features of the method.The application of the POD reduction to the cubicquintic Ginzburg-Landau equation and the waveguide array mode-locking model are then explored in Sections 4 and 5, respectively.Section 6 summarizes the results of the paper.

Proper Orthogonal Decomposition
Here we provide a short review of the POD method and describe its application to the derivation of reduced models for nonlinear evolution equations (see [40]).Abstractly, the POD method developed here can be represented by considering the following partial differential equation (PDE) system: where U is a vector of physically relevant quantities, and the subscripts t and x denote partial differentiation.The function N(•) captures the space-time dynamics that is specific to the system being considered.Along with this PDE are some prescribed boundary conditions and initial conditions.As a specific solution method for (1), we consider the separation of variable and basis expansion technique.
In particular, standard eigenfunction expansion techniques assume a solution of the form where the φ n (x) is an orthogonal set of eigenfunctions, and we have assumed that the PDE is now described by a scalar quantity u(x, t).The φ n (x) can be any orthogonal set of functions such that where δ jk is the Dirac function, and the notation (φ j , φ k ) = φ j φ * k dx gives the inner product.For a given physical problem, one may be motivated to use a specified set of eigenfunctions such as those special functions that arise for specific geometries or symmetries.More generally, for computational purposes, it is desired to use a set of eigenfunctions that produce accurate and rapid evaluation of the solutions of (1), that is, a Fourier mode basis and the associated fast Fourier transform.In the method developed here, optimal POD basis functions are generated from a singular value decomposition of the representative dynamics of the governing equations, thus, recasting the dynamics of the system into the best low-dimensional framework.
The POD method is related to the singular value decomposition (SVD).Computationally, the SVD is implemented as a built-in routine in many scientific software packages, such as MATLAB or NumPy.To generate a complete set of POD modes, a data set is compiled and represented as the matrix X.Each row of the matrix consists of a sample solution taken at a specific value of time, and the number of rows in the matrix is the number of samples taken at evenly spaced values in time.Therefore, if the data consists of m samples with n points per sample, then X ∈ C m×n .The SVD factorizes the matrix X into three matrices where U ∈ C m×m , V ∈ C n×n , and Σ ∈ R m×n , and the asterisk denotes the conjugate transpose.In a matrix form, the factorization in (4) is expressed as The matrix Σ is a diagonal matrix with nonnegative elements σ j .For m > n, j is an index that takes the values j = 1, . . ., n; otherwise, it takes the values j = 1, . . ., m.The σ j are referred to as the singular values of X and are ordered such that The matrices U and V are composed of the eigenvectors u i (rows of U) and φ j (rows of V * ) of the covariance matrices XX * and X * X, respectively.As a result, the singular value decomposition allows the decomposition of the kth row of X into assuming that m > n.Hence, the SVD returns a complete orthonormal set of basis functions for the rows of the data matrix X.The elements of this basis are the vectors φ j and are referred to as the POD modes.The key idea of this paper is embodied in (6).Specifically, the POD-Galerkin method attempts to provide an accurate approximation of the σ j u k j with a system of ordinary differential equations.
One way to reduce the dimensionality of the matrix X is to use a subset of the POD basis.The relative importance of the jth POD mode φ j in the approximation of the matrix X is determined by the relative energy E j of that mode, defined as [40,52] where the total energy is normalized such that n j=1 E j = 1.If the sum of the energies of the retained modes is unity, then these modes can be used to completely reconstruct X.Typically, the number of modes required to capture all of the energy is very large and does not result in a significant dimensionality reduction.In practice; however, it has been found that the matrix can be accurately approximated by using POD modes whose corresponding energies sum to almost all of the total energy.Then the POD basis that is used in the approximation is a truncated basis, where the POD modes with negligible energy are neglected.In practice, the truncated basis with energies that sum to 99% of the total energy is a plausible truncation.The advantage of using a truncated set of POD modes rather than any other set of modes is that the representation of the data generated by the POD modes is guaranteed to have a smaller least squares error than the representation of the data generated by any other orthonormal set of the same size [40].

Soliton Dynamics: A Simple Example Application
To give a specific demonstration of this technique, consider the nonlinear Schrödinger (NLS) equation with the boundary conditions u → 0 as x → ±∞.Note that, for this example, time t and space x are represented in a standard way versus the traditional moving frame reference of optics.If not for the nonlinear term, this equation could be solved easily in closed form.However, the nonlinearity mixes the eigenfunction components in the expansion (2) making a simple analytic solution not possible.It now remains to consider a specific spatial configuration for the initial condition.For the NLS, there are the set of special initial conditions corresponding to the soliton solutions

International Journal of Optics
where N is an integer.We will consider the soliton dynamics with N = 1 and N = 2, respectively.In order to do so, the initial condition is projected onto the Fourier modes with the fast Fourier transform (denoted by a hat).Rewriting (8) in the Fourier domain, that is, Fourier transforming, gives the set of differential equations where the spatial Fourier mode (k is the wavenumber) mixing occurs due to the nonlinear mixing in the cubic term.This gives the system of differential equations to be solved for in order to evaluate the NLS behavior.The dynamics of the N = 1 and N = 2 solitons are demonstrated in Figures 1 and 2, respectively.During evolution, the N = 1 soliton only undergoes phase changes while its amplitude remains stationary.In contrast, the N = 2 soliton undergoes periodic oscillations.In both cases, a large number of Fourier modes, about 50 and 200, respectively, are required to model the simple behaviors illustrated.
The potentially obvious question to ask in light of our dimensionality reduction thinking is this: is the soliton dynamics really a 50 or 200 degree-of-freedom system as implied by the Fourier mode numerical solution technique?The answer is no.Indeed, the inverse scattering transform ensures that there is an orthogonal, low-dimensional basis in terms of nonlinear, soliton modes.Such a technique is analytically difficult to handle and requires in-depth knowledge of the technique.On the other hand, with the appropriate basis, that is, the POD modes generated from the SVD, it can be shown that the dynamics is a simple reduction to 1 or 2 modes, respectively.Indeed, it can easily be shown that the N = 1 and N = 2 are truly low dimensional by computing the singular value decomposition of the evolutions shown in Figures 1 and 2, respectively.
Figures 3 and 4 demonstrate the low-dimensional nature explicitly by computing the singular values of the numerical simulations along with the modes to be used in our new eigenfunction expansion.What is clear is that for both of these cases, the dynamics is truly low dimensional with the N = 1 soliton being modeled exceptionally well by a single POD mode while the N = 2 dynamics is modeled quite well with two POD modes.Thus, in performing the expansion (2), the modes chosen should be the POD modes generated from the simulations themselves.In the following subsections, the dynamics of the modal interaction for these two cases are derived, showing that quite a bit of analytic progress can then be made within the POD framework.

N = 1 Soliton Reduction.
To take advantage of the lowdimensional structure, we first consider the N = 1 soliton dynamics.Figure 1 shows that a single mode in the SVD dominates the dynamics.Thus, the dynamics is recast in a single mode so that where now there is no sum in (2) since there is only a singlemode φ(x).Plugging this expansion into the NLS equation ( 8) yields the following: The inner product is now taken with respect to φ which gives where and the inner product is defined as integration over the computational interval so that (φ, ψ) = φψ * dx.Note that the integration is over the computational domain and the * denotes the complex conjugate and where (φ, φ) = 1.
The differential equation for a(t) given by ( 13) can be solved explicitly to yield where a(0) is the initial value condition for a(t).To find the initial condition, recall that Taking the inner product with respect to φ(x) gives Thus, the one-mode expansion gives the approximate PDE solution This solution is the low-dimensional POD approximation of the PDE expanded in the best basis possible, that is, the SVD determined basis.
For the N = 1 soliton, the spatial profile remains constant while its phase undergoes a nonlinear rotation.The POD solution (18) can be solved exactly to give a characterization of this phase rotation.Figure 5 shows both the full PDE dynamics along with its one-mode approximation.

N = 2 Soliton
Reduction.The case of the N = 2 soliton is a bit more complicated and interesting.In this case, two modes clearly dominate the behavior of the system.These two modes are now used to approximate the dynamics observed in Figure 2. In this case, the two-mode expansion takes the form where the φ 1 and φ 2 are simply taken from the first two columns of the U matrix in the SVD.Inserting this approximation into the governing equation (8) gives Multiplying out the cubic term gives the equation All that remains is to take the inner product of this equation with respect to both φ 1 (x) and φ 2 (x).Recall that these two modes are orthogonal, thus, the resulting 2 × 2 system of nonlinear equations results in International Journal of Optics where and the initial values of the two components are given by This gives a complete description of the two-mode dynamics predicted from the SVD analysis.The 2 × 2 system ( 22) can be easily simulated with any standard numerical integration algorithm (e.g., fourth-order Runge-Kutta).Before computing the dynamics, the inner products given by α jk , β jkl and σ jkl must be calculated along with the initial conditions a 1 (0) and a 2 (0).
Note that the two-mode dynamics does a good job in approximating the solution as demonstrated in Figure 6.
However, there is a phase drift that occurs in the dynamics that would require higher precision in both taking time slices of the full PDE and more accurate integration of the inner products for the coefficients.Indeed, the most simple trapezoidal rule has been used to compute the inner products, and its accuracy is somewhat suspected.Higher-order schemes could certainly help improve the accuracy.Additionally, incorporation of the third mode could also help.In either case, this demonstrates sufficiently how one would in practice use the low-dimensional structures for approximating PDE dynamics.

Cubic-Quintic Ginzburg-Landau's Equation
As a more complicated example, we consider the master mode-locking model that describes the evolution of the nonlinear polarization rotation mode-locking system.The master mode-locking model is known as the cubic-quintic Ginzburg-Landau equation (CQGLE).The model describes the averaged pulse evolution in a ring cavity laser subjected to the combined effect of chromatic dispersion, fiber birefringence, self-and cross-phase modulations for the two orthogonal components of the polarization vector in the fast-and slow-fields, bandwidth-limited gain saturation, cavity attenuation, and the intensity-discriminating element (saturable absorber) of the cavity.In particular, the equation takes the form where Here u represents the complex envelope of the electric field propagating in the fiber.The independent variables Z and T denote the propagating distance (number of cavity roundtrips) and the time in the rest frame of the pulse, respectively.D characterizes the dispersion in the cavity and is positive for anomalous dispersion and negative for normal dispersion.We will restrict ourselves to the case of anomalous dispersion (D > 0) which is consistent with experiments presented in [3,4].The rest of the parameters are also assumed to be positive throughout this paper, which is usually the case for physically realizable laser systems [11].The parameter γ represents the self-phase modulation of the field which results primarily from the nonlinearity (Kerr) of the refractive index, while ν (quintic modification of the self-phase modulation), β (cubic nonlinear gain), and μ (quintic nonlinear loss) arise directly from the averaging process in the derivation of the CQGLE [10,11].
For the linear dissipative terms (the first three terms on the right-hand side of the CQGLE), δ is the distributed total cavity attenuation.The gain g(Z), which is saturated by the total cavity energy (L 2 -norm) u 2 , has two control parameters g 0 (pumping strength) and e 0 (saturated pulse energy).In what follows, we will assume without loss of generality that e 0 = 1 so that the cavity gain is totally controlled by g 0 .The parameter τ characterizes the parabolic bandwidth of the saturable gain.Unless specified otherwise, the parameters are assumed to be D = 0.4, γ = 1, ν = 0.01, τ = 0.1, δ = 1, β = 0.3, and μ = 0.02.These parameters are physically achievable in typical ring cavity lasers [11].
To obtain a low-dimensional model that is able to describe the multipulsing transition of the CQGLE as a function of the gain parameter g 0 , we consider a combination of POD modes φ j from different regions.The underlying idea is to use the dominating POD modes from different attractors of the CQGLE so that the resulting low-dimensional system International Journal of Optics In the bottom figure, the phase of the pulse evolution at x = 0 is plotted for the full PDE simulations and the one-mode analytic formula (circles) given in (18).carries as much information as possible and thus approximates the mode-locking dynamics better [40].We illustrate this by combining the first m modes from the single-pulse region with the first n modes from the double-pulse region.
Denote S as the union of the two sets of orthonormal basis mentioned above, that is, Here φ (1)  j and φ (2)  j are the POD modes computed from the single-pulse and double-pulse region, respectively.The original field u can be projected onto the elements of S as where S j represent the elements of the combined basis S, a j are the modal amplitudes, θ 1 is the phase of the first mode, and ψ j denote the phase differences with respect to θ 1 .One can obtain a low-dimensional system governing the evolutions of the modal amplitudes a k and the phase differences ψ k by substituting (28) into the CQGLE, multiplying the resulting equation by the complex conjugate of S k (k = 1, . . ., m + n), integrating over all T and then separating the result into real and imaginary parts.The reduced system has, however, a complicated form due to the fact the elements in the combined basis S are not all orthonormal to each other.To address this issue, one can orthogonalize the combined basis S prior to the projection to obtain a new basis {Φ j } m+n j=1 , which can be achieved by the Gram-Schmidt algorithm.The reduced model derived from the new basis Φ j has a simpler form than the one derived directly from S, since it does not have terms that are due to nonorthogonality of the basis functions.This reduced model will be called the (m + n) model, which indicates the number of modes taken from the single-pulse (m) and double-pulse (n) region.
We specifically consider the (1 + 3) model in this work [44].The orthonormal basis {Φ j } 4 j=1 obtained using the Gram-Schmidt algorithm is shown in Figure 7.By construction, the first mode Φ 1 in the new basis is identical to φ (1)  1 .The second mode Φ 2 contains a tall pulse which is identical to the first mode at T = −6 and a small bump in the origin.In general, this tall pulse can be located on either the left or the right of the origin, depending on the data set X for which the POD modes are computed from.The other two modes have complicated temporal structures, and their oscillatory nature is mainly due to the orthogonalization.whose basis functions are shown in Figure 7.This fourmode model is 7 dimensional in the sense that the dynamic variables are the four modal amplitudes a 1 , a 2 , a 3 , a 4 and the three phase differences ψ 2 , ψ 3 , ψ 4 .To effectively characterize any periodic behavior in the system, it is customary to use the new set of variables x j = a j cos ψ j , y j = a j sin ψ j (29) for j = 2, 3, 4, so that the formal projection of the field is given by Figure 8 shows the solution of the (1+3) model expressed in terms of the above variables and the reconstructed field u at different cavity gain g 0 .At g 0 = 3 (top panel), the (1 + 3) model supports a stable fixed point with x 1 = 2.2868 while the other variables are nearly zero; that is, the steady state content of u mainly comes from Φ 1 , and; thus, a single mode-locked pulse centered at the origin is expected.Increasing g 0 eventually destabilizes the fixed point, and a periodic solution is formed (middle panel) which corresponds to a limit cycle in the 7-dimensional phase space.The reconstructed field has a tall pulse centered at the origin due to the significant contribution of x 1 (and hence Φ 1 ).The small-amplitude side pulses are resulted from the linear combination of the higher-order modes present in the system, and their locations are completely determined by the data matrix X representing the double-pulse evolution of the CQGLE.Further increasing g 0 causes the limit cycle to lose its stability and bifurcate to another type of stable fixed point.Unlike the first fixed point, this new fixed point has significant contributions from both x 1 and y 2 which are associated to Φ 1 and Φ 2 .Since the other variables in the system become negligible when Z → ∞, a double-pulse solution is formed.The (1 + 3) model is able to describe both the single-pulse and double-pulse evolution of the CQGLE.The two key features of the periodic solution of the CQGLE are explicitly revealed: (i) the existence of side pulses and (ii) the small-amplitude oscillations of the entire structure.In fact, these two features are also observed in the full International Journal of Optics Figure 7: Basis functions used in the (1 + 3) model.These functions are obtained by applying the Gram-Schmidt algorithm to the set S = {φ (1)  1 , φ (2)  1 , φ (2)  2 , φ (2)   3 }.
simulations of the CQGLE.The (1 + 3) model also captures the amplification and suppression of the side pulses in the formation of the double-pulse solution in the CQGLE.We construct a global bifurcation diagram in Figure 9 to characterize the transition between the single-pulse and double-pulse solutions of the (1 + 3) model.The diagram shows the total cavity energy of the reconstructed field u as a function of the cavity gain g 0 , which is obtained using MATCONT [53].In terms of the dynamic variables, the energy of the field is given by In the diagram, the branch with the lowest energy corresponds to the single-mode-locked pulses of the (1 + 3) model.With increasing g 0 , the reconstructed mode-locked pulse readjusts its amplitude and width accordingly to accommodate the increase in the cavity energy.The branch of the single-pulse solutions is stable until g 0 = 3.181 where a Hopf bifurcation occurs.The fluctuation in the cavity energy as a result of the periodic solution is small at first but increases gradually with g 0 .At g 0 = 3.764, the periodic solution becomes unstable by a fold bifurcation of the limit cycle.
In this case, a discrete jump is observed in the bifurcation diagram, and a new branch of solutions is reached.This new branch represents the double-pulse solution of the system and has higher energy than the single-pulse branch and the periodic branch.The bifurcation diagram also illustrates that, in the region 1.681 ≤ g 0 ≤ 3.181, the single-pulse and the doublepulse solutions are stable simultaneously.This bistability is not restricted to the (1 + 3) model as it also happens in the CQGLE as well [10,17].Given a set of parameters, the final form of the solution is determined by the basin of attraction of the single-pulse and double-pulse branches rather than the initial energy content.When the initial condition is selected such that it is "close" to the double-pulse branch (characterized by the Euclidean distance between them), the result of linear analysis holds, and a double-pulse modelocked state can be achieved.For random initial data that is far away from both branches, the system tends to settle into the single-pulse branch as it is more energetically favorable [12].

Comparison of Different (m + n) Models.
To demonstrate that the (1 + 3) model is the optimal (m + n) model to characterize the transition from a single-mode-locked pulse into two pulses that occurs in the CQGLE, we consider the (1 + 2) model whose results are shown in the top row of Figure 10.This model can be derived by setting the variables x 4 and y 4 in the (1 + 3) model to zero.At g 0 = 3 a tall stable pulse is formed at T = −6 together with a small residual at the origin.The entire structure has close resemblance to Φ 2 from the Gram-Schmidt procedure (see Figure 7).The tall pulse agrees with the single-pulse solution of the CQGLE quantitatively.When the structure loses its stability, a periodic solution is formed which persists even at unrealistically large values of cavity gain such as g 0 = 7.
The middle row of the same figure shows the numerical simulations of the (1 + 4) model, which is derived from the combination of the first POD mode from the singlepulse region with the first four modes from the double-pulse region.This model is able to reproduce all the key dynamics observed during the multipulsing phenomenon shown in Figure 8.The difference between the (1 + 4) model and the (1 + 3) model is the locations at which the bifurcations occur.In the previous section we showed that the singlepulse solution described by the (1+3) model loses stability at g 0 = 3.181 (Hopf bifurcation), and the subsequent periodic solution bifurcates into the double-pulse solution at g 0 = 3.764.For the (1 + 4) model, these two bifurcations occur International Journal of Optics at g 0 = 3.173 and g 0 = 3.722, respectively, which are closer to the values predicted using the full simulations of the CQGLE (g 0 ≈ 3.15 and g 0 ≈ 3.24).The (2 + 3) model (bottom row of Figure 10) produces similar results as the (1 + 3) and the (1 + 4) models, with pulse transitions occurring at g 0 = 3.177 and 3.678.The comparison here shows that the third POD mode φ (2)  3 from the double-pulse region is essential for getting the right pulse transition.Once this mode is included, further addition of the POD modes has only a minor effect on the accuracy of the low-dimensional model.

Waveguide Array Mode-Locking Model
The waveguide array mode-locking model (WGAML) describes the mode-locking process in a fiber waveguide system that combines the saturable absorption of the waveguide with the saturating and the bandwidth-limited gain of erbium-doped fiber [12,14].The governing equations for such a system are given by International Journal of Optics 13 with where A 0 , A 1 , and A 2 are the envelopes of the electric field in waveguides 0, 1, and 2, respectively.The normalized parameters used in the following sections are D, β, C, γ 0 , γ 1 , γ 2 , τ, e 0 = (1, 8, 5, 0, 0, 10, 0.1, 1), (34) where D controls the strength and the type (normal or anomalous) of the chromatic dispersion, β is the strength of the Kerr nonlinearity, g is the saturating gain, τ is the bandwidth of the gain, γ j is a linear loss in the jth waveguide, and C controls the strength of the evanescent coupling between adjacent waveguides.For a derivation and description of the governing equations, see the work of Proctor and Kutz [14] as well as Kutz and Sandstede [12].As shown by Kutz and Sandstede [12], the WGAML undergoes a region of chaotic translating behaviors for sufficiently large values of g 0 .This creates major difficulties in applying reduced dimensional models based on the POD [45].Therefore, in this paper we restrict ourselves to the set of solutions that are even up to a translation in T.
In order to develop a low-dimensional model, we assume that where φ (i) j is the jth POD mode from the ith waveguide, and N is the total number of retained POD modes.Substituting (35) into the WGAML in (32) and taking an inner product with each of φ (i) j gives the following system of ordinary differential equations where u, v = ∞ −∞ u * vdT and The resulting system is a set of 3N nonlinear evolution equations for the coefficient of the POD modes.This reduced set of equations is quite generic and is valid for any set of orthonormal functions, whether obtained through the POD or some other means.In order for this reduced model to be an accurate approximation of the full evolution equations over a range of g 0 , the φ (i) j must be chosen carefully.The process through which the φ (i) j were obtained is described in Section 5.1.In Section 5.2, we use these modes to determine the bifurcations of the reduced model.

POD Mode Selection.
In order to accurately capture the dynamics of the WGAML, the POD modes must be drawn from a particular set of data.Unlike with the CQGLE case, we do not combine POD modes drawn from different regimes through an orthonormalization process.Instead, we combine individual runs obtained at different values of g 0 into one combined dataset.Each individual run consists of a short-time but high-resolution solution for Z ∈ [0, 10] with a sampling period of ΔZ = 5 × 10 −3 .The initial conditions used were a small perturbation away from the attractor for the fixed value of g 0 .By combining these individual runs into a single data set, the SVD provides the set of POD modes that best captures the energy of the global dynamics through the transition from one to two pulses.
Figure 11 shows the first six POD modes generated using the methodology described in the preceding paragraph for solutions of the even WGAML model.Only the modes of the 0th waveguide are shown.However, there are two additional sets of modes, one for each of the other two waveguides.The data set includes information from the single-pulse, the breather, the chaotic, and the double-pulse solutions.These six modes capture over 99% of the modal energy, and 99.9% may be captured by including an additional four modes.The resulting modes do not resemble the modes that would be International Journal of Optics acquired from any run with a fixed value of g 0 .Indeed, the POD modes appear to be a nontrivial combination of the modes from all the runs at different gain values.These modes would be impossible to predict a priori even with knowledge about the solutions of the system.

Low-Dimensional Dynamics.
In this section, the reduced model is used to study the multipulsing transition of the WGAML.The reduced model is obtained by the Galerkin projection of WGAML model in (32) onto the global POD modes in Figure 11.This produces a system of differential equations of the form shown in (36).To determine the validity of the reduced model, we compare the full version of the evolution equations and the POD model dynamics for fixed values of g 0 in the relevant ranges for single-pulse, breather, chaotic, and double-pulse solutions using low-amplitude noise as the initial condition in both the full and reduced dynamics.The evolution considered is for a long period of time so that the attractor is reached for each value g 0 considered.In this manner, the full and reduced dynamics of the WGAML starting from a cold-cavity configuration can be compared.In addition to comparing the solutions at a single value of g 0 , we also compare the branches of singlepulse and breather solutions in both models.This allows the multipulsing transition to be studied in both systems.
The first row of Figure 12 shows the single-pulse, breather, chaotic, and double-pulse solutions reconstructed using a six-mode POD model, ordered from left to right.To obtain these solutions, the reduced model was evolved forward for 1000 units in Z starting with a low-amplitude white-noise initial condition.These results compare favorably with the same four regimes of the even WGAML dynamics shown in the second row of Figure 12, whose solutions were obtained by evolving for 200 units in Z starting from low amplitude even white noise.The reduced model qualitatively captures the dynamics and the profile of the solution.The primary difference is the value of g 0 at which these solutions are obtained.The stationary solutions of the POD model lose stability at lower values of g 0 than the stationary solutions of the WGAML.Further, the breather solutions in the POD model remain stable for larger values of g 0 than those in the WGAML.Due to the vastly smaller dimension of the reduced model and the range of g 0 for which the POD generated differential equations are valid, some disparity between the models is inevitable.
The advantage of the reduced model over the full dynamics is that the bifurcation diagram, including the stability and bifurcations of periodic solutions, can be explicitly calculated and categorized in the reduced model.The bifurcation software, MATCONT [53], was used to track and compute the stability of the single-pulse and breather solutions of the WGAML.The chaotic and the double-pulse solutions were obtained using direct numerical simulations.While other solution branches may exist, they do not appear as attractors for white-noise initial conditions and are not discussed here.
The first row of Figure 13 shows the bifurcation diagram of the reduced model including solution branches for the stationary, breather, chaotic, and double-pulse solutions.The branches marked in solid (blue) curves are linearly stable stationary solutions, and the branches marked in dashed (red) curves are linearly unstable stationary solutions.For the periodic solutions, the branches marked by black x symbols and denote the extrema of where If no extrema exist, the mean value A is denoted with a black dot.This bifurcation diagram reveals the sequence of bifurcations responsible for the multi-pulse transition.At the lowest values of g 0 , the only stable solution is the quiescent (trivial) solution because the gain is insufficient to overcome the cavity losses.The first nontrivial solution is the single-pulse solution, which first appears from the saddle-node bifurcation that is labeled (a) in Figure 13 at g 0 = 0.7229.The value of g 0 represents the minimum gain needed for a single-pulse solution to exist in the WGAML.Although only the modes of the 0th waveguide are shown, these single-pulse solutions distribute energy in all three of the waveguides, and the presence of energy in the other two waveguides is critical for the stabilization of these solutions.For larger values of g 0 , there are two branches of single-pulse solutions that emanate from (a).The first branch is a high-amplitude stable branch of solutions, and the other is a low-amplitude unstable branch.Following the low-amplitude unstable branch of solutions, another saddlenode bifurcation, which is labeled (b) in Figure 13, occurs.This bifurcation creates a stable branch of low-amplitude stationary solutions.Because their basin of attraction is small, these solutions are unlikely to appear for white-noise initial conditions.They are, however, known to exist in the full WGAML dynamics [29].This extra branch of stationary solutions loses stability through a Hopf bifurcation at the point labeled (c) in Figure 13.Therefore, the branch of lowamplitude stationary solutions does not play a role in the multi-pulse transition process.
On the other hand, the branch of larger-amplitude solutions participates in the multi-pulse transition.This branch of single-pulse solutions can be shown to be stable.Because we have assumed even solutions, the zero eigenvalue associated with translational invariance does not exist in the lowdimensional system.On the other hand, the zero eigenvalue associated with phase invariance persists.Although the values of the eigenvalues differ, the stability results agree qualitatively with the full results obtained with the Floquet-Fourier-Hill method [29].For a range of g 0 , all other eigenvalues have negative real parts with the exception of the zero eigenvalue.As g 0 increases, this branch becomes unstable through a supercritical Hopf bifurcation, which is shown in Figure 13 at (d).The Hopf bifurcation occurs when the eigenvalues of the linearized operator include a complexconjugate pair of eigenvalues with zero real part.As a result, the single-pulse solution becomes unstable, but a stable limit-cycle is generated.These limit cycles are the breather solutions of the WGAML.
As g 0 is increased, the breather solutions themselves will lose stability through two bifurcations in rapid succession.The first of these bifurcations is a period-doubling bifurcation.After the period-doubling bifurcation, the solution is still periodic in Z but with a larger period.In the reduced model, the Floquet multipliers can be computed explicitly.One of the multipliers was found to be −1, so the bifurcation that occurs is indeed a period-doubling bifurcation.The new stable period-doubled limit cycles almost immediately lose stability through supercritical Neimark-Sacker bifurcation.The Neimark-Sacker bifurcation, also called a Torus bifurcation or secondary Hopf-bifurcation, occurs when a complex-conjugate pair of the Floquet multipliers leave the unit circle.In this case, the limit cycle loses stability, but a torus that encloses the original limit cycle becomes stable.The Neimark-Sacker bifurcation indicates the presence of quasi-periodic solutions in the POD model.The period doubling and the Neimark-Sacker bifurcations can be seen in the top right panel of Figure 13 when the pair of extrema splits at the points labeled (e) and (f).MATCONT is unable to track branches of quasiperiodic solutions, but as shown in Figure 13, further bifurcations occur, and the ODE exhibits chaos.At the same time, the double-pulse solutions are gaining stability, and the system completes the transition from single-to double-pulse solutions.
The bifurcation diagram of the ODE shows remarkable qualitative agreement with the bifurcation diagram of the full WGAML dynamics illustrated in the second row of Figure 13.For the full system, linear stability has only been determined for the stationary solutions [29].For the stationary solutions, the types of bifurcations agree completely between the reduced and full models.Additionally, the value of g 0 agrees relatively well for all of the bifurcations as shown in Table 1.Although explicit stability results do not exist for the breather solutions, the same qualitative sequence occurs.The breather solutions lose stability around g 0 = 2.9 and serves as a route to chaos in Z.

Conclusion
The POD method, used here to construct the low-dimensional models for two different mode-locking systems, results as a robust tool for the study of the dynamics of nonlinear evolution equations describing various optical systems.In many cases where the observed dynamics are coherent but not trivial and exhibit various phenomena, the POD is the natural choice to construct a reduced and tractable model.If the model correctly reproduces the observed dynamics, it can then be studied using simpler computational methods and in special cases analytically in order to identify the bifurcations responsible for observed phenomena.The POD methodology can be easily modified for changes in the model.Such robustness of the method, allows one to study the operating regimes of mode-locked lasers as a function of the parameters in the system, specifically the change in gain that occurs during the multi-pulse transition.
For the mode-locking dynamics governed by the CQGLE model, the reduced system obtained via the POD method is able to capture the pulse transition that occurs in the laser when the cavity gain g 0 is increased.It is shown that a fourmode reduction (also known as the (1+3) model) is adequate to capture the transition from a single-mode-locked pulse to two pulses through an intermediate periodic state observed in the CQGLE.It is found that the form of the singlepulse, double-pulse, and the periodic solutions observed in the full numerics of the CQGLE can be accurately represented by the fixed points and limit cycles of the (1 + 3) model.When g 0 is increased, the fixed point representing the single-mode-locked pulse bifurcates into a stable limit cycle (periodic solution) through a Hopf bifurcation and then into a new type of fixed point (double-pulse solution).The pulse profiles predicted by the (1 + 3) model agree well with those generated by the full system.The model also demonstrates a bistability between the single-pulse solution and the double-pulse solution which also occurs in the CQGLE.The bifurcation diagram can be used as a guideline for achieving different types of solutions by choosing the initial condition appropriately.
For the mode-locking dynamics governed by the WGAML, the reduced model obtained by projecting the even WGAML model onto six global POD modes qualitatively reproduces the dynamics over the relevant range of gain values for which the transition from a single pulse to a double pulse occurs.In particular, for low values of g 0 , the model completely reproduces the dynamics of the singlepulse solution, including its bifurcations in a region of low amplitude solutions discovered earlier [29].With increasing gain g 0 , the single-pulse solution bifurcates to a limit cycle via a Hopf bifurcation, matching the study of WGAML model [12].Our analysis indicates that as the limit cycle grows when g 0 is further increased, it will eventually become chaotic through a series of bifurcations, initiated by a period-doubling bifurcation and then almost immediately followed by a Neimark-Sacker (Torus) bifurcation.The chaos in the Z direction is terminated by the double-pulse solution gaining stability and becoming the attractor for low-amplitude solutions.Here, the WGAML was modeled as a three-waveguide array.However, using the same methodology, a system of N-waveguide arrays could be studied.Only minor modifications to (36) are required to model a system with a different number of waveguides.
The parameters in the two reduced models constructed here indicate the physical setup of the mode-locking system, and their values were considered as fixed.By allowing variation of some of the parameters, the reduced models can serve as a framework for investigating optimal mode-locking and studying the nontrivial behavior observed in the multipulsing transition [21,24].Specifically, by using the described methodology, the application of standard continuation methods (such as MATCONT) and stability analysis to the reduced model allows for the constructing of a bifurcation diagram for the reduced model.Such a diagram indicates the impact of the parameters on the multi-pulse transition.It can reveal operating regimes for which multipulse transition is abrupt or results with a more complex behavior.Furthermore, it can be used to determine the optimal working regime for the mode-locking laser (maximal suppression of the multi-pulse instability), by searching for a set of parameters for which there is a maximal range of input cavity gain that supports stable single-pulse solutions.
In addition, models based on the POD are sometimes able to predict dynamics in regions outside of the data set, where the POD modes were computed, the low-dimensional model can be used to find interesting regions within the system without having to develop an appropriate approach.Once these regions have been identified, by altering the used data set, a more accurate representation of the region can be developed.In summary, the POD reduction gives a completely new and insightful way to explore the dynamics and bifurcations of a given system.

Figure 1 :
Figure 1: Evolution of the N = 1 soliton.Although it appears to be a very simple evolution, the phase is evolving in a nonlinear fashion, and approximately 50 Fourier modes are required to model accurately this behavior.

InternationalFigure 2 :
Figure 2: Evolution of the N = 2 soliton.Here a periodic dynamics is observed, and approximately 200 Fourier modes are required to model accurately the behavior.

Figure 3 :
Figure 3: Projection of the N = 1 evolution to POD modes.The top two figures are the singular values σ j of the evolution demonstrated in (1) on both a regular scale and log scale.This demonstrates that the N = 1 soliton dynamics is primarily a single-mode dynamics.The first two modes are shown in the bottom panel.Indeed, the second mode is meaningless and is generated from noise and numerical roundoff.

Figure 4 :
Figure 4: Projection of the N = 2 evolution to POD modes.The top two figures are the singular values σ j of the evolution demonstrated in (2) on both a regular scale and a log scale.This demonstrates that the N = 2 soliton dynamics is primarily a two-mode dynamics as these two modes contain approximately 95% of the evolution energy.The first three modes are shown in the bottom panel.

Figure 5 :
Figure 5: Comparison of (a) full PDE dynamics and (b) one-mode POD dynamics.In the bottom figure, the phase of the pulse evolution at x = 0 is plotted for the full PDE simulations and the one-mode analytic formula (circles) given in(18).

4. 1 .Figure 6 :
Figure6: Comparison of (a) full PDE dynamics and (b) two-mode POD dynamics.In the bottom two figures, the phase of the pulse evolution at x = 0 (middle panel) is plotted for the full PDE simulations and the two-mode dynamics along with the nonlinear oscillation dynamics of a 1 (t) and a 2 (t) (bottom panel).

Figure 11 :
Figure 11: POD modes taken from the combined data set for the 0th waveguide that are capable of qualitatively reproducing the dynamics observed in the full WGAML model.

Figure 12 :Figure 13 :
Figure 12: Top: surface and pseudocolor plots of the single-pulse, breather, chaotic, and double-pulse solutions computed for the finitedimensional model at g 0 = 1.5, 2.5, 3.495, and 3.5, respectively.Bottom: the same plots for the even WGAML model taken at g 0 = 2.3, 2.5, 3.1, and 3.2, respectively.The reduced model accurately reproduces the four behaviors observed.

Table 1 :
The values of g 0 associated with each bifurcation in the ODE and PDE.The labels correspond to the bifurcations of the ODE in Figure13.Values in parenthesis are estimated from the same figure.