---
author:
- 'E. Ding'
- 'H. N. Chan'
- 'K. W. Chow'
- 'K. Nakkeeran'
- 'B. A. Malomed'
bibliography:
- 'ref.bib'
title: Exact states in waveguides with periodically modulated nonlinearity
---

Introduction
============

It is commonly known that optical spatial solitons arise in planar and bulk waveguides through the balance of the Kerr nonlinearity and transverse diffraction [@agrawal]. Modern fabrication technologies make it possible to create waveguides featuring spatially inhomogeneous nonlinearities that support novel classes of propagation patterns [@boris1]. In particular, spatially inhomogeneous waveguides with a defocusing nonlinearity, whose local strength grows toward the periphery, can support diverse species of fundamental and higher-order solitons, including vortices, necklace rings, vortex gyroscopes, *hopfions*, and complex hybrid modes [boris2,Lei\_Wu,wu, zhong,Radik,Yasha,hybrids]{}, as well as *localized dark solitons* [@Zeng]. Similar nonlinearity landscapes, featuring different growth rates of the local nonlinearity in opposite transverse directions, support strongly asymmetric bright solitons [@boris4]. Asymmetric solitons also appear spontaneously if the nonlinearity profile features a dual-well structure [@dual; @Nir]. Furthermore, a combination of the fast growing local strength of the defocusing nonlinearity with the usual $\mathcal{PT}$-symmetric gain-loss profile makes it possible to produce solitons that exhibit *unbreakable* $\mathcal{PT}$ symmetry [@unbreakable; @raju; @2D], which is essential for constructing robust solitons in such systems [@Demetri; @PTreview1; @PTreview2]. It is also relevant to mention that a combination of $\mathcal{PT}$-symmetric with competing nonlinearities supports spatiotemporal solitons [@ref1].

Considerable interest has also been drawn to models with uniform nonlinearity, either self-defocusing or focusing, and specially designed periodic potentials that support exact periodic wave solutions [Carr1,Carr2]{}. Although both the particular potentials and the corresponding exact periodic solutions are not generic, and the analysis of their stability can only be performed numerically, these models provide direct insight into the possibility to support periodic wave patterns by utilizing the interplay of periodic potentials and the ubiquitous cubic nonlinearity. Furthermore, nontrivial exact solutions serve as benchmarks which suggest the shape of generic solutions. The inverse problem, aimed at engineering waveguiding potentials adjusted to maintaining periodic waves with prescribed properties, is a physically relevant issue too [@Spain].

In this work, we introduce a model with a class of periodic modulations that represent spatially periodic *pseudopotentials* [@pseudo] induced by the local nonlinearity. This model admits exact solutions in the form of periodic wave patterns which, in the limiting case of an infinite modulation period, become bright solitons. Stability of these patterns is studied numerically. The same model can also be used to solve the inverse problem of engineering a nonlinearity-modulation profile needed to support a wave pattern with prescribed period and amplitude.

The Mathematical Model
======================

The light propagation in a planar waveguide with spatially modulated nonlinearity is described by the model that is based on the scaled nonlinear Schrödinger equation for the electromagnetic wave amplitude $%
\Psi (x,z)$, $$i\Psi _{z}+\Psi _{xx}+g(x)|\Psi |^{2}\Psi =0,  \label{eq:nls}$$where $x$ and $z$ are the transverse and longitudinal coordinates, respectively. The periodically-modulated nonlinearity profile is defined by $%
g(x)$, chosen as $$g(x)=\frac{\alpha }{\mathrm{dn}^{2}(x)}+\beta +\gamma ~\mathrm{dn}^{2}(x),
\label{eq:g}$$where $\alpha $, $\beta $, and $\gamma $ are real constants, and $\mathrm{dn}%
(x)$ is the standard Jacobi elliptic function with modulus $\sqrt{m}$ and period $2K$ ($K$ being the complete elliptic integral of the first kind). It is relevant to mention that, in the general case, the periodic inhomogeneity affects not only the local nonlinearity, but also the local refractive index, which would generate an additional term $U(x)\Psi $ in Eq. (\[eq:nls\]), with an effective spatially-periodic potential, $U(x)$. Nevertheless, specific experimental methods, such as resonant doping, make it possible to create waveguides in which the nonlinearity is affected by the periodic modulation, while the refractive index remains nearly constant [@boris1].

