1. Introduction

Waveguide arrays (WGAs) are of significant technological interest in the broader photonics community due to their inherently nonlinear response to an intense applied electric field [1–3]. The rich dynamics of WGAs are fundamentally based upon the mode-coupling dynamics observed in optical waveguides [4], which was first theoretically investigated by Jones in 1965 [5] and later characterized experimentally [6], including in gallium arsenide (GaAs). Christodoulides was the first to proposed a theoretical model that demonstrated the self-localization of light in nonlinear WGAs [1], thus showing the possibility of observing discrete solitons in the WGAs. This original study led to a number of theoretical investigations aimed at understanding the possible photonic applications enabled by WGAs [7–9]. WGAs further found application in optical beam switching [10–13], Kerr lensing and mode-locking [3,14], pulse reshaping [15–18] , and semiconductor lasing [19–23]. Semiconductor WGAs are of great interest in state of the art technologies due to their nonlinear response to the applied electric field [24]. Recently, results of an experimental study [25] showed that incorporating a waveguide array in a fiber can result in mode-locking, and it was demonstrated that WGAs can be utilized as a fast saturable absorber to generate a discrete version of Kerr-lens mode-locking [26] in a passive erbium-doped fiber laser. As such we investigate the use of WGAs with quadratically nonlinear materials since the nonlinear responses of the WGAs can be invoked with much lower intensities of the applied electric field. Specifically, we show that mode-locking can be achieved in quadratic WGAs, with robust and stable solutions investigated in detail.

The nonlinear mode-coupling and spatial self-focusing properties of the WGAs made them ideal photonic devices for a discrete form of Kerr lensing in a mode-locked laser cavity [3,14]. In such applications the nearest-neighbor coupling of electromagnetic energy in the WGAs are governed by the following nonlinear difference-differential equation [1]

It has shown that using the WGAs, ultra-short pulses can be generated in passively mode-locked fiber laser cavities for both the normal and anomalous dispersion regimes [27,28] by balancing the mode-locked cavities chromatic dispersion and self-phase modulation. In [29], an averaged waveguide array based mode-locking (ML) model has been developed to describe the dynamics in a laser cavity

The reduction of the number of waveguides to this simplest set, along with other details of the simplification, can bee seen in [27].

Recently [30], a theoretical approach was proposed for spatial mode-locking of light-bullets in a physically constructed planar waveguide structure, in a cavity-less configuration, by a fundamental balance of spatial diffraction and self-phase modulation. For this configuration, the model Eq. (4) was modified by adding three-photon absorption and replacing chromatic dispersion with diffraction. Achievement of mode-locking in the developed WGA made it a potential photonic device to produce localized cavity-like solitons [30].

Both in the laser cavity [27,28] and cavity-less configuration [30], the WGAs were assumed to be formed from isotropic (Kerr) materials that are referred to as $\chi ^{(3)}$ materials, which is the dominant application studied [31–34]. However, in many optical systems, anisotropic quadratic ($\chi ^{2}$) materials, such as lithium niobate (LiNbO$_3$), potassium niobate (KNbO$_3$) or potassium titanyl phosphate (KTP), are applied [2,35–44]. $\chi ^{2}$ is the second order susceptibility that describes second harmonic generation first experimentally observed by Franken et al. [45]. After pulse shaping was demonstrated in a medium with quadratic nonlinearity [46], optical solitons theoretically predicted in [47,48]. In [49] and [50], it was demonstrated that similar to the centro-symmetric $\chi ^{3}$ materials, modulational instability and the self-defocusing phenomena are observed in $\chi ^{2}$ materials. In 1995, optical solitons in a quadratic bulk material were observed by Torruellas et al. [36] , and existence of the solitons in $\chi ^{2}$ waveguide were observed experimentally by Schiek et al. [38] in 1996. These original results have since been corroborated and extended in many follow-up experiments [51–56], demonstrating that quadratic solitons can exist in both the spatial and the temporal domains in waveguides or bulk materials [39–42].

