Localization of nonlinear excitations in curved waveguides

Motivated by the examples of a curved waveguide embedded in a photonic crystal and cold atoms moving in a waveguide created by a spatially inhomogeneous electromagnetic field, we examine the effects of geometry in a ‘quantum channel’ of parabolic form. Starting with the linear case we derive exact as well as approximate expressions for the eigenvalues and eigenfunctions of the linear problem. We then proceed to the nonlinear setting and its stationary states in a number of limiting cases that allow for analytical treatment. The results of our analysis are used as initial conditions in direct numerical simulations of the nonlinear problem and in this case localized excitations are found to persist. We found also interesting relaxational dynamics. Analogies of the present problem in context related to atomic physics and particularly to Bose–Einstein condensation are discussed.


Introduction
Localization phenomena are widely recognized as key to understanding the excitation dynamics in many physical context such as light propagation, charge and energy transport in condensed-matter physics and biophysics, and Bose-Einstein condensation of dilute atomic gases [1,2]. Recent advances in micro-structuring technology have made it possible to fabricate various low-dimensional systems with complicated geometry. Examples are photonic crystals with embedded defect structures such as microcavities, waveguides and waveguide bends [2]- [5], narrow structures (quantum dots and channels) formed at semiconductor heterostructures [6]- [8], magnetic nanodisks, dots and rings [9]- [11].
On the other hand, it is well known that the wave equation subject to Dirichlet boundary conditions has bound states in straight channels of variable width [12]- [14], and in curved channels of constant cross-section [15,16]. Spectral and transport characteristics of quantum electron channels [17] and waveguides in photonic crystals [3] are in essential ways modified by the existence of segments with finite curvature. The two-dimensional Laplacian operator supported by an infinite curve which is asymptotically straight has at least one bound state below the threshold of the continuum spectrum, as was recently proved in [18]. The appearance of an effective attractive potential in the wave equation is due to constraining quantum particles from higher to lower dimensional manifolds [19]- [21]. Curvature-induced bound-state energies and corresponding wave functions were studied in [22].
Until recently, there have been few theoretical and numerical studies of the effect of curvature on properties of nonlinear excitations. Nonlinear whispering gallery modes for a nonlinear Maxwell equation in microdisks were investigated in [23]; the excitation of whisperinggallery-type electromagnetic modes by a moving fluxon in an annular Josephson junction was shown in [24]. Nonlinear localized modes in two-dimensional photonic crystal waveguides were studied in [25]. A curved chain of nonlinear oscillators was considered in [26] and it was shown that the interplay of curvature and nonlinearity leads to a symmetry breaking when an asymmetric stationary state becomes energetically more favourable than a symmetric stationary state. Propagation of Bose-Einstein condensates in magnetic waveguides was experimentally demonstrated recently in [27]; single-mode propagation was observed along homogeneous segments of the waveguide, while geometric deformations of the microfabricated wires led to strong transverse excitations.
Motivated by the experimental relevance of the above-mentioned geometric deformations, in the present work we aim at investigating nonlinear excitations in a prototypical setup incorporating such phenomena. As our case example, we will examine an infinitely narrow curved nonlinear waveguide (channel) embedded in a two-dimensional linear medium (idle region). While the geometry that we propose is novel and different from the photonic lattices/crystals presently used in experiments, we expect that the intensities and parameters used in the latter setting (see e.g., the recent works of [28,29] and references therein as an example) will be similar to the ones relevant for excitation of nonlinear waves in the present context.
Our presentation will proceed as follows: in section 2, we set up the mathematical model of interest and examine its general properties and equations of motion. In section 3, we study the linear case and present its explicit solutions for the bound states, as well as for the corresponding eigenvalues. In section 4, we focus on the nonlinear problem, while in section 5 we supplement our analysis with numerical results. Finally, in section 6, we summarize our findings and present our conclusions. In the appendices, we present a discussion of a special example for the confining potential as well as some additional technical details on the solution of the linear problem and its Green's function formulation.