The same model, with propagation distance $z$ replaced by time $t$, represents the scaled Gross-Pitaevskii equation for the mean-field wave function of an atomic Bose-Einstein condensate (BEC), for which the periodic nonlinearity modulation can be induced by means of the Feshbach resonance in a spatially non-uniform magnetic or optical field. In particular, the necessary periodic profile of the magnetic field can be accurately implemented by means of the known technique based on the use of appropriately designed magnetic lattices [@magnetic]. Furthermore, the spatially periodic distribution of the local nonlinearity coefficient in BEC has been experimentally realized by means of an optical lattice [Feshbach]{}. A particular anharmonic profile corresponding to Eq. (\[eq:g\]) can be effectively approximated by a superposition of several harmonics of its Fourier decomposition, represented by the respective optical lattices. In a real experiment, the setting also includes an overall parabolic trapping potential. However, in many situations the characteristic scale of this potential is much larger than the period of the spatial modulation, which makes it possible to neglect the trapping potential while analyzing the effects of periodic lattices [@Heidelberg].

We look for a stationary solution to Eq. (\[eq:nls\]) in the form of a $\mathrm{dn}$-wave": $$\Psi (x,z)\equiv \psi (x)e^{-i\Omega z}=\frac{A_{0}\;\mathrm{dn}(x)}{1+b\;%
\mathrm{dn}^{2}(x)}e^{-i\Omega z},  \label{eq:ansatz}$$where $A_{0}$ and $b$ are real constants, and $-\Omega $ is the propagation constant (or the chemical potential, in the BEC model). Note that this ansatz is nonsingular under the conditions $b>-1$ or $b<-1/(1-m)$. Substituting it into Eq. (\[eq:nls\]) gives rise to the following system of equations for the three ansatz parameters $A_{0}^{2}$, $b$, and $\Omega $: $$\left\{ \begin{aligned} 2+A_{0}^{2} \; \alpha -6b(m-1)-m+\Omega =0 \; , \\
A_{0}^{2} \; \beta +2b^{2}(m-1)+2b(\Omega +3m-6)-2=0 \; , \\ 6b+A_{0}^{2} \;
\gamma +b^{2}(2-m+\Omega )=0 \; . \end{aligned}\right.  \label{eq:para}$$One can calculate any three parameters from this system for given values of the others. In particular, this allows one to address the above-mentioned inverse problem, aimed at determining the nonlinearity modulation profile (see Eq. (\[eq:g\])) needed for maintaining a particular wave pattern.

Stability Analysis
==================

The stability of the $\mathrm{dn}$-wave denoted by Eq. (\[eq:ansatz\]) is investigated by means of the standard linearization procedure [linstab1,linstab2]{}. Substituting $\Psi (x,z)=\left[ \psi (x)+u(x,z)\right]
e^{-i\Omega z}$ into Eq. (\[eq:nls\]), with the small complex perturbation defined as $u\left( x,z\right) \equiv R(x,z)+iI(x,z)$, we arrive at the linearized system, $$\left\{ \begin{aligned} \partial _{z}R &=&\left( -\Omega -g(x)\psi
^{2}(x)-\partial _{x}^{2}\right) I \; ,\\ \partial _{z}I &=&\left( \Omega
+3g(x)\psi ^{2}(x)+\partial _{x}^{2}\right) R \; . \end{aligned}\right.
\label{eq:linear}$$The stability of the $\mathrm{dn}$-wave is determined by substituting $%
\left\{ R(x,z),I(x,z)\right\} =\left\{ P(x),Q(x)\right\} \exp (\lambda z)$ into the above equations. The resulting problem for stability eigenvalue $%
\lambda $ is solved numerically using the finite-difference method. In particular, *modulational instability* of periodic states [agrawal]{} is accounted for by eigenvalues with $\mathrm{Re}(\lambda )>0$. Generic examples of stable and unstable $\mathrm{dn}$-waves are presented in Fig. \[fig:stab\]. The stability, as predicted by the calculation of the eigenvalue spectra, is corroborated by direct simulations of Eq. ([eq:nls]{}), using the Fourier transform in $x$ and a fourth-order Runge-Kutta algorithm in $z$.

The dependence of the stability of the $\mathrm{dn}$-waves on the system’s parameters can be explored with the help of numerical-continuation techniques [@Alan; @cont; @Doedel]. In particular, it is important to know how the stability is affected by varying the nonlinearity-modulation period $%
2K$, which is determined by the squared modulus, $m$. We fix $A_{0}$, $b$, and $\beta $ in Eq. (\[eq:para\]) as $A_{0}^{2}=1$, $b=-0.6$, and $\beta
=-1$, and determine the other parameters, *viz*., $\alpha $, $\gamma $, and $\Omega $, for each value of $m$. The results are summarized in Fig. \[fig:branch\]. In this case, the $\mathrm{dn}$-wave is found to be stable in the region of $0\leq m\leq 0.725$, where none of the eigenvalues in the spectrum has a positive real part (see the left panel). For $m>0.725$, the long-period $\mathrm{dn}$-waves are destabilized by at least one eigenvalue with $\mathrm{Re}(\lambda )>0$. The strongest instability is found at around $m=0.852$, where the $\mathrm{dn}$-wave is quickly destroyed by the instability. The evolution of the wave profiles at the onset of the instability, as well as the strongest-instability point, are also shown in the left panel. The *duty cycle* ($\mathrm{DC}$) of the modulation profile, i.e., the share of the region carrying a self-focusing nonlinearity per one period of $g(x)$, is shown in the right panel of Fig. [fig:branch]{}. The nonlinearity is entirely self-defocusing at $m<0.453$ where $\mathrm{DC}\equiv 0$. Once $m$ exceeds this threshold, the $\mathrm{DC%
}$ first increases to a maximum of $14.62\%$ at $m=0.755$, and then approaches zero in the long-wave limit of $m\rightarrow 1$ where the modulation period $2K$ becomes infinite. We have found that the mean value of $g(x)$ is always negative in the entire range of $m$ values, i.e., the nonlinearity is self-defocusing on average. It is worthy to note that the maximum of the $\mathrm{DC}$ roughly coincides with the onset of instability.

In the long-wave limit of $m\rightarrow 1$, the modulation profile (Eq. ([eq:g]{})) assumes the localized shape: $$g(x)=\alpha \cosh ^{2}x+\beta +\gamma \;\mathrm{sech}^{2}x\;.
\label{eq:sech}$$The left panel of Fig. \[fig:branch\] suggests that the instability of the -wave vanishes in this limit, with the corresponding stable exact solution (see Eq. (\[eq:ansatz\])) being a bright soliton: $$\Psi (x,z)=\frac{A_{0}\;\mathrm{sech}(x)}{1+b~\mathrm{sech}^{2}(x)}%
e^{-i\Omega z}\;.  \label{eq:soliton}$$The stability of the soliton solution is analyzed here only for the case where $\alpha =0$, to ensure that the localized modulation profile (Eq. ([eq:sech]{})) is not singular at $|x|\rightarrow \infty $ (nevertheless, the self-defocusing singularity with $\alpha <0$ may readily support robust self-trapped modes [boris2,Lei\_Wu,wu,zhong,Radik,Yasha,hybrids,boris4,dual,Nir]{}). In this situation, i.e., with $\alpha =0$ and $m\rightarrow 1$, system (\[eq:para\]) yields $$\Omega =-1\;,\;A_{0}^{2}=\frac{2(1+4b)}{\beta }\;,\;\gamma =-\frac{3b\beta }{%
1+4b}\;.  \label{eq:reduced}$$

The stability condition for the soliton pinned to the spatially modulated nonlinearity profile given by Eq. (\[eq:sech\]) with $\alpha =0$, is $%
\gamma >0$, as in that case the soliton is pulled to the local maximum of self-attraction. Equations (\[eq:soliton\]) and (\[eq:reduced\]) admit $%
\gamma >0$ in two cases:$$\beta >0,~0<-b<1/4 \; ;  \label{beta>0}$$$$\beta <0,~1/4<-b<1 \; .  \label{beta<0}$$In the former case, the nonlinearity is globally self-focusing, while in the latter one a finite self-focusing region (“defect") is embedded in a defocusing background.

An example of the latter situation is shown in Fig. \[fig:sech\], where the stability of the pinned bright soliton is confirmed by both the eigenvalue spectrum and direct simulations. In this case, $\beta =-2.8$ and $%
\gamma =3.6$,$\ g(x)$ being positive at $|x|<0.74$. In fact, stable solitons can be produced in a wide range of parameter values, as shown in Fig. [fig:soliton\_stab]{}. Unstable solitons are only found in the case of $\beta
>0>\gamma $ (represented by red dashed curves in the top panels). In this case, the self-focusing is stronger farther from the center, hence the soliton is repelled by the effective nonlinear potential. An example of such a nonlinearity profile is shown in the subplot in the top right panel of Fig. \[fig:soliton\_stab\].

Collisions between solitons play an important role in the study of their dynamics. In the case corresponding to condition (\[beta&gt;0\]), one may consider the collision of a free bright soliton, with inverse width $\eta $ and velocity (slope) $c$,$$\Psi _{\mathrm{free}}\left( x,z\right) = \sqrt{\frac{2}{\beta}} \eta\;%
\mathrm{sech}\left( \eta \left( x-cz\right) \right) e^{\left( \frac{i}{2}%
cx+i\left( \eta ^{2}-\frac{c^{2}}{4}\right) z\right)} \; ,
\label{free-bright}$$with a pinned soliton. An example of such a collision is displayed in Fig. \[fig:collision\]. The incident soliton captures the pinned one, merging with it into a single soliton which continues to move with original velocity. This outcome may find applications to the design of soliton=based data-processing schemes.

In the case of $\gamma \ll \beta $, i.e., $-b\ll 1/4$ (see Eqs. ([eq:reduced]{}) and (\[beta&gt;0\])), the collision can be considered by means of the perturbation theory [@RMP], which uses the exact result for the collision-induced soliton’s shift generated by the solution of the nonlinear Schrödinger equation. In this case, one may expects that the incident soliton will pass through, while the pinned one will start oscillating around the attractive nonlinear defect.

A completely novel situation arises in the case of Eq. (\[beta&lt;0\]), when the model admits a freely moving dark soliton far from the defect. Its collision with the pinned bright soliton will be governed by the repulsive interaction, which may lead to various outcomes, such as rebound of the incident dark soliton and destruction of the bright one through its dislodgment from the pinned position. These possibilities call for systematic numerical simulations of the collisions, which is a subject for a separate work.

The Manakov System
==================

Lastly, Eq. (\[eq:nls\]) can be generalized to a coupled system $$\left\{ \begin{aligned} i\Psi _{z}+\Psi _{xx}+g(x)\left( |\Psi |^{2}+|\Phi
|^{2}\right) \Psi &=0 \;,\\ i\Phi _{z}+\Phi _{xx}+g(x)\left( |\Psi
|^{2}+|\Phi |^{2}\right) \Phi &=0 \end{aligned}\right.   \label{eq:coupled}$$that describes the copropagation of light modes with orthogonal polarizations in a bimodal waveguide, under the Manakov’s condition that the self-phase- and cross-phase-modulation coefficients are equal [Manakov,Menyuk]{}, as well as a binary Bose-Einstein condensate composed of two hyperfine atomic states [@binary-BEC] (in the latter case, the relative nonlinearity is very close to the Manakov’s point). In the long-wave limit similar to Eq. (\[eq:sech\]) where $g(x)=\beta +\gamma ~%
\mathrm{sech}^{2}(rx)$, with $\gamma >0$ and parameter $r$ which determines the width of the attracting region, an exact solution of the coupled equations can be found in the form of a stable *symbiotic* dark-bright soliton complex [@symbio1; @symbio2]: $$\begin{aligned} \Psi (x,z) &=A_{0}\tanh (rx)e^{-i\Omega _{1}z} \;, \\ \Phi
(x,z) &=A_{0}\; \mathrm{sech}(rx)e^{-i\Omega _{2}z}\;, \end{aligned}
\label{eq:ms}$$with $$A_{0}=\sqrt{\frac{2}{\gamma }}\;r\;,\;\Omega _{1}=-\frac{2r^{2}\beta }{%
\gamma }\;,\;\Omega _{2}=-r^{2}-\frac{2r^{2}\beta }{\gamma }\;.$$In this case, the dark component $\Psi $ cannot exist without the interaction with the bright counterpart $\Phi $, and the background supporting the dark component is modulationally stable when $\beta <0$. An example of a stable dark-bright complex is displayed in Fig. [fig:darkbright]{}.

Conclusion
==========

In this work, we have studied the one-dimensional model for the wave transmission in a medium with a periodically-modulated local nonlinearity that is based on the Jacobi elliptic function. The model, which can be realized in optics and BEC [@ref2], admits both exact periodic solutions and bright solitons (in the long-wave limit). Stable solutions of these types provide a benchmark suggesting the shape of generic solutions that can be found numerically in the same model. The model also allows for the prediction of the modulation profile needed to support a particular periodic wave form with prescribed period and amplitude. The numerical analysis of the modulational stability has demonstrated that the periodic patterns can be unstable for sufficiently large periods. However, stability is retrieved in the limit of an infinite period which corresponds to bright solitons. In addition, we have found an exact solution for dark-bright soliton bound states in a similar two-component model that applies to periodically inhomogeneous bimodal planar optical waveguides and binary BEC. As an extension of the analysis, it will be relevant to study the periodic solutions, soliton solutions, localized structures [@ref3], and in particular bright-antidark soliton complex supported by this coupled model in detail. Also, as mentioned above, it may be interesting to systematically simulate collisions of a moving free bright or dark soliton with the pinned one.

Partial financial support has been provided by the Research Grants Council (Hong Kong) contract HKU 17200815.