It is well-known that the pulse dynamics in multidimensional non-resonant $\chi ^{2}$ materials cannot be generally described by nonlinear Schrödinger (NLS) based equations [57–60]. Indeed, these dynamics are governed by generalized NLS systems with coupling to a mean term (hereafter denoted as NLSM systems which are sometimes referred to as Benney–Roskes or Davey–Stewartson type equations) [61,62]. Benney and Roskes [61] first obtained NLSM equations in water of finite depth $h$ and without surface tension in 1969. In 1974, Davey and Stewartson [62] reached an equivalent form of the NLSM equations by investigation of the evolution of a 3D wave packet in water of finite depth. The integrability of NLSM systems was studied in 1975 by Ablowitz and Haberman [63] in the shallow water limit. In 1977, the results of Benney and Roskes was extended to include surface tension by Djordevic and Reddekopp [64]. Ablowitz et al. [57–59] derived from first principles NLSM type equations describing the evolution of the electromagnetic field in quadratic nonlinear media. Recently, in [65–67], it was demonstrated that optical wave collapse can be arrested in the NLSM system by adding an external potential (lattice) to the model.

The NLSM system is physically derived from an expansion of the slowly-varying wave amplitude in the first and second harmonics of the fundamental frequency and, a mean term that corresponds to the zeroth harmonic. This system describes the nonlocal-nonlinear coupling between a dynamic field that is related with the first harmonic and a static field that is related with the zeroth harmonic [68]. The general NLSM system is given by [37,58,59]

Given the NLSM model, we construct a theoretical model by combining the averaged model Eq. (2) and the NLSM Eq. (5) to characterize ML in quadratically nonlinear waveguide arrays (QWGAs). Thus we develop a ML theory for the QWGAs, by including three-photon absorption, bandwidth-limited saturating gain, linear attenuation, coupling between the waveguides and quadratic electro-optic effects to model ML dynamics which is capable of producing stable, 2D-soliton like solutions. After steady-state solutions are obtained, the linear and nonlinear stability of the proposed model are examined.

The paper is outlined as follows: In Sec. 2, model equations for ML in QWGAs is presented and steady-state solutions (fundamental solitons) of the model are obtained numerically. In Sec. 3, the mode-locking dynamics of the solutions are explored by direct numerical simulations and linear stability of the solutions are investigated by linear spectra. Results of the study are discussed in Sec. 4.

2. Governing equations

Mode-locking in quadratically nonlinear waveguide arrays (ML-QWGA) is characterized by the following (2+1) dimensional model

Here $A_0$, $A_1$, and $A_2$ are the envelope of the electric fields in the $0th$, $1st$ and $2nd$ waveguides respectively. The model combines the physics of Eqs. (4) and (5) along with the effects of three photon absorption [69]. The parameter $D$ is the diffraction coefficient in each of the slab arrays. The parameter $\beta$ is the strength of self-phase modulation, $p$ is proportional to the probability of three-photon absorption, $\gamma _j$ is the coefficient of linear attenuation in each waveguide, and $C$ is proportional to the strength of the coupling between waveguides. The saturable gain coefficient $g(t)$ accounts for the depletion of minority charge carriers at high optical intensities. The $g(t)\nabla ^{2}$ term prohibits the growth of high frequency spatial modes via the bandwidth limiting parameter $\tau$. The evolution of $0th$ waveguide is coupled to the $\phi (x,y)$ field. The coupling parameter $\rho$ describes the combined optical rectification and electro-optic effects and the material anisotropy is denoted by the parameter $\nu$. The parameter $\alpha$ represents quadratic polarization effect.

2.1 Numerical solutions using the pseudo-spectral renormalization (PSR) method

To calculate the dissipative steady-state solutions of the (2+1)D waveguide arrays mode-locking model (6), we use a modification of the pseudo-spectral renormalization (PSR) method [70] that is based on spectral renormalization method developed by Ablowitz and Musslimani [71]. The PSR method is outlined as follows.

