Stability and instability of nonlinear defect states in the coupled mode equations—analytical and numerical study

# Stability and instability of nonlinear defect states in the coupled mode equations—analytical and numerical study

Roy H. Goodman                 Michael I. Weinstein Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 01710Department of Applied Physics and Applied Mathematics, Columbia University, New York, NY 10027
###### Abstract

Coupled backward and forward wave amplitudes of an electromagnetic field propagating in a periodic and nonlinear medium at Bragg resonance are governed by the nonlinear coupled mode equations (NLCME). This system of PDEs, similar in structure to the Dirac equations, has gap soliton solutions that travel at any speed between and the speed of light. A recently considered strategy for spatial trapping or capture of gap optical soliton light pulses is based on the appropriate design of localized defects in the periodic structure. Localized defects in the periodic structure give rise to defect modes, which persist as nonlinear defect modes as the amplitude is increased. Soliton trapping is the transfer of incoming soliton energy to nonlinear defect modes. To serve as targets for such energy transfer, nonlinear defect modes must be stable. We therefore investigate the stability of nonlinear defect modes. Resonance among discrete localized modes and radiation modes plays a role in the mechanism for stability and instability, in a manner analogous to the nonlinear Schrödinger / Gross-Pitaevskii (NLS/GP) equation. However, the nature of instabilities and how energy is exchanged among modes is considerably more complicated than for NLS/GP due, in part, to a continuous spectrum of radiation modes which is unbounded above and below. In this paper we (a) establish the instability of branches of nonlinear defect states which, for vanishing amplitude, have a linearization with eigenvalue embedded within the continuous spectrum, (b) numerically compute, using Evans function, the linearized spectrum of nonlinear defect states of an interesting multiparameter family of defects, and (c) perform direct time-dependent numerical simulations in which we observe the exchange of energy among discrete and continuum modes.

## 1 Introduction

In optical fiber communication systems, data is carried as pulses of light. Expensive and rate-limiting steps in these systems come in processing the data at so-called optical / electrical interfaces. Processing is typically done electronically; signals are converted to electronic form, read, processed, and retransmitted. A major goal, therefore, is to bypass the optical / electrical interface and implement “all-optical” processing by using the nonlinear optical properties of the medium. Hence, there has been great interest in finding novel materials and optical structures (arrangements of materials) which effect light in different ways. Among these are Bragg gratings and Bragg gratings with defects–one dimensional arrangements, as well as higher-dimensional structures such as photonic crytal fibers in which the grating structure is transverse to the direction of propagation [43], and other two and three-dimensional structures [16].

In an optical fiber Bragg grating, the refractive index of the glass varies periodically, with period resonant with the carrier wavelength of the propagating light; see figure 1.1. Light propagation through such fibers has several interesting properties that makes them useful in technological applications. The grating structure couples forward-moving light at the resonant wavelength to backward-moving light of the same wavelength. In the low-amplitude (linear) limit, this makes the fiber opaque to light in a certain range of wavelengths, the so-called photonic band-gap. When the amplitude of the light is increased, the band-gap is shifted due to the nonlinear dependence of refractive index on intensity. Thus, a range of wavelengths that are non-propagating at low intensity are shifted into the range of propagating wavelengths (the pass-band) at higher intensities. This is the mechanism by which localized pulses known as gap solitons exist; see subsection 2.2. In the regime of weak nonlinearity and Bragg resonance, Maxwell’s equations can be reduced, via multiple scale asymptotic methods, to the nonlinear coupled mode equations (NLCME), reviewed in section 2.1 [2, 15, 19, 32].111A special case of NLCME (equation (2.1) with the self-phase modulation terms omitted) is the massive Thirring system [25], a completely integrable PDE modeling the interaction of massive fermions. The methods described in this article apply equally well to the massive Thirring system..

Gap solitons may, in theory, propagate (in the stationary reference frame) at any speed between 0 and . Here, denotes the vacuum speed of light, the refractive index of the optical fiber core and is the speed of light in the fiber without the grating. In [29], we showed via modeling, analysis and numerical simulation, that propagating optical gap solitons, could be trapped at specially-constructed defects in the grating structure. Trapping and interaction of gap solitons in a number of related structures is considered in [46, 47, 14]. Recent experimental advances have made possible the slowing of propagating gap soliton light pulses from  [23] to  [52].

The focus of the present paper is on the stability and dynamics of such trapped light. Localized defects in a grating appear as spatially localized potentials in the NLCME model. In the linear (low light intensity) limit, the resulting linear coupled mode equations with potentials have spatially localized linear defect eigenstates; see subsection 2.3. As the intensity is increased from zero, nonlinear defect modes, “pinned” at the defect location, bifurcate from the zero state at the linear eigenfrequencies; see subsection 2.4.

