# Dark Solitons and Grey Solitons in Waveguide Arrays with Long-Range Linear Coupling Effects
## 2 Model and Basic Equations  
In [17], we designed a scheme for an optics experiment to apply our model of long-range linear coupling waveguide arrays. The model is built based on an AlGaAs single-mode waveguide structure. The real-scale linear coupling parameter between different waveguides can be fitted with an exponential decay, hence the long-range linear coupling effect is introduced in the system.  
In our scheme, DSs can be described by the following equation, which is adapted from the discrete nonlinear Schrodinger equation [18; 19]:  
\[i\frac{\partial}{\partial z}u_{n}=\gamma|u_{n}|^{2}u_{n}-\sum_{m}C_{mn}u_{m}. \tag{1}\]  
Here, \(\gamma\) is the fixed nonlinear parameter of the system, where \(\gamma=1\) or \(\gamma=-1\) indicates that the system features self-focusing or self-defocusing nonlinearity, respectively. \(u_{n}\) is the field amplitude of the \(n\)-th mode. Because each waveguide is identical, for simplicity, the propagation constant \(\beta\) is absorbed into the phase of \(u_{n}\). \(z\) is the propagation distance along the waveguides, \(n\) is the number of waveguides, and the coefficient \(C_{mn}\) defines the coupling, which depends on the optical wavelength and the field overlap between m-th waveguide and n-th waveguide. Generally, Equation (1) can be expressed in matrix form as follows:  
\[i\frac{\partial}{\partial z}U=(C+V)U, \tag{2}\]  
where the matrix \(U\) and the elements of the matrix \(V\) are defined as \(U=\left(u_{1},\cdots,u_{N}\right)^{T}\) (where the superscript T indicates the transposition of the matrix and \(N\) is the number of waveguides) and \(V_{mn}=\gamma|u_{m}|^{2}\delta_{mn}\) (where \(\delta_{mn}\) is the Kronecker symbol), respectively.  
The total power of the guide mode in the system is given by \(P=\sum_{n}^{N}|u_{n}|^{2}\). In the model, we consider the matrix elements \(C_{mn}\) in Equation (1) to be given by  
\[C_{mn}=\begin{cases}c_{0}\exp(-j/d)&(j\
eq-1),\\ 0&(j=-1),\end{cases} \tag{3}\]  
where \(c_{0}\) is the control parameter, \(j=|m-n|-1\), and \(d\) is the normalized decay length, which describes the effective length of the coupling effects in the waveguides. When \(d\ll 1\), the system corresponds to a nearest-neighbour-coupled model. By contrast, when \(d\gg 1\), the system exhibits strong coupling effects. Equation (3) forms the linear coupling matrix \(C\), which represents the long-range linear coupling interactions among all lattice sites. In optics, when the separation between waveguides is sufficiently narrow, such higher-order cross-coupling effects can be induced [20].  
The relationship between the scaled parameters in Equation (1) and the real-scale parameters are given in [17] in detail.  
We assume that the soliton solutions for Equation (1) are written as  
\[u_{n}(z)=\phi_{n}e^{-i\mu z}, \tag{4}\]  
where \(\phi_{n}\) is the stationary solution and \(-\mu\) is the propagation constant, which is defined as  
\[\mu=\frac{U^{\mathrm{t}}(C+V)U}{P}. \tag{5}\]  
\(U\) is the solution of field amplitude in matrix form and \(U^{\mathrm{t}}\) is the conjugated matrix of \(U\). The stability of stationary solitons can be numerically determined by computing the eigenvalues for small perturbations or through direct simulations. The perturbed solution is given as \(u_{n}=e^{-i\mu z}(\phi_{n}+w_{n}e^{i\lambda z}+v_{n}^{*}e^{-i\lambda^{*}z})\). Substituting this solution into Equation (1) and linearizing yields the following eigenvalue problem:  
\[\left(\begin{array}{cc}C-\mu+2V&\Phi\\ -\Phi^{*}&-C+\mu-2V\end{array}\right)\left(\begin{array}{c}w\\ v\end{array}\right)=\lambda\left(\begin{array}{c}w\\ v\end{array}\right), \tag{6}\]  
where the elements of the matrix \(\Phi\) are defined as \(\Phi_{mn}=\gamma\phi_{m}^{2}\delta_{mn}\).  
The solution \(\phi\) is stable if all eigenvalues \(\lambda\) are real.