Substituting the ansatz $A_j = U_j\left ( {x,y} \right ){e^{i\mu t}}$ into each waveguide in the model (6), where $\mu$ is the propagation constant, the following nonlinear eigenequation system is obtained

Applying the Fourier and inverse Fourier transformation to the terms, $\nabla ^{2} U_0={U_{0_{xx}}}+{U_{0_{yy}}}$, we get following equations

We multiply first equation in the system (10) by $w_0(x,y)$ and integrate over the entire space ($x$, $y$) to obtain renormalization parameter $|\lambda |$ as a root of following system

The iteration continues until the relative error $\lambda _{error}=|\lambda _{n+1}/\lambda _{n}-1|<10^{-10}$. It has been demonstrated that this algorithm converges rapidly for a wide-range of nonlinear PDEs [70–75]. Thus, the steady-state solutions (fundamental solitons) of the waveguide array mode-locking model (6) are obtained from a convergent iterative scheme. The initial starting point $w_0(x,y)$ is typically chosen to be a Gaussian and $w_1(x,y)=w_2(x,y)=0$. The algorithm usually convergences to the solution lower than 100 iterations when suitable parameter values are chosen for the model (6).

3. Mode-locked dynamics

The PSR method is utilized to compute a steady-state solution of the (2+1)D ML-QWGA model (6). Once the solution is obtained, it can be used to examine stability properties of the solitons. Unless otherwise specified, we set parameters to

We specifically set the coupling constant to $\rho = 0.5$ and the asymmetry parameter to $\nu =1.5$, typical values corresponding to the propagation of focused beams in potassium niobate (KNbO$_3$) [37].

With these parameters, the numerical convergence to the steady-state solution is shown in Fig. 1. 256 points were used in each of the spatial directions $x$ and $y$ on a square domain of length 80. Convergence to the steady-state was accomplished after 73 iterations and the solutions in waveguides 1 and 2 have inherited their form from waveguide 0, as expected. It is observed that, as demonstrated in previous studies [65,68,76], diffraction behavior of quadratic materials is anisotropic and as a result, allows the existence of elliptic solitons in the nonlinear regime. In the other words, the ML-QWGA model generates astigmatic steady-state solutions that do not possess radial symmetry. To investigate the level of astigmatism for the steady-state solutions of the waveguides, the following formulation is defined as a measure of astigmatism

When a solution is radially symmetric, $e=1$, whereas being comparatively wider along the $x$ and $y$ axes correspond to $e<1$ and $e>1$ cases respectively. Thus the solutions are elliptical if $e\neq 1$.

 

Fig. 1. Steady-state solution for (a) waveguide 0, (b) waveguide 1, and (c) waveguide 2 with parameter set that is given in Eq. (15).

Download Full Size | PDF

Contour plots of steady-state solutions are depicted in Fig. 2 for varied $\rho$ and $\alpha$ parameters. The first row of Fig. 2 shows that when $\rho =\alpha =0$, solutions are radially symmetric in all waveguides. The astigmatism along the $y$-axis is emphasized when the steady state solutions are obtained with larger $\alpha$ values (see the second and third rows). On the other hand, the solutions become more astigmatic along $x$-axis with larger coupling parameter $\rho$ (for a fixed $\alpha$), and when $\rho$ is increased to 1, the solitons become comparatively wider along $x$-axis (see the last row). It is also observed that when the parmaters $\rho$ and $\alpha$ are increased, astigmatism level of the $2nd$ waveguide is pronounced the most (see the second column).

 

Fig. 2. Contour image of steady state solutions for waveguide 0 (the first column), waveguide 1 (the second column) and waveguide 2 (the third column) when $\rho =0, \alpha =0$ (the first row), $\rho =0, \alpha =0.5$ (the second row), $\rho =0, \alpha =1$ (the third row), $\rho =0.5, \alpha =1$ (the fourth row) and $\rho =1, \alpha =1$ (the fifth row). All other parameters are taken as in Eq. (15) and $e$ shows the degree of astigmatism for each solution.

Download Full Size | PDF