Trapping of a gap soliton by a defect can be understood as the resonant transfer of energy from an incoming soliton to a pinned nonlinear defect mode. In [29] we demonstrated through numerical experiments such resonant energy transfer / trapping, for sufficiently slow soliton pulses. The analogous question has also been studied for the nonlinear Schrödinger / Gross-Pitaevskii (NLS/GP) equations; see [28, 35, 36, 34] and references cited therein. For the purpose of comparison, we review the stability of and interactions between nonlinear defect modes of NLS/GP in section 3.

In order for the energy localized in a defect to remain spatially confined in a nonlinear defect mode, it is necessary that the mode be stable. Thus, in this paper we consider the stability of nonlinear defect modes for a large class of defects introduced in [29] by

• studying the linearized spectral problem about different families (branches in the global bifurcation diagram) of nonlinear defect modes; see sections 4 and 5 for analytical perturbation theory and numerics, and

• studying time-dependent simulations of the initial value problem for NLCME. We consider defects which support multiple nonlinear linear defect modes. As the time-evolution proceeds these nonlinear bound state families compete for the energy localized in the defect; see section 6.

The linearized stability of a nonlinear defect mode is governed by a spectral problem of the form:

 Σ3Hψ = β ψ,  Σ∗3=−Σ3,  H∗=H; (1.1)

see section 4. Eigenvalues are complex values of , for which (1.1) has a nontrivial solution, giving rise to a time-dependent solution of the linearized dynamics . The self-adjoint operator, , can be expressed as , where corresponds to the linear (zero intensity) coupled mode equation operator and tends to zero quadratically in the amplitude (-norm) of the nonlinear defect mode about which we have linearized.

The continuous spectrum of typically consists of two symmetric semi-infinite real intervals with an open gap centered at . A priori, since is not self-adjoint, discrete eigenvalues may lie anywhere in the complex plane, subject to constraints inherited by (1.1) from NLCME, a Hamiltonian system. In particular, if is an eigenvalue, then so are and . Thus, a necessary condition for linearized stability of the flow is that the spectrum of be real.

Now has stable spectrum and may have real discrete eigenvalues in the spectral gap or real embedded eigenvalues within the continuous spectrum. Numerous different instability-onset scenarios corresponding to different types of bifurcations arise as deforms to , as the intensity ( norm) of a nonlinear defect mode is increased along the bifurcation branch. We focus on one such scenario here and in more detail in section 4.2. The other scenarios are discussed briefly in section 4.3.

The first scenario occurs when has a symmetric pair of real embedded eigenvalues within the continuous spectrum. Generically these perturb, for positive and arbitrary, to two pairs of complex conjugate eigenvalues and therefore yield instability. The mechanism is related to the notion of Krein signature. A detailed perturbation calculation demonstrating this instability is presented in section 4.2, showing that the instability is of order , or equivalently, of order equal to the fourth power of the defect mode amplitude. In the particular examples we study numerically, these instability rates are observed, but are seen in some cases to be quite small.

In section 6 we explore the temporal dynamics of the initial value problem for NLCME via direct numerical simulation. We present simulations with various defects supporting one, two and three modes which illustrate both the spectral instability scenarios explored in section 4 and the temporal dynamics of energy exchange among defect modes and radiation modes. The mechanisms are similar to, but more complicated than, those studied by Soffer and Weinstein [67, 68] for the nonlinear Schrödinger / Gross-Pitaevskii (NLS / GP) equation, who consider NLS / GP with a “defect potential” supporting two bound states and thus, two branches of nonlinear defect modes (ground and excited), which compete for energy confined by the potential; see also [71]. In the NLS/GP problem, resonance of an embedded eigenvalue, associated with the linear excited bound state, with the continuous spectrum leads to energy transfer from the excited state to the ground state and to radiation modes. The asymptotic distribution of excited state energy is also studied; see [67, 68, 74].

We recap the structure of the paper. In section 2, we introduce NLCME, its gap solitons, and two families of defects and their linear and nonlinear bound states. We also summarize the results of a previous paper in which these defects are used to trap gap solitons. In section 3 we review some analogous results for the NLS equation with a localized defect. Section 4 discusses the properties of the linearized operator, outlining several types of instability with particular attention to the bifurcations of embedded frequencies. Section 5 provides a brief introduction to the Evans function, a useful analytical tool for exploring the discrete spectrum of the linearized operator, as well numerical results using the Evans function to study the spectral stability of nonlinear defect modes. Section 6 shows time-dependent simulations confirming the stability predictions of the previous section. Section 7 contains a short summary and discussion. The appendices contain a list of symbols, a derivation of NLCME with potentials from the problem of wave propagation in a Bragg grating with defects, and a detailed description the numerical measures that are taken to ensure accurate computations of both the nonlinear defect modes and the discrete spectrum of the linearization.

