## 3 Existence: Nonlinear Modes Bifurcating from the Linear Limit  
As is customary, we focus on stationary soliton solutions of (1) of the form \(\Psi(x,y,z)=\psi(x,y)e^{i\mu z}\). Thus, one obtains the following stationary equation for \(\psi(x,y)\).  
\[\psi_{xx}+\psi_{yy}+U(x,y)\psi+\sigma|\psi|^{2}\psi=\mu\psi \tag{3}\]  
Figure 1: The plots show the spatial distribution of real (\(V\), left panel) and imaginary (\(W\), right panel) parts of the potential \(U\) with \(x_{0}=y_{0}=1.5\).  
Figure 2: The top left plot shows the spectrum of Schrödinger operator associated with the potential \(U\) in the complex plane (see also the text). Plots of the magnitude of the normalized eigenvectors for the three discrete eigenvalues \(\mu_{0}\) are shown in the other three plots.  
In [19], it is discussed that a continuous family of solitons bifurcates from each of the linear solutions in the presence of nonlinearity. In order to see this, let \(\mu_{0}\) be a discrete simple real eigenvalue of the potential \(U\) (such as one of the positive real eigenvalues in the top left plot of Figure 2). Now, following [19], expand \(\psi(x,y)\) in terms of \(\epsilon=|\mu-\mu_{0}|<<1\) and substitute the expression  
\[\psi(x,y)=\epsilon^{1/2}\left[c_{0}\psi_{0}+\epsilon\psi_{1}+\epsilon^{2}\psi _{2}+\dots\right] \tag{4}\]  
into Equation (3). This gives the equation for \(\psi_{1}\) as  
\[L\psi_{1}=c_{0}\left(\rho\psi_{0}-\sigma|c_{0}|^{2}|\psi_{0}|^{2}\psi_{0}\right) \tag{5}\]  
where \(\rho=\text{sgn}(\mu-\mu_{0})\) and  
\[|c_{0}|^{2}=\frac{\rho\langle\psi_{0}^{*},\psi_{0}\rangle}{\sigma\langle\psi_ {0}^{*},|\psi_{0}|^{2}\psi_{0}\rangle}. \tag{6}\]  
Here, \(\psi_{0}^{*}\) plays the role of the adjoint solution to \(\psi_{0}\).  
Thus, in order to find solutions of (3) for \(\sigma=\pm 1\), we perform a Newton continuation in the parameter \(\mu\) where the initial guess for \(\psi\) is given by the first two terms of (4). The bottom left panel of Figure 3 shows how the (optical) power \(P(\mu)=\int\int|\psi|^{2}dxdy\) of the solution grows as a function of increasing \(\mu\) for \(\sigma=1\), or as a function of decreasing \(\mu\) for \(\sigma=-1\) (from the linear limit). The first branch begins at the first real eigenvalue of \(U\) at \(\mu_{0}\approx 0.286\), the second branch at \(\mu_{0}\approx 0.487\), and the third branch begins at \(\mu_{0}\approx 0.785\). Plots of the solutions and their corresponding time evolution and stability properties are shown in the next section.  
As a general starting point comment for the properties of the branches, we point out that all the branches populate both the gain and the loss side. In the branch starting from \(\mu_{0}\approx 0.286\), all four \"wells\" of the potential of Figure 1 appear to be populated, with the lower intensity \"nodes\" being more populated and the higher intensity ones less populated. The second branch starting at \(\mu_{0}\approx 0.487\), as highlighted also in [19], possesses an anti-symmetric structure in \(x\) (hence the apparent vanishing of  
Figure 3: The bottom left plot shows the power of the solution \(\psi\) plotted in terms of the continuation parameter \(\mu\). The curves begin at the lowest power (i.e., at the linear limit) at the discrete real eigenvalues of approximately 0.286 (branch 1), 0.487 (branch 2), 0.785 (branch 3). Each power curve is drawn with its corresponding stability noted: a blue solid curve denotes a stable solution and a red dashed curve denotes an unstable solution. The other three plots track the maximum real part of eigenvalues \(v\) as a function of the continuation parameter \(\mu\): the red dashed line represents the max real part of eigenvalues that are real (exponential instability) while the blue dotted line tracks the max real part for eigenvalues that have a nonzero imaginary part (quartets); this case corresponds to oscillatory instabilities.  
the density at the \(x=0\) line). Both in the second and in the third branch, the higher intensity nodes of the potential appear to bear a higher intensity.