3.1 Nonlinear evolution of the steady-state solutions

To investigate the dynamics of the steady-state solutions in waveguide arrays, we directly simulate Eq. (6) for long times. A finite-difference discretization scheme is used in the spatial domain and the solution is advanced in time with a fourth-order Runge-Kutta method for each waveguide. In Fig. 3, we plot the evolution of the steady-state solution, that is obtained in Fig. 1, from $t = 0$ to $t = 100$ with a numerical time-step of $dt = 0.01$ for the cross-section $y=0$. As can be seen from Fig. 3, mode-locked solutions can be achieved for the considered parameter regime that is given in Eq. (15).

 

Fig. 3. Nonlinear evolution of the steady-state solution in (a) $0th$, (b) $1st$, and (c) $2nd$ waveguide that is obtained in Fig. 1. The time step $dt = 0.01$ and $t\in [0,100]$. The cross-section of the solution for $y=0$ is shown.

Download Full Size | PDF

Additionally, we have observed that steady-state solutions of waveguide arrays can stay mode-locked for the parameter space that is depicted in Fig. 4. In this analysis, we have changed parameters in Eq. (15) one by one, when all other parameters are fixed, to find a stable mode-locking interval for the considered parameter and repeat the same procedure for each parameter. For instance, when all the parameters are fixed except $\rho$, we find steady-state solutions can stay mode-locked for $\rho \in [0,1.8]$. It has seen that the model requires a more careful selection of parameter values $C, E_0$ and $\tau$. It should also be noted that, although the $\nu$ parameter has a semi-infinite stability interval ($\nu >0$) mathematically , it is not practically applicable to real optical systems, since $\rho$ and $\nu$ parameters are fixed values depending on the type of material that is considered. The steady-state solutions give a great deal of insight into the underlying stable structures in the system. However, what is most important is the full dynamics of the ML-QWGA model without using the steady-state solutions as initial conditions. Specifically, whether or not the fundamental solutions spontaneously arise from white-noise initial solutions. Relaxing the steady-state assumption and repeating the computation with a very low amplitude white-noise initial condition, the system again evolves to a steady-state as shown in Fig. 5 for the parameters considered in Fig. 3 with $g_0=6$. We have observed that the stability interval of the parameters for white-noise initial condition are narrower than that of steady-state solutions. In other words, obtaining fundamental solitons and using them as initial conditions extends the stability region of the ML-QWGA model. For instance, stable mode-locking can be achieved for $g_0 \in [4.7,8.1]$ and $\rho \in [0,1.8]$ with the steady-state solution when it can be achieved for $g_0 \in [5.8,7.8]$ and $\rho \in [0,0.6]$ with white-noise initial condition.

 

Fig. 4. Stability intervals for the parameters of ML-QWGA model Eq. (6) when steady-state solutions are used as initial conditions of the evolution. The red dashes show the fixed parameter values that are given in Eq. (15).

Download Full Size | PDF

 

Fig. 5. Nonlinear evolution of the mode-locked pulses starting from white-noise initial condition in (a) $0th$, (b) $1st$, and (c) $2nd$ waveguide for the parameters given in Eq. (15) with $g_0=6$. The time step $dt = 0.01$ and $t\in [0,100]$. The cross-section of the solution for $y=0$ is shown.

Download Full Size | PDF

To examine the impact of the quadratic electro-optic and polarization effects on the dynamics of WGAs, we compare the stability regions of the parameters in ML-QWGA model Eq. (6) while $\rho =0, \nu =0, \alpha =0$ (without quadratic effects) and $\rho =0.5, \nu =1.5, \alpha =1$ (with quadratic effects) in Table 1. From Table 1 it can be seen that the stability regions of all parameters ($D, C, \beta, \gamma _0,\gamma _1,\gamma _2,E_0,\tau, p, g_0, \mu$) are extended when the model include quadratic effects ($\rho =0.5, \nu =1.5, \alpha =1$). Thus the existence of quadratic electro-optic effects in the model Eq. (6) improves the stability of generated pulse in waveguide arrays.