## 2 NLCME and NLCME with defects

### 2.1 The Nonlinear Coupled Mode Equations

Consider an electromagnetic wave-packet, a pulse which is spectrally concentrated about a carrier wavelength, propagating in a waveguide whose refractive index varies periodically in the direction of propagation. If carrier wavelength is in resonance with the waveguide periodicity (Bragg resonance), then backward and forward waves couple strongly.

Envelope equations for these forward and backward wave amplitudes satisfy the nonlinear coupled mode equations. In appendix B we show that localized defects in the periodic structure give rise to the following modification of the NLCME, which we call by the same name, as derived in [29, 73]:

 i∂TE++i∂ZE++κ(Z)E−+V(Z)E++(|E+|2+2|E−|2)E+=0i∂TE−−i∂ZE−+κ(Z)E++V(Z)E−+(|E−|2+2|E+|2)E−=0. (2.1)

The coefficient functions (“potentials”) and are determined from the refractive index profile. Here, is the “slow” coordinate along the direction of propagation, is a slow time variable and the full electric field is given by

 E=E+(Z,T)ei(z−t)+E−(Z,T)e−i(z+t)+c.c.,

i.e. and are complex envelopes of rapidly varying electromagnetic fields, assumed to be small amplitude and linearly polarized. The length and time variables are chosen such that the wavenumber and frequency are both one, i.e. if in dimensional variables , the carrier wave has wavelength and period , then we choose and . The condition that the grating is uniform (exactly periodic) away from the defect region implies that , constant, and as . Defining , equation (2.1) may be rewritten as

 (i∂T+iσ3∂Z+V(Z)+κ(Z)σ1)→E+N(→E,→E∗)→E=0, (2.2)

Where and are the Pauli matrices and , the superscript asterisk represents complex conjugation, and represents the nonlinear term of (2.1),

 N(→E,→E∗)=(|E+|2+2|E−|200|E−|2+2|E+|2). (2.3)

A detailed derivation, including careful accounting of the nondimensionalization and scalings, is given in [29, 32, 73].

The NLCME system (2.1) has two conserved quantity, the total intensity and the Hamiltonian (energy) . The total intensity, , is defined as

 I(→E)=∫∞−∞(|E+|2+|E−|2)dZ, (2.4)

i.e. the square of the norm of the solution. The Hamiltonian is given by

 H(→E)=∫∞−∞(iE+∂ZE+−iE−∂ZE−+κ(Z)(E+E∗−+E∗+E−)+V(Z)(|E+|2+|E−|2)+2|E+|2|E−|2+12|E+|4+12|E−|4) dZ (2.5)