System and equations of motion
We consider a nonlinear Schrödinger model which is described by the Hamiltonian where ψ( r, t) is the complex amplitude function, r = (x, y), ∇ 2 = ∂ 2 x + ∂ 2 y , the coefficient A characterizes the nonlinearity of the medium, e.g., the nonlinear corrections to the refractive index of the photonic band-gap materials, self-interaction of the quasi-particles in the quantum channel, or interaction between cold atoms in Bose-Einstein condensates. V(x, y) is a scaled confining potential with a characteristic depth D (which has the dimension [L −2 ]) and width w. For example, in quantum channels formed in semiconductors, the depth D = 2m /h 2 , where is the energy difference between the quantum channel and the passive region and m is the effective mass of carriers; in photonic crystals with embedded defects D = 4π 2 (n 1 − n 0 )/λ 2 n 0 , where λ is the wave length, n 0 and n 1 are the refractive index of photonic crystal and waveguide, respectively (see e.g. [25]). To simplify calculations and separate the effects of finite curvature of waveguides from other effects, we restrict ourselves to the case when w 2 D 1, and wD = ν is finite. In this limit, the confining potential V(x, y) may be approximated as where the function y = f(x) gives the shape of the channel. Note that our approximation of the confining potential V(x, y) in the form (2) is consistent if the curvature κ of the channel satisfies the inequality: κ w 1.
From the Hamiltonian, we obtain the equation of motion in the form where the function F(|ψ| 2 ) is given by Equation (3) has as integrals of motion the Hamiltonian (1) and the (L 2 ) norm (referred to e.g., as the number of atoms in BEC or power in nonlinear optics) It is well known (see e.g., [30]) that if the tangent and the normal to the plane curve at a given point O ≡ (x 1 , y 1 ) are considered as the axes of a Cartesian coordinate system, then the equation of the curve in a neighbourhood of this point has the form where s is the arclength and κ the curvature of the curve at this point. Thus, if the curvature of the plane curve is not too large, one can represent it as a parabola In this case, in the frame of reference where the point O is chosen as the origin of coordinates, equation (3) takes the form where R = 1/κ is the maximum radius of curvature of the curve. Note that cold atoms in BEC experiments may experience the potential when they interact with two coherent light beams of which one is a diverging spherical wave and the other is a plane wave (see appendix A for details). In this case, equation (8) is the Gross-Pitaevski equation for BEC in a curved optical guide. It is convenient to use the parabolic coordinates The coordinate lines are two orthogonal families of confocal parabolas. These lines are given by The variable u is allowed to range from −∞ to ∞, whereas v is positive. Introducing equations (10) into equation (8) and using the properties of the delta-function (see e.g., [31]), we obtain Using the Fourier transform with respect to t, where the bars denote the Fourier-transformed quantities, one can represent equation (12) in the form Equation (12) can, in turn, be expressed in the form of the integral equation where the Green function G(u, v; u , v ) satisfies the equation and has the form (see appendix B for details) Here, where H n (z) is the Hermite polynomial [32], a n = (ω 1/2 /R π) 1/4 (2 n n!) −1/2 is the normalization constant, and with V(a, x) and U(a, x) being the Weber parabolic cylinder functions [32].
It can then be seen from equations (15) and (17) that the Fourier-transformed channel wave functionφ satisfies the equation or equivalently the equation The wave functionψ(u, v, ω) may be represented in terms of the channel wave functionφ(u, ω) as follows:

Linear case: A = 0
In the linear case (A = 0), equation (22) assumes the form Taking into account that the set of functions F n (u) is complete and orthonormal, one can rewrite equation (24) as where Thus, we can conclude that the solution of the linear eigenvalue problem can be presented in the form Equation (28) determines the frequency (ω n = λ 2 n ) of the nth eigenstate and equation (27) with N n being a normalization constant and δ kn being the Kronecker delta, yields its amplitude.
Introducing equations (26) and (27) into equation (23), we obtain that the eigenfunction , which corresponds to the eigenvalue given by equation (28), can be expressed as where θ(v) is the Heaviside step function. For even values of n, n = 2m (m = 0, 1, 2, . . .), equation (28) always has a solution. For νR → 0, the eigenvalue is given by For n = 2m + 1, m = 0, 1, 2, . . . , the bound state exists only for νR 1 and near the lower bound the energy of the bound state is given by In the limit of large radius of curvature R and moderate n, i.e., νR n + 1 2 , we obtain from equation (28) that the eigenvalues λ n are determined by Thus, the bound state energy decreases when the curvature of the chain increases and in the limit R → ∞ we obtain the straight-line result: λ = ν/2. It is interesting to return to the Cartesian coordinates and to consider the shape of the bound state wave function. Let us consider the case n = 0. It is seen from equations (B.3) and (B.9) that where the function v(x, y) is given by equation (11). When R → ∞ (straight waveguide) the wave function is localized in the y-direction only. However, for finite R, the function is localized both in x-and y-directions. The localization length in the x-direction is proportional to R. The expression of equation (33) will also be used as a starting point in our direct numerical simulations of equation (12); see section 5 below.

Nonlinear case
In the following, we restrict ourselves to the case of the waveguides with small or moderate curvature (1/R < 1). We use Darwin's expansion of the parabolic cylinder functions [32]: Using this approximation, equation (22) may be represented in the form: or taking into account that in the equivalent form, Thus, eliminating the waves in the linear medium in which the waveguide is embedded and applying the inverse Fourier transformation with respect to equation (13), leads to the following equation for the waveguide function: Thus, the dynamics of the system is described by the pseudo-differential or, in other words, the nonlocal, in time and space, equation. The nonlocal character of the nonlinear waveguide dynamics is due to the existence of two paths for the excitation energy transfer: directly along the waveguide and via the linear medium in which the waveguide is embedded.

Stationary states
We are interested here in the stationary solutions of equation (38) Here, the shape function (u) satisfies the equation Let us consider in more detail the case of repulsive excitations (A > 0). This case is the most interesting from the point of view of the interplay between the nonlinearity and curvature, because, for the straight waveguide (R → ∞), equation (40) has no localized solutions.
When the nonlinearity parameter A is sufficiently large, one can use the so-called Thomas-Fermi approximation [33] in which one neglects the term ∂ 2 u in equation (36), and one finds a density profile where u 0 = R ν 2 4λ 2 − 1. By using equations (39) and (12), it is straighforward to show that in terms of the channel wave function (20) the Hamiltonian (1) and the norm (5) may be written as Inserting equation (41) into equations (42) and (43), we obtain the following expressions for the energy and the norm of the nonlinear excitations: From equations (44) and (45)

Non-stationary dynamics
When R → ∞ equation (38) assumes the form which is the nonlinear Hilbert-NLS equation introduced in [34].
In the limit of weak nonlinearity (A 1) and small curvature (R 1) one can significantly simplify equation (38) by using the ansatz and assuming that ϕ(u, t) depends slowly on u and t. Inserting equation (47) into equation (38), we obtain Considering the scaling Thus, the behaviour of nonlinear excitations in a curved waveguide is equivalent to the behaviour of nonlinear excitations in the parabolic potential whose curvature coincides with the curvature of the waveguide. In the linear case A = 0, the corresponding eigenvalue problem gives the eigenvalues presented in equation (32). Furthermore, it is interesting that this reduction gives rise to an effective one-dimensional Gross-Pitaevskii (GP) equation, analogous to the one used in the study of Bose-Einstein condensates in cigar-shaped traps [35]. While we do not pursue this analogy in detail in the present proof-of-principle paper, we note that it would naturally be of interest to examine excitations known in the context of the GP equation, such as e.g., dark solitons [36] and their dynamics, in the present setting.