Table 1. Stability intervals for the parameters of ML-QWGA model Eq. (6) when $\rho =0, \nu =0, \alpha =0$ (without quadratic effects) and $\rho =0.5, \nu =1.5, \alpha =1$ (with quadratic electro-optic effects).

View Table

The behaviors of the mode-locking states in WGAs are significantly altered by the applied electric field (or the saturating gain $g(t)$), thus it is important to examine the impact of the gain on the pulses. $g(t)$ can be easily tuned in practical systems by the variation of the gain parameter $g_0$ and the saturated gain $E_0$. In Fig. 6, the peak amplitude of the steady state solutions are displayed as a function of $g_0$ and $E_0$.

 

Fig. 6. Peak amplitude of the pulse as a function of (a) the gain parameter $g_0$, and (b) the saturated gain $E_0$. Stable and unstable regions are shown by solid blue and red dotted lines, respectively. The insets in both panels highlight the bounds between stability and instability regions.

Download Full Size | PDF

Figure 6 shows that as the saturating gain increases via $g_0$ and $E_0$, the peak amplitude of the steady state solutions increase. Further, the stability (blue line) and instability (dotted red line) regions of solutions are determined in Fig. 6. It is seen that steady state solutions are stable (or mode-locked) when $g_0 \in [4.7,8.1]$ and $E_0 \in [0.94,1.18]$. When the gain parameter $g_0$ is larger than $8.1$, the pulses become unstable. From the previous studies, it is known that if the gain is increased, a Hopf bifurcation can occur, leading to the formation of breathing light bullet solutions [30,77].

To investigate the quadratic electro-optic and polarization effects on the dynamics of waveguides, peak amplitude of the field in waveguide 0 versus propagation time is plotted in Fig. 7 for varied values of $\rho, \nu$ and $\alpha$. It is observed that as the optical rectification parameter $\rho$ and anisotropy parameter $\nu$ are increasing, the peak amplitude of solution decreasing (see panel (a) anb (b)). On the other hand, increased quadratic polarization coefficient $\alpha$ amplifies peak amplitude of the solution (see panel (c)). In addition, we observe that peak amplitude of fields in the $1th$ and $2nd$ waveguides change in accordance with the peak amplitude of first waveguide when $\rho, \nu$ and $\alpha$ are varied. The dynamical behavior of the fields are considerably altered by the quadratic polarization parameter $\alpha$ that raises the intensity of the pulses in waveguides as increased and, the coupling parameter $\rho$ has a stronger influence than anisotropy parameter $\nu$ on the magnitude of considered fields.

 

Fig. 7. Evolution of the peak amplitude in the $0th$ waveguide ($|A_0|$) from $t=0$ to $t=100$. Steady-state solution is obtained (a) for varied values of $\rho$ when $\nu =1$ and $\alpha =1$; (b) for varied values of $\nu$ when $\rho =1$ and $\alpha =1$, and (c) for varied values of $\alpha$ when $\rho =1$ and $\nu =1$.

Download Full Size | PDF

3.2 Linear spectra of the steady-state solutions

To examine linear stability of the ML-QWGA model (6), we calculate the spectrum of linearization of the model about the steady-state solutions that are obtained with the PSR method. By denoting

Separating the steady-state solution $u_j$ and the perturbations $\tilde {u}_j$ into real and imaginary parts as follows

If the real part of the $\lambda$ is positive, the fundamental soliton is unstable. The eigenvalues of ${\mathbf A}$ can be calculated numerically with finite difference discretization of the spatial domain. Note that the matrix coefficients of ${\mathbf A}$ are given by

The $*$ notation denotes the integral $a_0* b_0=\int _{-\infty }^{\infty }{a_0(\tau ) b_0(\tau ) d\tau }$, which results from the nonlocal behavior given by the saturated gain dynamics [78].