In the absence of a defect ( constant, , NLCME also conserves momentum, as a consequence of Noether’s theorem.

The primary focus of this paper is the standing wave solutions of (2.2) of the form

 →E(Z,T)=→E(Z)e−iωT,

for which (2.2) reduces to the nonlinear eigenvalue problem

 (ω+ iσ3∂Z+V(Z)+κ(Z)σ1)→E+N(→E,→E∗)→E=0,→E∈L2×L2. (2.6)

.

### 2.2 Constant coefficient NLCME—The spectral gap and “gap solitons”

NLCME with a uniform grating has , constant, and , and has been studied extensively. In the linearized system, setting the cubic terms to zero, one may look for solutions of the form

 (E+E−)(Z,T)=ei(kZ−ωT)(e+e−)

and find the dispersion relation , so that for frequencies , is purely imaginary, i.e. there exists a spectral gap . This results in the reflection of light at the resonant frequencies in the gap.

The fully nonlinear coupled mode system was shown in [2, 15] to support uniformly propagating bound states called gap solitons, parameterized by a velocity and a detuning parameter satisfying and :

 E±=±Δ∓1α√κ∞2 (sinδ) ei(η+σ)sech(θ∓iδ/2); (2.7)

where:

 Δ=(1−v1+v)14;α=√2(1−v2)3−v2;eiη=(−e2θ+e−iδe2θ+eiδ)2v3−v2; θ=κ∞√1−v2(sinδ)(z−vt);σ=κ∞√1−v2(cosδ)(vz−t).

In the Lorenz-shifted reference frame, the term is responsible for an internal oscillation with frequency which is inside the spectral gap, the nonlinearity having served to effectively shift the gap. This is also known as self-induced transparency. A comprehensive introduction to gap solitons is given in the review paper by de Sterke and Sipe [19]. Gap solitons may propagate in the laboratory reference frame at any velocity below the speed of light (here by the leading to B.6) and thus are intriguing as possible components of all-optical communications systems. Zero-speed waves, for example, are candidates for bits in an optical buffer or memory device.

### 2.3 Defect modes of the linear coupled mode equations

In the zero-intensity limit, we ignore the nonlinear term in (2.6), giving a linear eigenvalue problem

 ( ω∗I + hV,κ ) →E∗ ≡ (ω∗I+ iσ3∂Z+V(Z)+κ(Z)σ1)→E∗=0, →E∗∈L2×L2. (2.8)

Hereafter, frequencies and vectors with subscript asterisks will refer to solutions the linear eigenproblem (2.3), and those without will refer to solutions of (2.6).

Note that is self-adjoint on and therefore has real spectrum. Moreover, we assume and rapidly as , and therefore this spectral problem has two branches of continuous spectrum and a finite number of discrete frequencies given by

 spec(ω∗I+hV,κ) = speccontinuous∪specdiscrete={ω∗∈R:|ω∗|≥κ∞}∪{ωj∗:j=1,…,N},

where the discrete frequencies satisfy .

We consider two families of defects in the current study which we refer to as “even defects” and “odd defects”.

#### (a) Even defects

If we specify to be even and , then it is simple to show that if is a real eigenvalue of (2.3), then so is . A simple argument eliminates the possiblity of nontrivial null eigenstates. Thus, system (2.3) must have an even number of discrete eigenvalues with corresponding eigenfunctions. As we do not know of any for which the eigenvalue problem (2.3) may be solved in closed form, we choose a relatively simple form of ,

 κ(Z)=1−bsech(kZ). (2.9)

Thus, and the linear operator has continuous spectrum . We have found via numerical simulation that as is increased with held fixed, the existing frequencies migrate toward the origin, and new discrete frequencies bifurcate in pairs from opposite endpoints of the band gap at . Such emergence from the edges has been studied for related problems, for example, in [42, 58].

#### (b) Odd defects

In [29], we construct a three-parameter family 222By a rescaling of the space and time coordinates, we may set , so in reality this is a two-parameter family, but the above form makes it easier to create different defects with the same limiting grating strength . of defects of the form

 κ(Z)=√ω2+n2k2tanh2(kZ);V(Z)=ωnk2sech2(kZ)2(ω2+n2k2tanh2(kZ)). (2.10)

with , , and nonzero. 333In the special case , then and . The expressions for the eigenmodes are also altered. Standing wave solutions exist of the form:

 →E(0)∗=(eiΘ−se−iΘ)sechn(kZ)e−iω0∗T,  ω0∗=ω, s=signω (2.11) Θ=∫Z0V(ζ) dζ=12arctannktanh(kZ)ω. (2.12)

For , the defect supports a total of distinct eigenvalues, where is the greatest integer less than or equal to . Its eigenvalues are

 ω±j∗=±√ω2+(2nj−j2)k2;j=1,…,n−1. (2.13)

If and are held constant, while is increased, the “ground state” eigenvalue remains constant, while the other eigenvalues move from the edge of the spectral gap toward . Each time passes through an integer value, two new eigenvalues are created in edge bifurcations from the two ends of the continuous spectrum.

Reference [31] contains the general formula for the eigenfunctions. The solutions can ultimately be expressed in terms of hypergeometric functions of , which simplify to Legendre functions of when , and can be expressed as algebraic combinations of and . The general formula is quite complicated.

In much of what follows we specialize to the case , in which case there are three linear defect modes. The first eigenfunction, with frequency , is given by (2.11). The remaining two, with frequencies , are given (in a non-normalized form) by

 →E(±1)∗=((k+i(ω±1∗+ω)tanhkZ)eiΘ(Z)s(k−i(ω±1∗+ω)tanhkZ)e−iΘ(Z))sechkZ. (2.14)

### 2.4 Defect modes of nonlinear coupled mode equation

Consider equation (2.6), the nonlinear eigenvalue problem solved by the nonlnear defect modes. In the small amplitude (linear) limit, we expect solutions to be well-approximated by those of the linear eigenvalue problem (2.3), whose explicit eigenstates are displayed in section (2.3).

In [29], we used perturbation analysis to show that in analogy with the NLS/GP case [60], there exist nonlinear defect mode solutions of (2.6)

 →E(j)(Z) = α( →E(j)∗(Z)+O(|α|2) ) and (2.15) ω =ωj∗+q|α|2+O(|α|4),α∈C,α→0 (2.16) q =−⟨Ej∗,N(Ej∗,E∗j∗)Ej∗⟩⟨Ej∗,Ej∗⟩ (2.17)

Notice that , due to the focusing character of the nonlinearity. Thus, as the amplitude is increased, the nonlinear frequency is shifted toward the left edge of the spectral gap.

We aim to compute the nonlinear modes and their associated frequencies as accurately as possible. This is important so that errors in the numerical solutions contribute as little as possible to errors in their numerically calculated linearized spectrum and determination of stability. Details of the numerical calculation are provided in appendix C. Briefly, the derivatives are discretized using Fourier transforms, while the infinite interval is truncated using a nonuniform-grid method and the resulting algebraic equations derived are solved using MINPACK routines [53].

Figure 2.1 shows a typical example of the dependence of the defect mode’s frequency on its amplitude. The frequency of each of the three defect modes shifts to the left as the intensity is increased; . At small intensities the frequency depends linearly on the intensity, and at larger intensities, the curves steepen near the left band edge. We believe that these branches terminate at the left band edge, but the nonlinear modes decay very slowly as the edge is approached, and thus become difficult to calculate numerically.

### 2.5 Summary of previous numerical experiments

In [29], we explored the behavior of gap solitons incident on localized defects of the form (2.10). We found parameter regimes in which the soliton—or at least the energy it carries—can be captured. We hypothesized the existence of a nonlinear resonance mechanism between the soliton, and a nonlinear defect mode. By (2.7), the gap soliton oscillates with an internal frequency of about  as it propagates. If there exists a nonlinear defect mode of the same frequency and smaller -norm, then the solitary wave may transfer its energy to the defect mode. In particular, as is increased from to this internal frequency decreases from to . This is the thickest curve in figure 2.1; a family of states of the translation invariant NLCME, which bifurcates from the zero state at the band edge frequency.

For sufficiently small (frequency near the right band edge, marked (a) in the figure), there exists no nonlinear defect mode of the same frequency and lower amplitude to which the gap soliton can transfer its energy. In this case, the defect behaves as a barrier. For solitary waves above a critical velocity, the pulse passes by the defect with almost all its amplitude and its original velocity almost unchanged. Below the critical velocity, the solitary wave is reflected, again almost elastically. If, by contrast, there exists a nonlinear defect mode resonant and of with the solitary wave and of smaller amplitude, then the pulse may be trapped, for example points (b) and (d). In this case there exists a critical velocity, below which solitary waves are transmitted (inelastically) and below which they are largely trapped. Behavior at the point marked (c) is similar to that at (a)–the defect mode that bifurcates from is a larger-amplitude state than the gap soliton, so it cannot be excited, and the defect mode bifurcation from is not resonant with this frequency.

Similar behavior, namely the dichotomy between the nonresonant and elastic transmit/reflect behavior and the resonant and inelastic capture/transmit behavior, has also been seen for nonlinear Schrödinger solitons [28, 44]. Estimates for in several related problems and an explanation for the existence of a critical velocity are given in [27] and references therein. This trapping phenomenon was subsequently seen by Dohnal and Aceves for a two-dimensional generalization of equation (2.1[1, 21].

Trapping could be more accurately described as a transfer of energy from the traveling soliton mode to a stationary nonlinear defect mode. This is pictured in figure 2.2 reprinted from [29]. This figure describes the evolution of a pulse captured by a wide defect supporting 5 linear defect modes. Energy initially moves into the mode associated with the frequency but is subsequently transferred to the mode associated with frequency . A part of the energy is transferred to radiation modes and is dispersed to infinity. It is thus important for us to understand the stability of these nonlinear modes, to assess their suitability as long-time containers for electromagnetic energy, as well as the mechanisms through which energy moves among discrete and continuum modes.

## 3 Dynamics and energy transfer for NLS/GP

Many of the phenomena and mechanisms present in the dynamics of NLCME (2.1) are present in the related and more studied nonlinear Schrödinger / Gross-Pitaevskii (NLS/GP) equation:

 i∂tΦ=−ΔΦ+V(x)Φ+g|Φ|2Φ, (3.1)

where is complex-valued, , . corresponds to a repulsive or defocusing nonlinear potential, while corresponds to an attractive or focusing nonlinear potential. Before embarking on a study of NLCME, we briefly review results for NLS/GP.

To fix ideas, we take to be a smooth potential well (), decaying to zero rapidly as . NLS/GP is a Hamiltonian system with conserved Hamiltonian energy

 H[U]=∫Rn(|∇U|2+V(x)|U|2+g2|U|4)dx. (3.2)

Solitary standing wave solutions are solutions of the form: , from which we obtain the nonlinear eigenvalue problem for

 ΩΨ=−ΔΨ+V(x)Ψ+g|Ψ|2Ψ;Ψ∈L2 (3.3)

In the low intensity (linear) limit, (3.3) reduces to the the linear eigenvalue problem for the Schrödinger operator .

The operator has continuous spectrum covering the non-negative half-line: and, in general, a finite number of negative discrete eigenvalues , with corresponding eigenfunctions ; recall subscript asterisks denote solutions in the linear limit. The eigenfunction corresponding to is the ground state, a minimizer of the (linearized) energy , subject to the constraint: . Nonlinear defect mode families, , bifurcate from the zero state at the linear eigenfrequencies [60]. The nonlinear ground state, can also be characterized variationally as a minimum of subject to the constraint . For small -norm we have for and

 Ψαj(x)=α( ψj∗(x) + O(|g| |α|2) ) Ωj(α)=Ωj∗ + gc2j∗|α|2 + O(g2|αj|4),  αj∈C

where is a positive constant.

Soffer and Weinstein [67, 68] studied the dynamics of the initial value problem for NLS/GP (3.1) in the case where the potential well, , supports exactly two eigenstates, “linear defect modes”, and . By the above discussion, there exist two bifurcating branches of nonlinear bound states . These are used to parameterize general small amplitude solutions of NLS/GP:

 Φ(x,t) ∼ Ψα0(t)(x) + Ψα1(t)(x) + Φdispersive(x,t)

and it is shown that generic finite energy solutions of the initial value problem converge as to a nonlinear ground state: , locally in . This ground state selection phenomenon has been experimentally observed; see  [49].

The very large time dynamics, which involve energy leaving the excited state and being transferred to the ground state and radiation modes is governed by nonlinear master equations

 dP0(t)dt ∼ΓP21(t)P0(t) dP1(t)dt ∼−2ΓP21(t)P0(t). (3.4)

Here, and is non-negative and generically strictly positive, provided an arithmetic condition implying coupling to of the discrete to continuum radiation modes at second order in :

 2Ω1∗−Ω0∗>0. (3.5)

The expression for is a nonlinear analogue of Fermi’s golden rule (see also [66]), which arises in the calculation of the spontaneous emission rate of an atom to its ground state [17].

The resonance condition (3.5) also appears in a natural way in the linearized (spectral) stability analysis of nonlinear defect modes. In particular, consider the linearization about a nonlinear excited state. This yields a linear spectral problem of the form: , where is a two by two self-adjoint matrix operator and is the standard Pauli matrix. Now , where tends to zero quadratically in the nonlinear excited state amplitude. Since is spatially localized, the continuous spectrum of is given by two semi-infinite intervals, the complement of an open symmetric interval about the origin. Now under the resonance condition (3.5), has an embedded eigenvalue within the continuous spectrum. The perturbation theory of embedded eigenvalues is a fundamental problem in mathematical physics; see, for example, [59, 65]. These embedded eigenvalues can be shown generically to perturb for arbitrarily small to complex eigenvalues, corresponding to instabilities [18].

In subsequent sections we explore the analogous picture for NLCME. In particular, in section 4, we obtain an analogous (more complicated) resonance condition, yielding embedded eigenvalues for a spectral problem, perturbing to instabilities and analogous dynamics of energy transfer among modes.

We conclude this section by considering numerical simulations for NLS/GP, with a simple family of potentials that support two linear defect modes is

 VL(x)=−2sech2(x−L)−2sech2(x+L).

More detailed numerical simulations and their interpretation in light of [67, 68] is given in [62].

A potential well has a single localized eigenfunction with frequency . For all values of , the potential supports exactly two stationary solutions. For sufficiently large, the normalized ground state is given by , with frequency where is positive and exponentially small in . The excited state is , with frequency . As is decreased, increases and increases toward . For , condition (3.5) is satisfied.

We have simulated the initial-value problem for this problem with initial conditions given by a linear combination of the two eigenfunctions, in both the resonant and non-resonant cases. The results of these simulations is shown in figure 3.1 and the behavior in the two cases is strikingly different. Subfigure (a) shows the case =1, in which the ground state grows slightly, while the excited state decays. In addition to this transfer of energy, there is a fast oscillation in the amplitudes of the two modes. Subfigure (b) shows the nonresonant case . Here there is a much larger amplitude oscillation between the two modes, but none of the one-way energy transfer. In section 6, we perform analogous numerical experiments, which confirm analytical / numerical predictions of section 5.

## 4 Linearized stability analysis of nonlinear defect modes

We now study the linearization of the variable-coefficient NLCME (2.1) about a nonlinear defect mode. In the introduction, we reviewed one type of instability that may arise. The analysis of subsection 4.2 is focused on this scenario—instabilities arising from perturbations of eigenvalues, embedded within the continuous spectrum. We first discuss the conditions under which, in the limit of vanishing-amplitude defect modes, the linearization contains discrete eigenvalues embedded in the continuous spectrum, and then develop a time-dependent perturbation theory that shows when these embedded modes may lead to instability. In section 4.3, we discuss other types of instability that may arise.

### 4.1 Conditions for embedded eigenvalues in the linearization

Letting be a solution to the nonlinear eigenvalue problem (2.6) constrained to have intensity , we linearize about by letting

 E+ =(E++y1(Z,T))e−iωT E∗+ =(E∗++y3(Z,T))eiωT E− =(E−+y2(Z,T))e−iωT E∗− =(E∗−+y4(Z,T))eiωT.

Letting yields the eigenvalue problem for and an function

 β→y(Z)=−Σ3⋅(V(Z)+ω+i(σ300−σ3)∂Z+(σ100σ1)κ(Z)+W(Z))→y (4.1)

where is a Pauli matrix. The matrix multiplication operator is given by

 W=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝2|→E|22E+E∗−E2+2E+E−2E∗+E−2|→E|22E+E−E2−E∗+22E∗+E∗−2|→E|22E∗+E−2E∗+E∗−E∗−22E+E∗−2|→E|2.⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ (4.2)

An eigenfunction is a square-integrable solution of (4.1), with corresponding discrete eigenvalue . Thus if any non-trivial solutions exist with , the nonlinear defect mode is unstable.

As , the spectrum of the linearization about a given defect mode approaches

 speclinearized={±(ωj∗−ω∗):ω∗∈speclinear}

In particular, from the form of equation (4.1) in the limit of , the branches of continuous spectrum must be given by

 spec1={β:|β|≥κ∞−|ωj∗|}

and

 spec2={β:|β|≥κ∞+|ωj∗|}

so that spectrum has multiplicity two for . In addition the discrete spectrum is given by

 specdiscrete={β±j,k=±(ωj∗−ωk∗):1≤k≤N}.

This implies that the spectral gap is the interval

 gap=(−κ∞+|ωj∗|,κ∞−|ωj∗|) (4.3)

This is summarized in figure 4.1. If for some ,

 |β±j,k|≥|κ∞|−|ωj∗|. (4.4)

then the discrete eigenvalue is embedded in the continuous spectrum. We will refer to the boundary between the band gap and the multiplicity-one continuous spectrum as the “primary band edge” and the boundary between multiplicity-one and multiplicity-two continuous spectrum as the “secondary band edge.”

Example 1 For “even” defects of the form (2.9) with exactly two discrete eigenvalues , therefore there will exist an embedded frequency if , i.e. if . Since in this example, the condition for an embedded eigenvalue is .

A similar calculation shows that for “odd” defects of the form (2.10), the linearization about the “excited states,” always produces embedded frequencies, while the linearization about the “ground state” has an embedded frequency if . Under perturbation, the frequencies in the gap will generically remain real at small amplitude, while the embedded frequencies will be shown to move off into the complex plane, giving growing modes.

### 4.2 Perturbation theory of embedded eigenvalues—time dependent approach

In example 1, we observe that there are cases when the linear spectral problem (4.1), in the limit of vanishing nonlinear defect mode amplitude, has eigenvalues embedded in the continuous spectrum. The perturbation theory of embedded eigenvalues in continuous spectrum is an important fundamental question in mathematical physics [59, 65]. In the self-adjoint case, such embedded eigenvalues perturb to resonances, time-decaying states; self-adjointness precludes instability. In the present non-self-adjoint setting, instability cannot be precluded.

We sketch a time-dependent approach to the theory of embedded eigenvalues [65, 67, 68] to study instabilities which arise. A time-independent approach can also be applied; see [18]. In particular, we sketch a proof of the following

###### Proposition 1.

Consider the linear spectral problem (4.1), associated with a branch of nonlinear bound states , which bifurcates from the zero solution at frequency . If the linearized operator corresponding to the zero amplitude limit (the operator on the right hand side of (4.1) with ) has an embedded eigenvalue in its continuous spectrum (see condition (4.4) ), then for arbitrary sufficiently small nonlinear bound states of order , the linearized evolution equation has an exponentially growing solution, with exponential rate , where (generically ) is of order and is given by the expression (4.13).

Recall that the nonlinear coupled mode equations are a system of semilinear evolution equations for complex wave amplitudes and . Using (2.12), the changes of variables and simplify system (2.2) to

 (i∂T+iσ3∂Z+(0~κ(Z)~κ∗(Z)0))(~E+~E−)+N(~E+,~E−,~E∗+,~E∗−)(~E+~E−)=0, (4.5)

The linearized evolution equation, governing the perturbation:

 ϕ = →y = (y+,y−,y∗+,y∗−)e−iωT (4.6)

about a fixed nonlinear defect mode has the form

 i∂tϕ = Σ3 H ϕ (4.7)

Here, , with and self-adjoint operators.

 H0 = (h000h∗0), (4.8)

where the matrix is naturally expressed in terms of blocks, where

 h0 = (ωj−i∂Z−~κ(Z)−~κ∗(Z)ωj+i∂Z) (4.9)

and the matrix is given above in (4.2).

Note that if solves the spectral problem , then solves the adjoint spectral problem, . If we introduce, , the normalized adjoint eigenvector:

 ~ψ = Σ3 ψ(Σ3ψ,ψ),   such that (~ψ,ψ) = 1,

where denotes the inner product on -valued vector functions.

We assume that there is a spectral decomposition associated with the operator . If , then can be decomposed in terms of its bound state and continuous spectral parts:

 f = Pb f + Pc f,

where

 Pb f = ∑k ( ~ψk,f ) ψk Pc f = ∫σc(Σ3H) ( ~ψλ,f ) ψλ dλ, (4.10)

where

 ~ψk = Σ3 ψk( Σ3ψk , ψk ), such that( ~ψk , ψk ) = 1

An analogous spectral decomposition holds for the operator .

Consider the situation, in which has a simple eigenvalue embedded in the continuous spectrum and corresponding eigenvector, .

What are the dynamics for the perturbed evolution equation ?

Following the analysis of  [65], we consider the initial value problem with data given by the unperturbed state, :

 Hψ0 = λψ0,  λ0∈σc(Σ3H0) (4.11)

For small (small amplitude nonlinear defect states), we seek a solution of the perturbed evolution equation

 i∂tϕ = Σ3( H0 + W ) ϕ. (4.12)

We show the existence of an exponential instability with growth proportional to where

 Γ=−π|(ψ~λ0,Wψ0)|2(Σ3ψ0,ψ0). (4.13)

Here is a generalized eigenfunction corresponding to the shifted frequency  (4.16). If is an embedded eigenvalue of with corresponding eigenfunction having negative “Krein signature”, , then .

We seek

 ϕ(t) = a(t)ψ0 + ϕ1(t),   ( ~ψ0,ϕ1(t) ) = 0 (4.14)

Let denote the spectral projection (with respect to the unperturbed operator ) onto the part of the spectrum which is (i) continuous, (ii) bounded away from infinity and (iii) from thresholds (endpoints of branches of continuous spectra). In [65] it is shown that the full dynamics are subordinate to the coupled dynamics of and , which are approximately governed by the system

 i∂ta(t) = ~λ0a(t) + ( ψ0,Wϕd(t) )(Σ3ψ0,ψ0) i∂tϕd(t) = Σ3H0ϕd(t) + a(t)P#cΣ3Wψ0, (4.15)

where

 ~λ0 = λ0 + (ψ0,Wψ0)(Σ3ψ0,ψ0 ). (4.16)

Terms neglected in arriving at (4.15) can be treated perturbatively  [65]. We consider the system (4.15) with initial data

 a(0)=1,   ϕd(0)=0

corresponding to the unperturbed embedded state.

We determine the asymptotic (in time) dynamics of (4.15) by solving for as a functional of and then substituting the result into the equation for . To facilitate this, we first extract the rapidly varying time dependence of by setting

 A(t) = ei~λ0ta(t) (4.17)

Then, the system (4.15) can be equivalently rewritten as

 ∂tA(t) = −i ei~λ0t ( ψ0,Wϕd(t) )(Σ3ψ0,ψ0 ) (4.18) i∂tϕd(t) = Σ3H0ϕd(t) + e−i~λ0tA(t) P#cΣ3Wψ0, (4.19)

Now solving for with the initial condition , gives by duHamel’s principle

 ϕd(t) = −ie−i~λ0t ∫t0 e−i(Σ3H0−~λ0I)(t−s) A(s) P#cΣ3Wψ0 ds (4.20)

Substitution of (4.20) into (4.18) yields the closed nonlocal equation for :

 (Σ3ψ0,ψ0) ∂tA(t) = −( Wψ0 , ∫t0 e−i(Σ3H0−~λ0I)(t−s) A(s) P#cΣ3Wψ0 ds ) (4.21)

Denote by