Research Article  Open Access
Eli Shlizerman, Edwin Ding, Matthew O. Williams, J. Nathan Kutz, "The Proper Orthogonal Decomposition for Dimensionality Reduction in ModeLocked Lasers and Optical Systems", International Journal of Optics, vol. 2012, Article ID 831604, 18 pages, 2012. https://doi.org/10.1155/2012/831604
The Proper Orthogonal Decomposition for Dimensionality Reduction in ModeLocked Lasers and Optical Systems
Abstract
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 modelocked laser onto a lowdimensional 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 modelocking and the waveguide array modelocking 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.
1. Introduction
Ultrashort lasers have a large variety of applications, ranging from small scale problems such as ocular surgeries and biological imaging to largescale 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 modelocking [1, 2]. One of the most widely used modelocked 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–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 modelocked 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).
Modelocking can be achieved through a variety of mechanisms, for example, nonlinear polarization rotation [1] or waveguide arrays [2]. In the context of modelocking with nonlinear polarization rotation, the master modelocking equation proposed by Haus [3, 8–11], which is the complex GinzburgLandau equation with a bandwidthlimited gain, was the first model used to describe the pulse formation and modelocking dynamics in the ring cavity laser. The GinzburgLandau 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 patternforming systems. In the context of modelocking with waveguide arrays, the nonlinear mode coupling inherent in semiconductor waveguide arrays is exploited to produce the saturable absorber needed for modelocking [12–14]. The governing equation in this context is the waveguide array modelocking model (WGAML).
In both of these systems, a large number of stationary solutions are supported by this equation [14–16]. These stationary solutions can be categorized as singlepulse, doublepulse, and in general pulse solution depending on the strength of the cavity gain. The phenomenon for which a single modelocked 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–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 singlepulse to a doublepulse 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 modelocking 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–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 singlepulse 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–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 lowerdimensional 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–37]. In the context of modelocking, the method (or the associated method of moments) was successful in characterizing the singlepulse dynamics of the governing systems [37–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 KarhunenLoéve (KL) expansion, is a singularvalue decomposition (SVD) based technique often used to generate a lowrank, orthogonal basis that optimally (in an 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–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 modelocked pulse to the doublepulse solution observed in numerical simulations of both systems of modelocking. 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 soliton dynamics. This illustrates the key features of the method. The application of the POD reduction to the cubicquintic GinzburgLandau equation and the waveguide array modelocking model are then explored in Sections 4 and 5, respectively. Section 6 summarizes the results of the paper.
2. 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 is a vector of physically relevant quantities, and the subscripts and denote partial differentiation. The function captures the spacetime 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 is an orthogonal set of eigenfunctions, and we have assumed that the PDE is now described by a scalar quantity . The can be any orthogonal set of functions such that where is the Dirac function, and the notation 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 lowdimensional framework.
The POD method is related to the singular value decomposition (SVD). Computationally, the SVD is implemented as a builtin 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 . 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 samples with points per sample, then . The SVD factorizes the matrix into three matrices where , and , 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 . For , is an index that takes the values ; otherwise, it takes the values . The are referred to as the singular values of and are ordered such that . The matrices and are composed of the eigenvectors (rows of ) and (rows of ) of the covariance matrices and , respectively. As a result, the singular value decomposition allows the decomposition of the th row of into assuming that . Hence, the SVD returns a complete orthonormal set of basis functions for the rows of the data matrix . The elements of this basis are the vectors and are referred to as the POD modes. The key idea of this paper is embodied in (6). Specifically, the PODGalerkin method attempts to provide an accurate approximation of the with a system of ordinary differential equations.
One way to reduce the dimensionality of the matrix is to use a subset of the POD basis. The relative importance of the th POD mode in the approximation of the matrix is determined by the relative energy of that mode, defined as [40, 52] where the total energy is normalized such that . If the sum of the energies of the retained modes is unity, then these modes can be used to completely reconstruct . 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].
3. 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 as . Note that, for this example, time and space 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 where is an integer. We will consider the soliton dynamics with and , 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 ( 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 and solitons are demonstrated in Figures 1 and 2, respectively. During evolution, the soliton only undergoes phase changes while its amplitude remains stationary. In contrast, the 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.
(a)
(b)
(a)
(b)
The potentially obvious question to ask in light of our dimensionality reduction thinking is this: is the soliton dynamics really a 50 or 200 degreeoffreedom 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, lowdimensional basis in terms of nonlinear, soliton modes. Such a technique is analytically difficult to handle and requires indepth 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 and 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 lowdimensional 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 soliton being modeled exceptionally well by a single POD mode while the 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.
(a)
(b)
(c)
(a)
(b)
(c)
3.1. Soliton Reduction
To take advantage of the lowdimensional structure, we first consider the 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 . 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 . Note that the integration is over the computational domain and the denotes the complex conjugate and where .
The differential equation for given by (13) can be solved explicitly to yield where is the initial value condition for . To find the initial condition, recall that Taking the inner product with respect to gives Thus, the onemode expansion gives the approximate PDE solution This solution is the lowdimensional POD approximation of the PDE expanded in the best basis possible, that is, the SVD determined basis.
For the 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 onemode approximation.
(a)
(b)
(c)
3.2. Soliton Reduction
The case of the 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 twomode expansion takes the form where the and are simply taken from the first two columns of the 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 and . Recall that these two modes are orthogonal, thus, the resulting system of nonlinear equations results in where and the initial values of the two components are given by This gives a complete description of the twomode dynamics predicted from the SVD analysis.
The system (22) can be easily simulated with any standard numerical integration algorithm (e.g., fourthorder RungeKutta). Before computing the dynamics, the inner products given by , and must be calculated along with the initial conditions and .
Note that the twomode 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. Higherorder 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 lowdimensional structures for approximating PDE dynamics.
(a)
(b)
(c)
(d)
4. CubicQuintic GinzburgLandau’s Equation
As a more complicated example, we consider the master modelocking model that describes the evolution of the nonlinear polarization rotation modelocking system. The master modelocking model is known as the cubicquintic GinzburgLandau 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 crossphase modulations for the two orthogonal components of the polarization vector in the fast and slowfields, bandwidthlimited gain saturation, cavity attenuation, and the intensitydiscriminating element (saturable absorber) of the cavity. In particular, the equation takes the form where Here represents the complex envelope of the electric field propagating in the fiber. The independent variables and denote the propagating distance (number of cavity roundtrips) and the time in the rest frame of the pulse, respectively. 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 () 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 selfphase modulation of the field which results primarily from the nonlinearity (Kerr) of the refractive index, while ν (quintic modification of the selfphase 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 righthand side of the CQGLE), δ is the distributed total cavity attenuation. The gain , which is saturated by the total cavity energy (norm) , has two control parameters (pumping strength) and (saturated pulse energy). In what follows, we will assume without loss of generality that so that the cavity gain is totally controlled by . The parameter characterizes the parabolic bandwidth of the saturable gain. Unless specified otherwise, the parameters are assumed to be , , , , , , and . These parameters are physically achievable in typical ring cavity lasers [11].
To obtain a lowdimensional model that is able to describe the multipulsing transition of the CQGLE as a function of the gain parameter , we consider a combination of POD modes from different regions. The underlying idea is to use the dominating POD modes from different attractors of the CQGLE so that the resulting lowdimensional system carries as much information as possible and thus approximates the modelocking dynamics better [40]. We illustrate this by combining the first modes from the singlepulse region with the first modes from the doublepulse region.
Denote as the union of the two sets of orthonormal basis mentioned above, that is, Here and are the POD modes computed from the singlepulse and doublepulse region, respectively. The original field can be projected onto the elements of as where represent the elements of the combined basis , are the modal amplitudes, is the phase of the first mode, and denote the phase differences with respect to . One can obtain a lowdimensional system governing the evolutions of the modal amplitudes and the phase differences by substituting (28) into the CQGLE, multiplying the resulting equation by the complex conjugate of (), integrating over all 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 are not all orthonormal to each other. To address this issue, one can orthogonalize the combined basis prior to the projection to obtain a new basis , which can be achieved by the GramSchmidt algorithm. The reduced model derived from the new basis has a simpler form than the one derived directly from , since it does not have terms that are due to nonorthogonality of the basis functions. This reduced model will be called the () model, which indicates the number of modes taken from the singlepulse () and doublepulse () region.
We specifically consider the () model in this work [44]. The orthonormal basis obtained using the GramSchmidt algorithm is shown in Figure 7. By construction, the first mode in the new basis is identical to . The second mode contains a tall pulse which is identical to the first mode at 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 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.
(a)
(b)
(c)
(d)
4.1. Classification of Solutions of the () Model
In this section, we study the dynamics governed by the () model, 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 , , , and the three phase differences , , . To effectively characterize any periodic behavior in the system, it is customary to use the new set of variables for , so that the formal projection of the field is given by
Figure 8 shows the solution of the () model expressed in terms of the above variables and the reconstructed field at different cavity gain . At (top panel), the () model supports a stable fixed point with while the other variables are nearly zero; that is, the steady state content of mainly comes from , and; thus, a single modelocked pulse centered at the origin is expected. Increasing eventually destabilizes the fixed point, and a periodic solution is formed (middle panel) which corresponds to a limit cycle in the 7dimensional phase space. The reconstructed field has a tall pulse centered at the origin due to the significant contribution of (and hence ). The smallamplitude side pulses are resulted from the linear combination of the higherorder modes present in the system, and their locations are completely determined by the data matrix representing the doublepulse evolution of the CQGLE. Further increasing 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 and which are associated to and . Since the other variables in the system become negligible when , a doublepulse solution is formed. The () model is able to describe both the singlepulse and doublepulse 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 smallamplitude oscillations of the entire structure. In fact, these two features are also observed in the full simulations of the CQGLE. The () model also captures the amplification and suppression of the side pulses in the formation of the doublepulse solution in the CQGLE.
(a)
(b)
(c)
We construct a global bifurcation diagram in Figure 9 to characterize the transition between the singlepulse and doublepulse solutions of the () model. The diagram shows the total cavity energy of the reconstructed field as a function of the cavity gain , 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 singlemodelocked pulses of the () model. With increasing , the reconstructed modelocked pulse readjusts its amplitude and width accordingly to accommodate the increase in the cavity energy. The branch of the singlepulse solutions is stable until 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 . At , 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 doublepulse solution of the system and has higher energy than the singlepulse branch and the periodic branch.
The bifurcation diagram also illustrates that, in the region , the singlepulse and the doublepulse solutions are stable simultaneously. This bistability is not restricted to the () 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 singlepulse and doublepulse branches rather than the initial energy content. When the initial condition is selected such that it is “close" to the doublepulse branch (characterized by the Euclidean distance between them), the result of linear analysis holds, and a doublepulse modelocked state can be achieved. For random initial data that is far away from both branches, the system tends to settle into the singlepulse branch as it is more energetically favorable [12].
4.2. Comparison of Different () Models
To demonstrate that the () model is the optimal () model to characterize the transition from a singlemodelocked pulse into two pulses that occurs in the CQGLE, we consider the () model whose results are shown in the top row of Figure 10. This model can be derived by setting the variables and in the () model to zero. At a tall stable pulse is formed at together with a small residual at the origin. The entire structure has close resemblance to from the GramSchmidt procedure (see Figure 7). The tall pulse agrees with the singlepulse 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 .
The middle row of the same figure shows the numerical simulations of the () model, which is derived from the combination of the first POD mode from the singlepulse region with the first four modes from the doublepulse region. This model is able to reproduce all the key dynamics observed during the multipulsing phenomenon shown in Figure 8. The difference between the () model and the () model is the locations at which the bifurcations occur. In the previous section we showed that the singlepulse solution described by the () model loses stability at (Hopf bifurcation), and the subsequent periodic solution bifurcates into the doublepulse solution at . For the () model, these two bifurcations occur at and , respectively, which are closer to the values predicted using the full simulations of the CQGLE ( and ). The () model (bottom row of Figure 10) produces similar results as the () and the () models, with pulse transitions occurring at and 3.678. The comparison here shows that the third POD mode from the doublepulse 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 lowdimensional model.
5. Waveguide Array ModeLocking Model
The waveguide array modelocking model (WGAML) describes the modelocking process in a fiber waveguide system that combines the saturable absorption of the waveguide with the saturating and the bandwidthlimited gain of erbiumdoped fiber [12, 14]. The governing equations for such a system are given by with where , , and are the envelopes of the electric field in waveguides 0, 1, and 2, respectively. The normalized parameters used in the following sections are where controls the strength and the type (normal or anomalous) of the chromatic dispersion, β is the strength of the Kerr nonlinearity, is the saturating gain, is the bandwidth of the gain, is a linear loss in the th waveguide, and 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 . 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 .
In order to develop a lowdimensional model, we assume that where is the th POD mode from the th waveguide, and is the total number of retained POD modes. Substituting (35) into the WGAML in (32) and taking an inner product with each of gives the following system of ordinary differential equations where and The resulting system is a set of 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 , the must be chosen carefully. The process through which the were obtained is described in Section 5.1. In Section 5.2, we use these modes to determine the bifurcations of the reduced model.
5.1. 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 into one combined dataset. Each individual run consists of a shorttime but highresolution solution for with a sampling period of . The initial conditions used were a small perturbation away from the attractor for the fixed value of . 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 singlepulse, the breather, the chaotic, and the doublepulse 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 acquired from any run with a fixed value of . 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.
5.2. LowDimensional 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 in the relevant ranges for singlepulse, breather, chaotic, and doublepulse solutions using lowamplitude 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 considered. In this manner, the full and reduced dynamics of the WGAML starting from a coldcavity configuration can be compared. In addition to comparing the solutions at a single value of , 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 singlepulse, breather, chaotic, and doublepulse solutions reconstructed using a sixmode POD model, ordered from left to right. To obtain these solutions, the reduced model was evolved forward for 1000 units in starting with a lowamplitude whitenoise 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 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 at which these solutions are obtained. The stationary solutions of the POD model lose stability at lower values of than the stationary solutions of the WGAML. Further, the breather solutions in the POD model remain stable for larger values of than those in the WGAML. Due to the vastly smaller dimension of the reduced model and the range of 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 singlepulse and breather solutions of the WGAML. The chaotic and the doublepulse solutions were obtained using direct numerical simulations. While other solution branches may exist, they do not appear as attractors for whitenoise 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 doublepulse 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 symbols and denote the extrema of where . If no extrema exist, the mean value is denoted with a black dot. This bifurcation diagram reveals the sequence of bifurcations responsible for the multipulse transition. At the lowest values of , 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 singlepulse solution, which first appears from the saddlenode bifurcation that is labeled (a) in Figure 13 at . The value of represents the minimum gain needed for a singlepulse solution to exist in the WGAML. Although only the modes of the 0th waveguide are shown, these singlepulse 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 , there are two branches of singlepulse solutions that emanate from (a). The first branch is a highamplitude stable branch of solutions, and the other is a lowamplitude unstable branch. Following the lowamplitude unstable branch of solutions, another saddlenode bifurcation, which is labeled (b) in Figure 13, occurs. This bifurcation creates a stable branch of lowamplitude stationary solutions. Because their basin of attraction is small, these solutions are unlikely to appear for whitenoise 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 multipulse transition process.
On the other hand, the branch of largeramplitude solutions participates in the multipulse transition. This branch of singlepulse 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 FloquetFourierHill method [29]. For a range of , all other eigenvalues have negative real parts with the exception of the zero eigenvalue. As 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 singlepulse solution becomes unstable, but a stable limitcycle is generated. These limit cycles are the breather solutions of the WGAML.
As is increased, the breather solutions themselves will lose stability through two bifurcations in rapid succession. The first of these bifurcations is a perioddoubling bifurcation. After the perioddoubling bifurcation, the solution is still periodic in 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 perioddoubling bifurcation. The new stable perioddoubled limit cycles almost immediately lose stability through supercritical NeimarkSacker bifurcation. The NeimarkSacker bifurcation, also called a Torus bifurcation or secondary Hopfbifurcation, occurs when a complexconjugate 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 NeimarkSacker bifurcation indicates the presence of quasiperiodic solutions in the POD model. The period doubling and the NeimarkSacker 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 doublepulse solutions are gaining stability, and the system completes the transition from single to doublepulse 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 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 and serves as a route to chaos in .

6. Conclusion
The POD method, used here to construct the lowdimensional models for two different modelocking 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 modelocked lasers as a function of the parameters in the system, specifically the change in gain that occurs during the multipulse transition.
For the modelocking 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 is increased. It is shown that a fourmode reduction (also known as the () model) is adequate to capture the transition from a singlemodelocked pulse to two pulses through an intermediate periodic state observed in the CQGLE. It is found that the form of the singlepulse, doublepulse, 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 () model. When is increased, the fixed point representing the singlemodelocked pulse bifurcates into a stable limit cycle (periodic solution) through a Hopf bifurcation and then into a new type of fixed point (doublepulse solution). The pulse profiles predicted by the () model agree well with those generated by the full system. The model also demonstrates a bistability between the singlepulse solution and the doublepulse 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 modelocking 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 , 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 , the singlepulse 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 is further increased, it will eventually become chaotic through a series of bifurcations, initiated by a perioddoubling bifurcation and then almost immediately followed by a NeimarkSacker (Torus) bifurcation. The chaos in the direction is terminated by the doublepulse solution gaining stability and becoming the attractor for lowamplitude solutions. Here, the WGAML was modeled as a threewaveguide array. However, using the same methodology, a system of 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 modelocking 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 modelocking 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 multipulse 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 modelocking laser (maximal suppression of the multipulse instability), by searching for a set of parameters for which there is a maximal range of input cavity gain that supports stable singlepulse 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 lowdimensional 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.
Acknowledgments
J. N. Kutz acknowledges support from the National Science Foundation (NSF) (DMS0604700) and the US Air Force Office of Scientific Research (AFOSR) (FA9550090174).
References
 H. A. Haus, “Modelocking of lasers,” IEEE Journal on Selected Topics in Quantum Electronics, vol. 6, no. 6, pp. 1173–1185, 2000. View at: Publisher Site  Google Scholar
 J. N. Kutz, “Modelocked soliton lasers,” SIAM Review, vol. 48, no. 4, pp. 629–678, 2006. View at: Publisher Site  Google Scholar
 K. Tamura, H. A. Haus, and E. P. Ippen, “Selfstarting additive pulse modelocked erbium fibre ring laser,” Electronics Letters, vol. 28, no. 24, pp. 2226–2228, 1992. View at: Google Scholar
 H. A. Haus, E. P. Ippen , and K. Tamura, “Additivepulse modelocking in fiber lasers,” IEEE Journal of Quantum Electronics, vol. 30, pp. 200–208, 1994. View at: Google Scholar
 A. Chong, J. Buckley, W. Renninger, and F. Wise, “Allnormaldispersion femtosecond fiber laser,” Optics Express, vol. 14, no. 21, pp. 10095–10100, 2006. View at: Google Scholar
 A. Chong, W. H. Renninger, and F. W. Wise, “Properties of normaldispersion femtosecond fiber lasers,” Journal of the Optical Society of America B, vol. 25, no. 2, pp. 140–148, 2008. View at: Publisher Site  Google Scholar
 T. Brabec, C. Spielmann, P. F. Curley, and F. Krausz, “Kerr lens mode locking,” Optics Letters, vol. 17, pp. 1292–1294, 1992. View at: Google Scholar
 H. A. Haus, J. G. Fujimoto, and E. P. Ippen, “Structures for additive pulse modelocking,” Journal of the Optical Society of America B, vol. 8, p. 2068, 1991. View at: Google Scholar
 W. H. Renninger, A. Chong, and F. W. Wise, “Dissipative solitons in normaldispersion fiber lasers,” Physical Review A, vol. 77, no. 2, Article ID 023814, 2008. View at: Publisher Site  Google Scholar
 A. Komarov, H. Leblond, and F. Sanchez, “Quintic complex GinzburgLandau model for ring fiber lasers,” Physical Review E, vol. 72, no. 2, Article ID 025604, pp. 1–4, 2005. View at: Publisher Site  Google Scholar
 E. Ding and J. N. Kutz, “Operating regimes, splitstep modeling, and the Haus master modelocking model,” Journal of the Optical Society of America B, vol. 26, no. 12, pp. 2290–2300, 2009. View at: Publisher Site  Google Scholar
 J. N. Kutz and B. Sandstede, “Theory of passive harmonic modelocking using waveguide arrays,” Optics Express, vol. 16, no. 2, pp. 636–650, 2008. View at: Publisher Site  Google Scholar
 D. D. Hudson, K. Shish, T. R. Schibli et al., “Nonlinear femtosecond pulse reshaping in waveguide arrays,” Optics Letters, vol. 33, no. 13, pp. 1440–1442, 2008. View at: Publisher Site  Google Scholar
 J. L. Proctor and J. N. Kutz, “Theory and simulation of passive modelocking by use of waveguide arrays,” Optics Letters, vol. 30, no. 15, pp. 2013–2015, 2005. View at: Publisher Site  Google Scholar
 N. N. Akhmediev, A. Ankiewicz, and J. M. SotoCrespo, “Multisoliton solutions of the complex ginzburglandau equation,” Physical Review Letters, vol. 79, no. 21, pp. 4047–4051, 1997. View at: Google Scholar
 J. M. SotoCrespo, N. N. Akhmediev, V. V. Afanasjev, and S. Wabnitz, “Pulse solutions of the cubicquintic complex GinzburgLandau equation in the case of normal dispersion,” Physical Review E, vol. 55, no. 4, pp. 4783–4796, 1997. View at: Google Scholar
 A. Komarov, H. Leblond, and F. Sanchez, “Multistability and hysteresis phenomena in passively modelocked fiber lasers,” Physical Review A, vol. 71, no. 5, Article ID 053809, pp. 1–9, 2005. View at: Publisher Site  Google Scholar
 C. Y. Wang, W. Zhang, K. F. Lee, and K. M. Yoo, “Pulse splitting in a selfmodelocked Ti:sapphire laser,” Optics Communications, vol. 137, no. 1–3, pp. 89–92, 1997. View at: Google Scholar
 M. Lai, J. Nicholson, and W. Rudolph, “Multiple pulse operation of a femtosecond Ti:sapphire laser,” Optics Communications, vol. 142, no. 1–3, pp. 45–49, 1997. View at: Google Scholar
 M. Horowitz, C. R. Menyuk, T. F. Carruthers, and I. N. Duling, “Theoretical and experimental study of harmonically modelocked fiber lasers for optical communication systems,” Journal of Lightwave Technology, vol. 18, no. 11, pp. 1565–1574, 2000. View at: Publisher Site  Google Scholar
 B. G. Bale, K. Kieu, J. N. Kutz, and F. Wise, “Transition dynamics for multipulsing in modelocked lasers,” Optics Express, vol. 17, no. 25, pp. 23137–23146, 2009. View at: Publisher Site  Google Scholar
 Q. Xing, L. Chai, W. Zhang, and C. Y. Wang, “Regular, perioddoubling, quasiperiodic, and chaotic behavior in a selfmodelocked Ti:sapphire laser,” Optics Communications, vol. 162, no. 1, pp. 71–74, 1999. View at: Publisher Site  Google Scholar
 J. M. SotoCrespo, M. Grapinet, P. Grelu, and N. Akhmediev, “Bifurcations and multipleperiod soliton pulsations in a passively modelocked fiber laser,” Physical Review E, vol. 70, no. 6, Article ID 066612, pp. 1–11, 2004. View at: Publisher Site  Google Scholar
 P. Grelu, F. Belhache, F. Gutty, and J. M. SotoCrespo, “Phaselocked soliton pairs in a stretchedpulse fiber laser,” Optics Letters, vol. 27, no. 11, pp. 966–968, 2002. View at: Google Scholar
 S. Namiki, E. P. Ippen, H. A. Haus, and C. X. Yu, “Energy rate equations for modelocked lasers,” Journal of the Optical Society of America B, vol. 14, no. 8, pp. 2099–2111, 1997. View at: Google Scholar
 T. Kapitula, J. N. Kutz, and B. Sandstede, “Stability of pulses in the master modelocking equation,” Journal of the Optical Society of America B, vol. 19, no. 4, pp. 740–746, 2002. View at: Google Scholar
 T. Kapitula, N. Kutz, and B. Sandstede, “The Evans function for nonlocal equations,” Indiana University Mathematics Journal, vol. 53, no. 4, pp. 1095–1126, 2004. View at: Publisher Site  Google Scholar
 J. M. SotoCrespo, N. N. Akhmediev, and V. V. Afanasjev, “Stability of the pulselike solutions of the quintic complex GinzburgLandau equation,” Journal of the Optical Society of America B, vol. 13, no. 7, pp. 1439–1449, 1996. View at: Google Scholar
 C. R. Jones and J. N. Kutz, “Stability of modelocked pulse solutions subject to saturable gain: computing linear stability with the FloquetFourierHill method,” Journal of the Optical Society of America B, vol. 27, no. 6, pp. 1184–1194, 2010. View at: Publisher Site  Google Scholar
 S. K. Turitsyn and V. K. Mezentsev, “Dynamics of selfsimilar dispersionmanaged soliton presented in the basis of chirped GaussHermite functions,” JETP Letters, vol. 67, no. 9, pp. 640–646, 1998. View at: Publisher Site  Google Scholar
 S. K. Turitsyn, T. Schäfer, and V. K. Mezentsev, “Selfsimilar core and oscillatory tails of a pathaveraged chirped dispersionmanaged optical pulse,” Optics Letters, vol. 23, no. 17, pp. 1351–1353, 1998. View at: Google Scholar
 Y. S. Kivshar, T. J. Alexander, and S. K. Turitsyn, “Nonlinear modes of a macroscopic quantum oscillator,” Physics Letters, Section A, vol. 278, no. 4, pp. 225–230, 2001. View at: Publisher Site  Google Scholar
 A. Bondeson, M. Lisak, and D. Anderson, “Soliton perturbations: a variation principle for the soliton parameters,” Physica Scripta, vol. 20, no. 34, pp. 479–485, 1979. View at: Google Scholar
 D. Anderson, M. Lisak, and A. Berntson, “A variational approach to nonlinear evolution equations in optics,” Pramana, vol. 57, no. 56, pp. 917–936, 2001. View at: Google Scholar
 B. A. Malomed, “Resonant transmission of a chirped soliton in a long optical fiber with periodic amplification,” Journal of the Optical Society of America B, vol. 13, no. 4, pp. 677–686, 1996. View at: Google Scholar
 D. Anderson and M. Lisak, “Bandwidth limits due to mutual pulse interaction in optical soliton communications systems,” Optics Letters, vol. 11, p. 174, 1986. View at: Google Scholar
 W. Chang, A. Ankiewicz, J. M. SotoCrespo, and N. Akhmediev, “Dissipative soliton resonances,” Physical Review A, vol. 78, no. 2, Article ID 023830, 2008. View at: Publisher Site  Google Scholar
 E. Ding and J. N. Kutz, “Stability analysis of the modelocking dynamics in a laser cavity with a passive polarizer,” Journal of the Optical Society of America B, vol. 26, no. 7, pp. 1400–1411, 2009. View at: Publisher Site  Google Scholar
 B. G. Bale and J. N. Kutz, “Variational method for modelocked lasers,” Journal of the Optical Society of America B, vol. 25, no. 7, pp. 1193–1202, 2008. View at: Publisher Site  Google Scholar
 P. Holmes, J. Lumley, and G. Berkooz, Eds., Turbulence, Coherent Structures, Dynamical Systems and Symmetry, Cambridge University Press, Cambridge, UK, 1996.
 L. N. Trefethen and D. Bau, Numerical Linear Algebra, SIAM, 1997.
 D. S. Broomhead and G. P. King, “Extracting qualitative dynamics from experimental data,” Physica D, vol. 20, no. 23, pp. 217–236, 1986. View at: Google Scholar
 L. Sirovich, “Turbulence and the dynamics of coherent structures part i: coherent structures,” Quarterly of Applied Mathematics, vol. 45, p. 561, 1987. View at: Google Scholar
 E. Ding, E. Shlizerman, and J. N. Kutz, “Modeling multipulsing transition in ring cavity lasers with proper orthogonal decomposition,” Physical Review A, vol. 82, Article ID 023823, 2010. View at: Google Scholar
 M. O. Williams, E. Shlizerman, and J. N. Kutz, “The multipulsing transition in modelocked lasers: a lowdimensional approach using waveguide arrays,” Journal of the Optical Society of America B, vol. 27, no. 12, pp. 2471–2481, 2010. View at: Publisher Site  Google Scholar
 M. Ilak and C. W. Rowley, “Reducedorder modeling of channel flow using traveling POD and balanced POD,” in Proceedings of the 3rd AIAA Flow Control Conference, pp. 864–874, June 2006. View at: Google Scholar
 M. Ilak and C. W. Rowley, “Modeling of transitional channel flow using balanced proper orthogonal decomposition,” Physics of Fluids, vol. 20, no. 3, Article ID 034103, 2008. View at: Publisher Site  Google Scholar
 E. A. Christensen, M. Brøns, and J. N. Sørensen, “Evaluation of proper orthogonal decompositionbased decomposition techniques applied to parameterdependent nonturbulent flows,” SIAM Journal on Scientific Computing, vol. 21, no. 4, pp. 1419–1434, 1999. View at: Google Scholar
 P. del Sastre and R. Bermejo, Numerical Mathematics and Advanced Applications, Springer, 2006.
 A. Liakopoulos, P. A. Blythe, and H. Gunes, “A reduced dynamical model of convective flows in tall laterally heated cavities,” Proceedings of the Royal Society of London, vol. 453, no. 1958, pp. 663–672, 1997. View at: Google Scholar
 L. Sirovich, R. Everson, E. Kaplan, B. W. Knight, E. O'Brien, and D. Orbach, “Modeling the functional organization of the visual cortex,” Physica D, vol. 96, no. 1–4, pp. 355–366, 1996. View at: Google Scholar
 A. Chatterjee, “An introduction to the proper orthogonal decomposition,” Current Science, vol. 78, no. 7, pp. 808–817, 2000. View at: Google Scholar
 A. Dhooge, W. Govaerts, and Y. A. Kuznetsov, “MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs,” ACM Transactions on Mathematical Software, vol. 29, no. 2, pp. 141–164, 2003. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2012 Eli Shlizerman et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.