The linear spectrum of solutions can be computed by evaluating the matrix ${\mathbf A}$. The spectrum of steady-state solution, that is obtained in Fig. 1 with parameter set in Eq. (15), is plotted in Fig. 8. It can be seen that spectrum of the solution is in the left-half plane and none of the eigenvalues has a positive real part, thus the considered steady-state solution is linearly stable in this parameter regime (see panel (a)). Also, we have observed that value of the propagation constant $\mu$ is critical for a stable linear stability spectrum and we determine stability interval of $\mu$ versus coupling coefficient $\rho$ and quadratic polarization parameter $\alpha$ (see panels (b) and (c)).

 

Fig. 8. (a) Linear spectrum of the steady-state solution that is obtained in Fig. 1 with parameters in Eq. (15). (b) Linear stability interval of the propagation constant $\mu$ versus optical rectification parameter $\rho$; and (c) Linear stability interval of $\mu$ versus quadratic polarization parameter $\alpha$; when all other parameters are fixed to values in Eq. (15).

Download Full Size | PDF

Furthermore, starting from a white-noise initial condition in each waveguide, we simulated linearized equation system (18) over a long distance. If the amplitude of pertubations ($\tilde {u}_0, \tilde {u}_1, \tilde {u}_2$) grow significantly in a finite time (or distance), the steady-state solution is considered to be linearly unstable. Otherwise, when the perturbation terms decay to zero, the original steady-state solution is linearly stable. Thus, it is mode-locked.

Linear evolution of the perturbations’ peak amplitude are displayed in Fig. 9 when the steady-state solution is obtained with parameters given in Eq. (15). It can be seen that peak amplitudes of perturbations ($\tilde {u}_0, \tilde {u}_1, \tilde {u}_2$) are decreasing significantly in a short time, and they decay to zero eventually ($t=30$). This result reveals that the steady-state solution in Fig. 1 is linearly stable. It is also observed that the stability regions shown in Fig. 8 is consistent with the evolution of linearized equation system (18).

 

Fig. 9. Linear evolution of the perturbations’ peak amplitude in (a) $0th$, (b) $1st$, and (c) $2nd$ waveguides for the steady-state solution that is shown in Fig. 1. The time step $dt = 0.01$ and $t\in [0,30]$.

Download Full Size | PDF

It is pointed out that while the steady-state solutions can stay mode-locked nonlinearly for $\rho \in [0, 1.8]$ and $\alpha \in [0.4, 4]$ (see Fig. 4), they have linearly stable spectra for $\rho \in [0, 1.6]$ and $\alpha \in [0.4, 2.2]$. This result reveals that obtaining linearly stable modes are more demanding than obtaining nonlinearly stable modes for the considered model.

4. Conclusion

We have proposed an ML-QWGA model by the addition of quadratic electro-optic effects and three-photon absorption to the averaged model (4) to describe ML in QWGA. Steady-state solutions of the ML-QWGA model are obtained by effective numerical methods. It has observed that these solutions are astigmatic due to anisotropic diffraction behavior of quadratic media. By direct simulation of the model, we have seen that both starting from steady-state solutions and white-noise initial conditions, stable mode-locking is achieved for a stability region of parameters. Stability intervals of all parameters in the ML-QWGA model have been determined comprehensively and it has seen that using steady-state solutions as initial conditions extends the stability region of the model. A linear stability analysis of the steady-state solutions are also thoroughly explored via the spectrum of linearization of the ML-QWGA model and the results show that when $\rho \in [0, 1.6]$ and $\alpha \in [0.4, 2.2]$ none of the eigenvalues has a positive real part, thus showing that the solutions in this region are linearly stable.

In addition, it has seen that pulse dynamics in the QWGA are significantly altered by the quadratic polarization parameter $\alpha$ that amplify amplitude of the fields as increased and, the coupling parameter $\rho$ has a stronger influence than anisotropy parameter $\nu$ on the magnitude of the fields. In conclusion, we have built the ML-QWGA model as an extension of averaged mode-locking model (4) and demonstrated the possibility of mode-locking of elliptic steady-state solutions in the QWGAs for a wide range of parameters.