Numerical results
We start by demonstrating the results of the linear case of equations (28) and (33). The eigenvalue (energy) of the linear case as a function of curvature is shown in figure 1, while the lowest energy, bound-state wavefunction of the linear problem is given in figure 2.
In order to demonstrate that this linear bound state persists in the nonlinear limit we have performed full dynamical evolution simulations of equation (8), with an initial condition of the form expressed in equation (33), and demonstrated in figure 2. We note in passing that similar results have been obtained with other initial conditions for e.g., a Thomas-Fermi initial profile  of the form of equation (41). In particular, we show typical numerical simulation results in figure 3 for R = 10, ν = A = 1. Notice that the δ function was represented as with = 0.05. The contour plot of figure 3 shows the result after a numerical evolution of 100 time units of equation (8).
The dynamical development indicates that after an initial transient, the original linear profile slightly reshapes itself into the nonlinear solution depicted in figure 3. In the process, some radiation waves ('phonons') are shed, that are absorbed by the absorbing boundary conditions used in a layer close to the end of the domain (our computational box is of size 25 × 25).
Beyond the proof-of-principle simulations for various initial conditions theoretically derived in sections 3 and 4, we also attempted to examine the dynamics of the nonlinear excitations of the channel. This was done using the following numerical protocol: after obtaining a quasirelaxed nonlinear localized mode for the channel of the form y = x 2 /(2R), we moved the  figure 4, the process of relaxation to this new equilibrium is shown as a function of time for a very long dynamical simulation of t up to 1000 time units. In this run, we observe (after an initial transient) a slow relaxation towards the new minimum of the potential well. Notice that, despite the Hamiltonian nature of the model, the excitation of an 'internal mode' of the nonlinear wave [38] can be dissipated due to mechanisms of coupling to extended wave, phonon modes, such as the ones reported in [39].

Summary
In summary, we have shown that • the finite curvature of the waveguide provides a stabilizing effect on otherwise unstable localized states of repelling excitations; • the binding energy of both linear and nonlinear localized excitations decreases when the curvature of the waveguide increases. • such linear bound states as the ones found here persist in the nonlinear dynamical problem as localized excitations. These have been found to be robust for different initial conditions and are centred at the minimum (point of largest curvature) of the parabola. When the mode is initialized away from this minimum, it slowly relaxes to it.
While our motivation was originally provided by the embedding of a waveguide in a twodimensional photonic crystal (as well as from more general geometric considerations), interesting analogies have arisen through our investigation, which warrant further study. In particular, we note the analogy in the weakly nonlinear regime of the equation for the waveguide/quantum channel with that of the Bose-Einstein condensate behaviour in the presence of a magnetic trap in atomic physics. Studying the transmission and reflection of linear and nonlinear waves in bent waveguides and in particular the role of the localized states, is another vital issue in mass, charge and energy transport in systems with complicated geometry. Just as bent regions of chains of nonlinear oscillators [26,40] provide energy funnelling, so can the curved regions of waveguides. Another topic worthy of further investigation is a more detailed numerical study of the stability of the nonlinear localized modes we have identified. It would also be of interest to examine the transmission of linear states through the bend and their interaction with localized states, e.g., in the spirit of [41]. Finally, another direction that could be of further interest is the examination of a case of a finite (rather than infinitesimal) width channel (which can be computationally achieved, e.g., by allowing the parameter of equation (51) to vary towards larger values). The examination of thresholds for genuinely two-dimensional instabilities, such as e.g. the transverse or the snaking instability, would be of particular interest within the latter context.
Such studies are currently in progress and will be reported in future publications.
Here, the potential U f ( r ) has the form Let us consider the case where there are two interfering coherent beams of wavelength λ and frequency ω. One of them is a diverging spherical wave which is emanated by a point source situated at the point (0, y s , 0). Its electric field strength is The second beam is a plane wave which propagates along the y-axis In Let us assume now that the sources have a broad emission band which is described by the Gaussian function where k 0 is the wave vector of the carrier wave, the parameter gives the width of the emission band and the parameter A characterizes the total intensity of emission. The total field intensity at the point r is 2) and considering the case when y s ≈ y p , we obtain that the spatially dependent part of the potential U( r ) has the form which coincides with the potential which was used in equation (9).

Appendix B. Green's function
To find an expression for the Green function G(u, u ; v, v ) we expand this function in terms of the eigenfunctions of the equation