Discrete vortex switching in nonlinear optical waveguides

Introduction to the discrete nonlinear Schrödinger equation model

In this work we are studying optical waveguides, we begin by considering the electric field of a monochromatic wave propagating in a waveguide along the z-axis, having some direction of polar-ization eï :

(1) E = EHx, yL CosHk_zz - wtL eï.

Fig.1 A single square waveguide with the variation of the strength of the electric field illustrated in the direction of propagation and the transverse direction at a fixed time and height.

From Maxwell’s equations it follows that to determine E(x,y), how the field varies in the plane perpendicular to the direction of propagation, one has to solve Helmholtz equation subject to the boundary condition of the waveguide [1]. If we think of a square shaped waveguide with walls made of an ideal conductor, where there can be no electric field, the solution will be sines and cosines vanishing at the boundary as in Fig.1. This case is also illustrated in Fig.2 by the green solution.

Fig.2 Three adjacent waveguides where the variation of the electric field strength in a direction perpendicular to the direction of propagation is illustrated. The blue shaded regions represent waveguides and an electric field has been launched in the center waveguide. The green solution corresponds to a situation where the unshaded region is an ideal conductor and the red overlapping

solution to the waveguides being enclosed by a dielectric.

In the case of a finite dielectric the red solution in Fig.2 shows an overlap between adjacent waveg-uides. In this work we’re not interested in the variation of the field in the transverse direction in a specific waveguide, instead we take an average of the waveguide and call that a site for which we try to model the field as it envelops in the direction of propagation z. Because of the overlap illus-trated in Fig.2 there must be a coupling between adjacent sites when considering a system of adjacent (coupled) waveguides.

The model used in this work for studying systems of coupled nonlinear waveguides is the discrete nonlinear Schrödinger equation (DNLS) [2],[3]:

(2) ä ¶zym= -ym+1- ym-1- ym¤

2 ym

Where ¶z denotes derivation with respect to the direction of propagation and ymHzL is a complex

representation of the electric field in waveguide m. Note that there is no time dependence in y_m and that the DNLS models the envelope of the electric field of a single mode. The y_m±1 terms are a linear nearest neighbor coupling resulting form the overlap previously discussed. The nonlinear part models the interaction between the electric field and a nonlinear dielectric function of the material in the waveguide. The types of systems considered here are ones with periodic boundary conditions, see Fig.3, where y_N+1= y1 for a system with N number of waveguides.

Fig.3 Illustration of a six site system with periodic boundary condition, each black line representing a waveguide and the evolution of ym(z) represented by the oscillating curves. The red arrows

illus-trate a power flow between the sites, forming a vortex solution.

Because of the coupling between neighboring sites and the periodic boundary condition it is possi-ble to obtain solutions where a power flow (poynting power [4]) circulates from site to site over the whole system, a discrete vortex. The direction of this flow (clockwise or counter clockwise) is related to the topological charge (TC), formally defined as the sum of the phase differences between adjacent sites in units of 2p. Depending on your sign convention a solution with TC=+1 would have a power flow in the opposite direction from a TC=-1 solution. Because of this discrete-ness this is an interesting system for information processing, communicating via the TC [2][3]. This work investigates how the power flow, and therefore the TC, behaves and how one could flip the TC in a controlled fashion in the six site system.

Inference from Maxwell’s equations

The DNLS models the propagation of a single mode of the electric field in a waveguide with a nonlinear dielectric function. To see this, we start with Maxwell’s equations without sources in SI units,

(3) Ñ × D = 0 Ñ × B = 0

Ñ ´ E = -¶tB Ñ ´ H = ¶tD

Where the displacement field D = eE (= e₀E+ P), P is the charge polarization and the magnetic flux H= B/m. For the moment, the dielectric function e is assumed to be linear.Taking the curl of Fara-day’s law above and expressing Ñ ´ B in terms of the displacement field via the Ampéré-Maxwell law one obtains the following wave-equation:

Where the displacement field D = eE (= e0E+ P), P is the charge polarization and the magnetic flux

H= B/m. For the moment, the dielectric function e is assumed to be linear.Taking the curl of Fara-day’s law above and expressing Ñ ´ B in terms of the displacement field via the Ampéré-Maxwell law one obtains the following wave-equation:

(4) IÑ2

-me ¶_t2M E = 0

Having wave propagation in a waveguide in mind, we choose z to be the direction of propagation and look for a solution of the form

(5) E = 1 2 IE¦Hx, yL×e iHkzz-wtL_{+ E} ¦ * Hx, yL×e-iHkzz-wtLM

Where E¦(x,y) is the variation of the electric field in the plane perpendicular to the z-axis, having

E_¦ constant recovers an ordinary plane wave solution.

Writing the displacement field in terms of the charge polarization D = e0E+P and expanding P in a

Taylor series to take into account nonlinear effects, assuming that the charge polarization is parallel to the electric field, results in

(6) IÑ2_-me 0¶t 2_{M E = m ¶} t 2 P = m ¶_t2Je0cH1LE + e0cH3LE 3 N

Where cHiL is the electric susceptibility of order i, and the expansion has been terminated at order 3. The material used in the waveguide has a Kerr type nonlinearity [2], therefore there are no even terms in the expansion. This is because such materials have inversion symmetry, meaning that Eq.6 must be invariant under E® -E, resulting in all even cHiL=0. The usual notation is to denote the dielectric function e = e0I1 + cH1LM, taking care of the linear part as in Eq.4. Using this Eq.6 can be

rewritten as

(7) IÑ2

-me ¶_t2M E = me₀cH3L¶_t2E3.

This looks familiar except for the perturbation in the right hand side of Eq.7. We now assume that this change will perturb the solution Eq.5 in the following way:

(8) E =

1

2 IyHrL×e

iHkzz-wtL_{+ y}- HrL×e-iHkzz-wtLM

Where yHr) now depends on the direction of propagation z. Assuming a specific direction of polariza-tion Eq.7 reduces to a scalar equapolariza-tion, and upon substituting the perturbed solupolariza-tion Eq.8 the follow-ing equation for the evolution of yHr) is deduced:

(9) IÑ2 -me ¶_t2M E - me₀cH3L¶_t2E3= I¶z 2 +Ñ_Þ2+mew2M 1 2 Iye iHkzz-wtL_{+ y}- e-iHkzz-wtLM -me₀ cH3L¶_t2 1 8I3 y¤ 2 yeiHkzz-wtL_{+ 3} y¤2_y- e-iHkzz-wtL_{+ y}3_e3 iHkzz-wtL_{+ y}-3_e-3 iHkzz-wtLM = 0

The generation of third harmonics is assumed to be negligible, thus the two last terms are thrown away. Focusing on y and y- separately now and operating with ¶_z2 assuming a slowly varying field in z such that ¶_z2y<< kz¶zy results in, for y:

(10) 1 2I2 äkz ¶zy - kz 2 yM eiHkzz-wtL₊IÑ Þ 2 +mew2M 1 2 yeiHkzz-wtL_{- me} 0cH3L¶t 23 8 y¤ 2 yeiHkzz-wtL_{= 0}

Noting that k_z2= mew2 from the original solution and canceling the factors of eiHkzz-wtL_{, we get}

(11) 2 äk_z¶zy + ÑÞ 2 y + me₀w2cH3L 3 4 y¤ 2 y = 0

for the evolution of y in a single waveguide.

To turn this continuous nonlinear Schrödinger equation for a single waveguide to a discrete DNLS modelling a discrete number of coupled sites, the Ñ_Þ2 is replaced by a discrete difference operator between distinct sites (waveguides) in the xy-plane. By noting that the second derivative of a func-tion can be written as 8fHx + dL + fHx - dL - 2 fHxL<•d2 as d®0 it is natural to then replace Ñ_Þ2 by yIr - dM + yIr + dM - 2 yHrL º y-1 + y1 -2 y0, where d is a constant vector between the sites. In

other words a coupling operator between the zero:th waveguide and its nearest neighbors. Inserting the difference operator in Eq.11 we now write down the equation for the evolution of ym(z), the

for the evolution of y in a single waveguide.

To turn this continuous nonlinear Schrödinger equation for a single waveguide to a discrete DNLS modelling a discrete number of coupled sites, the Ñ_Þ2 is replaced by a discrete difference operator between distinct sites (waveguides) in the xy-plane. By noting that the second derivative of a func-tion can be written as 8fHx + dL + fHx - dL - 2 fHxL<•d2 as d®0 it is natural to then replace Ñ_Þ2 by yIr - dM + yIr + dM - 2 yHrL º y-1 + y1 -2 y0, where d is a constant vector between the sites. In

other words a coupling operator between the zero:th waveguide and its nearest neighbors. Inserting the difference operator in Eq.11 we now write down the equation for the evolution of y_m(z), the excitation in waveguide m: (12) 2 äkz¶zym+ ym+1+ ym-1- 2 ym+ 3 4 me₀w2cH3L ym¤ 2 ym= 0

In order to turn Eq.12 into the DNLS Eq.2 one simply makes the substitution ym-> ym e -äz

kz to get rid

of the linear y_m term and then rescaling z ® az and y_m ® by_m, putting a = 1•H2 k_z) and b2=1/(3 4 me0w 2_cH3L_{) to get} (13) ä ¶zym+ ym+1+ ym-1+ ym¤ 2 ym= 0

Derivation from a stationary action principle

The DNLS can be derived from a least action principle [4], consider the following action integral, where the notation ¶_zy = y is used:

(14) S =à âz LHy, äy-, ¶zy, ¶zäy-, zL

L =â m=1 N _ä 2Iy-m y_m- ymym * M + Hymym+1 * + y_m* ym+1L + 1 2 ym¤ 4 where ym and ä ym *

are treated as independent variables. Requiring a stationary action gives

(15) dS = dà âz LHy, äy-, ¶zy, ¶zäy-, zL =

à dz : ¶ L ¶ ym dy_m+ ¶ L ¶ äy_m* däy_m* + ¶ L ¶ y_m d y_m+ ¶ L ¶ ä y_m* dä y_m*> = = Integration by parts, assuming vanishing variation at the endpoints¤ =

à dz : ¶ L ¶ ym - â â z ¶ L ¶ y_m dy_m- ä ¶ L ¶ y_m* - â â z ¶ L ¶ y_m* däy_m*> = 0

Taking the variations to be arbitrary, one obtains the Euler-Lagrange equations of motion

(16) ¶ L ¶ ym -â â z ¶ L ¶ y_m = 0, ¶ L ¶ y_m* -â â z ¶ L ¶ y_m* = 0.

Applying Eq.16 to the Lagrangian Eq.14 results in

(17) ä ¶zym= -ym+1- ym-1- ym¤2ym ä ¶zym * = +y_m+1* + y_m-1* + ym¤ 2 y_m*

Applying the continuous transformation y_m™ yme äe

to the Lagrangian L leaves it unchanged. Consider the change in the action under such an infinitesimal transformation ym ® ym(1+äe), then

(18) dS =â m=1 N à a b dz: ¶ L ¶ ym äey_m+ ¶ L ¶ y_m* I-äeym * M + ¶ L ¶ y_mIäe ymM + ¶ L ¶ y_m* I-äe ym * M> = =â m=1 N à a b dz ¶ L ¶ ym -â â z ¶ L ¶ y_m IäeymM + ¶ L ¶ y_m* -â â z ¶ L ¶ y_m* I-äeym * M> + + ¶ L ¶ y_mIäeymM_a b + ¶ L ¶ y_m* I-äeym * M a b = 0

By observing that the terms in the integral are the Euler-Lagrange equations, and setting the expres-sion to zero because the action is invariant, it is deduced that

(19) â m=1 N _{¶ L} ¶ y_m IäeymM + ¶ L ¶ y_m* I-äeym * M = Constant

since the endpoints of integration (a,b) are arbitrary. Putting in the expression for the Lagrangian Eq.14 in Eq.19 one obtains

(20) â m=1 N _äy m * 2 IäeymM + -äy_m 2 I-äeym * M = -e â m=1 N ym¤ 2 = Constant Norm.

It turns out that all DNLS equations share this symmetry and the norm is conserved in this way via Noether’s theorem [4].

Shifting gears form the Lagrangian viewpoint to the Hamiltonian formalism, we begin by defining dynamical variables (ym , äy_m*) and their canonically conjugate pairs (pm , pm

* ) where (21) pm= ¶ L ¶ y_m and p_m* = ¶ L ¶ ä y_m*

Performing a Legendre transformation and thereby defining the Hamiltonian H produces

(22) H ºâ m=1 N y_mpm+ ä ym * p_m* - L = - â m=1 N ymym+1 * + y_m* ym+1+ 1 2 ym¤ 4

Where the last equality follows from Eq.14. Considering the differential of H;

(23) dH ºâ m=1 N HymL dpm+HpmL d ym +Iä ym * M dpm * + Hpm *_{L dä y} m * - ¶ L ¶ ym dy_m- ¶ L ¶ äy_m* däy_m* - ¶ L ¶ y_m d y_m- ¶ L ¶ ä y_m* dä y_m*

Canceling the terms with dotted differentials via the Euler-Lagrange equations Eq.14 and the definition of conjugate variables Eq.23 simplifies to

(24) dH =â m=1 N HymL dpm+Iä ym * M dpm * -HpmL dy_m-Hpm * L däym * .

Now we must proceed with extreme caution. Because of the linear relation between dynamical variables and their conjugate pairs (Eq.21), p_m = ä y_m*/2 and p_m* = -y_m/2, they can not be varied independently. Therefore we recast Eq.24 in the following form by substitution:

(25)

dH =

â

m=1 N

Hy

m

L d

ä

2

y

_m*

+

Iä y

_m*

M d -

1

2

y

_m

-

ä

2

y

_m*

dy

_m

- -

1

2

y

_m

däy

_m*

=

â

m=1 N

Iä y

m

M dy

m *

-

Iä y

_m*

M dy

_m

And finally Hamilton’s equations of motion can be read off:

(26)

¶ H

¶ y

_m

= - ä y

_m*

,

¶ H

¶ y

_m*

= ä y

_m

Which is the same pair of equations as Eq.17.

By the structure of Eq.26 it is sensible to change the previous definition of conjugate pairs and say äy_m* = pm the momentum conjugate to ym. By introducing Poisson brackets

(27)

8AHy

m

, äy

m *

_{L, BHy}

m

, äy

m *

L< = â

j=1 N

¶ A

¶ y

j

¶ B

¶ äy

*_j

-

¶ B

¶ y

j

¶ A

¶ äy

*_j

for any functions A,B of our phase space variables, we can establish the canonical Poisson bracket commutation relations (28)

8y

m

, äy

k *

_{< = d}

mk

,

8y

m

, y

m

< = 8äy

m *

, äy

_m*

< = 0.

This system is not integrable unless the number of sites N£2 [5], but there exists an instructive canonical transformation provided by [4]:

(29)

y

_m

=

N

m

e

-äqm

_{, äy}

m *

= ä

N

m

e

äqm Or by inverting Eq.29 (30)

N

m

= y

m 2

e

2 äqm

, q

m

=

1

2 ä

ln

y

_m*

y

_m

To show that this is canonical we first confirm that the Poisson bracket commutations are pre-served, by substitution of Eq.29 into Eq.27 we have

(31)

8q

m

, N

k

< =

â

j=1 N

¶ q

_m

¶ y

j

¶ N

_k

¶ äy

*_j

-

¶ N

k

¶ y

j

¶ q

_m

¶ äy

*_j

=

â

j=1 N

ä

1

2 äy

_m*

d

_jm

I2 y

_k

e

2 äqk

_d

jkM =

y

k

y

_m*

e

2 äqk

_d

mk

= d

mk

.

And checking the Jacobian determinant between the old and new variables

(32)

¶

Iy

_{m ,}

äy

_m*

M

¶

Hq

_m

, N

m

L

=

- ä

N

_m

e

-äqm 1 2 Nm

e

-äqm

-

N

m

e

äqm ä 2 Nm

e

äqm

= + 1.

This proves that Eq.29 is canonical [6]. And the Equations of motion are

(33)

¶ H

¶ N

_m

= q

_m

, -

¶ H

¶ q

_m

= N

_m

By looking at the transformation we see that N_m= ym¤2 and that qm plays the role of a phase angle.

In fact, if there would be no coupling between waveguides in our Lagrangian, these would be action-angle variables. Expressing our Hamiltonian Eq.22 in terms of the new variables we obtain

By looking at the transformation we see that N_m= ym¤ 2

and that q_m plays the role of a phase angle. In fact, if there would be no coupling between waveguides in our Lagrangian, these would be action-angle variables. Expressing our Hamiltonian Eq.22 in terms of the new variables we obtain

(34)

H

Hq, NL = â

m=1 N

- 2

N

_m

N

_m+1

cos

Hq

_m+1

- q

_m

L -

1

2

N

_m2

We see that this canonical transformation did not produce any cyclic coordinates as it would if there were no coupling between sites. But it offers intuitive guidance, for instance, from Eq.33 we have

(35)

N

m

= 2

N

m

N

m+1

sin

Hq

m+1

- q

m

L + 2 N

m-1

N

m

sin

Hq

m-1

- q

m

L

That is, the power flow to the m:th waveguide depends on the norm and the phase difference to neighboring sites. One can now define the current, Jm, as the power flow from the (m+1):st to the

m:th site (from Eq.35) as

(36)

J

m

= 2

N

m

N

m+1

sin

Hq

m+1

- q

m

L.

To sum up this section, the DNLS can be derived from a stationary action principle and has two symmetries corresponding to conservation of norm and Hamiltonian, it is therefore integrable for systems with dimension less than three [5][7]. By use of a canonical transformation a continuity equation is established and a power flow between sites defined.

Solutions and numerical investigations

There is a simple solution to the DNLS Eq .2

(37) ym= y0e

äHkz+mqL

where y₀, k and q are real. Assuming there are N waveguides with periodic boundary conditions means that eäNq=1 and q must therefore take values q_n=2pn/N. By plugging in Eq.37 in Eq.2 the following dispersion relation is obtained k = y₀2+ 2 cos qn. By Eq.36 we have that the power flow

Jm= 2 y0 2

sin qn is constant between each site for this solution. qn is the phase difference between

neighboring sites, therefore we define the topological charge TC=n=ÚArgHy_m* ym+1), summing over

all the sites [3]. The topological charge is strongly related to J_m for this solution. If n=0 we have a plane wave solution on all sites and J_m=0. For nonzero n we have a flow of power between sites and we call this a discrete vortex, discrete because we are dealing with discrete sites, vortex because Jm is circulating over the sites.

The stability of these solutions has been studied by considering small perturbations to the solution [3]

(38) y_mp =Hy0+ pmL e

äHkz+mqnL

and looking at the linearized evolution of the perturbation p_m around the original solution. Applying a perturbation of the form

(39) y_mp =Iy0+ ume

lz

+ v_m* el- zM eäHkz+mqnL and substituting into the DNLS (Eq.2) yields

(40) B-kIy0+ ume lz + v_m* el- zM + äIlume lz + l- v_m* el- zM +Iy0+ um-1e lz + v_m-1* el- zM e-äqn₊Iy 0+ um+1e lz + v_m+1* el- zM eäqn +Iy0+ ume lz + v_m* el- zM ¡y0+ ume lz + v_m* el- z¥2F eäHkz+mqnL_{= 0} – -kume lz - kv_m* el- z+ älume lz + äl- v_m* el- z +um-1e lz e-äqn_{+ v} m-1 * el- ze-äqn_{+ u} m+1e lz eäqn_{+ v} m+1 * el- zeäqn +y₀2I2 ume lz + 2 v_m* el- z+ u_m* el- z+ vme lz_{M + OIu} m 2 , v_m2M = 0

where the original solution has been canceled out and OIu_m2, v_m2M stand for terms above linear order in the um and vm modes. Collecting the terms involving l and l- one obtains the linearized

eigen-value equation for the perturbation modes:

(41) l - terms : älum+ Bum+ y0 2 v_m+ um-1e -äqn_{+ u} m+1e äqn_{= 0} l- - terms : -älvm+ Bvm+ y0 2 u_m+ vm-1e äqn_{+ v} m+1e -äqn_{= 0}

Where B = y₀2- 2 cosHq_n) in accordance with the notation used in [3]. Because this is a Hamiltonian system the eigenvalues has the property that if l is an eigenvalue, then so is -l and ±l- [4][8]. Because of this, and the choice of the form of pm (cf.Eq.39) perturbations will not grow exponentially

and the system will be stable if the eigenvalues l lie on the imaginary axis in the complex plane. This eigenvalue problem is solved numerically in this work, but it can be solved analytically [3]. By use of a Fourier transform

(42) u_m=â s U_seämqs vm=â s Vse ämqs

where qs =2ps/N, s=1,2,.. and N is the number of sites, Eq.41 is rewritten as

(43) â s I-älVs+ BVs+ y0 2 Us+ 2 VsCosHqn- qsLM e ämqs_{= 0} â s IälUs+ BUs+ y0 2 V_s+ 2 UsCosHqn+ qsLM e ämqs_{= 0.}

Looking at a specific mode HU_s, V_s), the characteristic polynomial for l becomes

(44) det -äB - 2 äCosHqn- qsL - l -äy0

2

äy₀2 äB + 2 äCosHqn+ qsL - l

= 0

Þ l = -ä2SinHqnL SinHqsL ± ä2SinHqs•2L 2 CosHqnL I2 CosHqnL Sin2Hqs•2L - y0 2_{M .}

The 4-site system

For systems of size 4k, (k=1,2,..), special charge flipping solutions of the form [9]

(45) y2 m-1=H-1L m a eäIa2z+aM y2 m= H-1L m b eäIb 2 z+ bM

satisfy the DNLS Eq.2 .These are constant norm charge flipping vortex breathers and the stability of these solutions has been studied by [3]. They are called breathers because the current J_m flows sinusoidally, “breathing”. Below we have a numerical example of this solution for initial conditions, according to Eq.45 , [y1, y2, y3, y4]=[-0.5, -0.7, 0.5, 0.7]. The solution is illustrated in Fig.4 and Fig.5.

0 10 20 30 40 50 60 z 0.2 0.4 0.6 0.8 1.0 y1¤2 , 0 10 20 30 40 50 60 z 0.2 0.4 0.6 0.8 1.0 y2¤2 , 0 10 20 30 40 50 60 z 0.2 0.4 0.6 0.8 1.0 y3¤2 , 0 10 20 30 40 50 60 z 0.2 0.4 0.6 0.8 1.0 y4¤2

0 10 20 30 40 50 60 z 0.2 0.4 0.6 0.8 1.0 y1¤2 , 0 10 20 30 40 50 60 z 0.2 0.4 0.6 0.8 1.0 y2¤2 , 0 10 20 30 40 50 60 z 0.2 0.4 0.6 0.8 1.0 y3¤2 , 0 10 20 30 40 50 60 z 0.2 0.4 0.6 0.8 1.0 y4¤2

Fig.4 Evolution of the norm in the four sites, in accordance with Eq.45.

10 20 30 40 50 60 -1.0 -0.5 0.5 1.0 Topological Charge

Current from m+1:st to m:th site

Fig.5 The current J_m (same for all sites) flows sinusoidally and the TC behaves accordingly. Starting from the original solution y_m= y0eäHkz+mqnL, where qn=p/2 (TC=1), these new solutions can

be seen as perturbations corresponding to eigenvectors with eigenvalue l=0, solving the linearized equation for the perturbation (Eq.41) yields p_m = p₀[0, e±äp•2, 0, e±äp•2]. For instance starting with the original solution with y₀=1/2 and solving for p₀2=6/25 gives the solution above with a phase-shift b = p•2 - ArctanJ2 6 /5), a=p/2 (cf.Eq.45).

We now ask ourselves if similar perturbations give rise to charge flipping solutions in the six-site system even though no analytical expression (Eq.45) is available. In particular we would like to investigate if by means of a small perturbation corresponding to a small eigenvalue one could achieve a singe flip of the topological charge.

The 6-site system

For the 6-site system with topological charge n=1, solutions ym = y0 eäHkz+mqnL are stable with respect

to small perturbations for y₀<1/2 [3], as can be read of Eq.44.

Eigenvectors corresponding to eigenvalue l=0 are of the form p_m= p₀ eia. These do not produce charge flipping solutions since they correspond to global changes in phase and amplitudes to all sites. This is clear since the phase on perturbed site m at z=0 is Arg( (y0 +p0 eäa) eämqn) and the

change in phase is Arg( y₀ + p₀ eäa ), independent of m. Similarly the perturbed amplitude is Abs( (y0 +p0 eäa) eämqn )= Abs( y0 +p0 eäa ) independent of site m.

The solution to the linearized equation for a small perturbation (Eq.41) with y0=1/16 and TC=1 is

For the 6-site system with topological charge n=1, solutions ym = y0 eHkz+mqnL are stable with respect

to small perturbations for y0<1/2 [3], as can be read of Eq.44.

Eigenvectors corresponding to eigenvalue l=0 are of the form pm= p0 eia. These do not produce

charge flipping solutions since they correspond to global changes in phase and amplitudes to all sites. This is clear since the phase on perturbed site m at z=0 is Arg( (y0 +p0 eäa) eämqn) and the

change in phase is Arg( y₀ + p₀ eäa ), independent of m. Similarly the perturbed amplitude is Abs( (y₀ +p₀ eäa) eämqn _{)= Abs( y}

0 +p0 e äa

) independent of site m.

The solution to the linearized equation for a small perturbation (Eq.41) with y0=1/16 and TC=1 is

listed in Table.1 below.

Table.1 Eigenvalues and eigenvectors for y₀=1/16 and TC=1.

Eigenvectors corresponding to the complex conjugate of an eigenvalue are phase shifted versions of the ones listed below for this set up. Note that the eigenvectors are approximations valid for small

y0. l¤ pm 0 _poeäa 0.0039... poe-ä 2 pm 3 1.0039 ... po e -äpm 3-ä 2 p 3 1.9960... 2.9909 ... po eäpHm+1L poe ä2 pm 3

We will now consider eigenvectors corresponding to eigenvalues l such that l¤ is small. For the smallest nonzero eigenvalue with p_m= poe

-ä2 pm

3 , a numerical investigation reveals that the norm ym¤2 behaves according to Fig.6.

0 10 20 30 40 50 z 0.0030 0.0035 0.0040 0.0045 0.0050 y1¤2 , 0 10 20 30 40 50 z 0.0030 0.0035 0.0040 0.0045 0.0050 y2¤2 , 0 10 20 30 40 50 z 0.005 0.010 0.015 0.020 y3¤2 p0= 8 7y0 p₀=y0 p0= 16 19y0 p0= 1 2y0 p₀=0

Fig.6 The evolution of y_m, m=1,2,3, for increasing perturbations. The phase space has been reduced such that ym= ym+3 up to a phase angle of p. Notice the different scales, when p0= y0,

y3¤2>4 y1¤2.

As can be seen in Fig.6 small oscillations have been excited corresponding to eigenvectors with higher eigenvalues. This is because the applied perturbations are large compared to y0, and have

excited other eigenmodes. The slopes in Fig.6 are slow sinusoidal oscillations as can be seen in Fig.7.

As can be seen in Fig.6 small oscillations have been excited corresponding to eigenvectors with higher eigenvalues. This is because the applied perturbations are large compared to y₀, and have excited other eigenmodes. The slopes in Fig.6 are slow sinusoidal oscillations as can be seen in Fig.7. 0 500 1000 1500 2000 z 0.0035 0.0040 0.0045 0.0050 y3¤2 0 10 20 30 40 50 z 0.0035 0.0040 0.0045 0.0050 y3¤2

Fig.7 The slowly oscillating norm for p0= 1

10y0. For such a small perturbation, the small oscillations

seen in Fig.6 have not been excited.

The small eigenvalue for this perturbation, l¤=0.0039.., implies that for small perturbations one would expect oscillations with period 2p/ l¤>1570, which can be seen in Fig.7. The slow oscillation does not have this exact period for larger values of p0 because the perturbation is large enough

compared to y₀ to excite other modes, in Fig.7 the perturbation, p₀= y0/10, is small enough to only

excite the eigenvalue mode.

The behavior of the norm seen in Fig.6-7 is naturally reflected in the current J_m. On the large order of things however, the current J_m is behaving similarly on all sites in the following way seen in Fig.8.

10 20 30 40 50 z -0.002 0.002 0.004 0.006 J p0=0 p0= 1 2y0 p0= 16 19y0 p0=y0 p₀=8 7y0

Fig.8 The current J_m (looks the same on all sites on this scale) under increasing perturbations. At p₀= y0 the current reverses its direction and the topological charge changes from +1 to -1, the

Hamiltonian has been doubled at this point. Note that because of the scale one cannot make out the small oscillations seen in Fig.4.

The applied perturbation does not send the original solution into another type of solution like Eq.45 produces the charge flipping solution for the 4-site system. Instead J_m is smoothly changing direc-tion to a soludirec-tion with TC=-1, the switch happens when the Hamiltonian is doubled. This could be expected since in the limit y₀•p₀®0 the perturbed solution tends to a vortex breather with norm p₀2 and TC=-1 (Eq.38). At p0=y0 the current Jm is oscillating around zero, and the TC flips around in

this intermediary state as illustrated in Fig.9.

10 20 30 40 50 z -0.0001 -0.00005 0.00005 0.0001 J1 10 20 30 40 50 z -0.0001 -0.00005 0.00005 0.0001 J2 10 20 30 40 50 z -0.0001 -0.00005 0.00005 0.0001 J3 10 20 30 40 50 z -1.0 -0.5 0.5 1.0 TC

Fig.9 The current Jm, m=1,2,3, at p0= y0 oscillating around 0. Since Jm=Jm+3 the total sum of the

10 20 30 40 50 z -0.0001 -0.00005 0.00005 0.0001 J1 10 20 30 40 50 z -0.0001 -0.00005 0.00005 0.0001 J2 10 20 30 40 50 z -0.0001 -0.00005 0.00005 0.0001 J3 10 20 30 40 50 z -1.0 -0.5 0.5 1.0 TC

Fig.9 The current Jm, m=1,2,3, at p0= y0 oscillating around 0. Since Jm=Jm+3 the total sum of the

currents Ú_mJm is very close to 0. The TC is flipping around 0.

The next eigenvalue produces a pm= p0 e -äpm

3-ä 2 p

3_{, a phase change of p/3 between sites instead of} 2p/3 in the last case. By Eq.38 the limiting behavior as y0<<p0 is a plane wave, TC=0. Applying this

perturbation to the system and solving numerically, in this case for y₀=1/8, yields a behavior demon-strated in Fig.10. 5 10 15 20 25 z 0.01 0.02 0.03 0.04 0.05 0.06 y1¤2 5 10 15 20 25 z 0.01 0.02 0.03 0.04 0.05 0.06 J1 p0=0 p0= 1 8y0 p0= 1 2y0 p0= 8 9y0 p0=y0 10 20 30 40 50 z -1.0 -0.5 0.5 1.0 TC TC at p0= 8 9y0

Fig.10 y₁¤2 and J₁, the graphs on other sites are similar with small phase shifts.The perturbation causes oscillations, increasing in amplitude with increasing perturbations. For small perturbations the current oscillates but is never negative, at p0=

8

9y0 the current takes on negative values and the

TC begins to flip over to 0. At y₀=p₀ the TC=0.

Fig.10 shows how a small p_m perturbs the current but not enough to reverse direction, it is not until p₀=y₀ that the TC becomes identically zero. At this point the Hamiltonian has almost been tripled. The eigenvalue of this perturbation is very close to 1, giving a period of oscillations of 2p which can be seen in the figure.

Fig.10 shows how a small p_m perturbs the current but not enough to reverse direction, it is not until p₀=y₀ that the TC becomes identically zero. At this point the Hamiltonian has almost been tripled. The eigenvalue of this perturbation is very close to 1, giving a period of oscillations of 2p which can be seen in the figure.

The previous examples were ones with small y0 (=1/16 or 1/8). Now we consider a solution with

larger y₀=0.4 and TC=1, still in the region of stability y₀<1/2. The eigenvector-perturbations p_m are no longer well approximated by the expressions in Table.1. Applying the perturbation corresponding to the smallest non-zero eigenvalue, l¤>0.17 and pm> p0 @

11 10, e -äa , eäa, 11 10 , e -äa

, eäa] with a> 2.1478, the following dynamic is obtained (Fig.11):

0 50 100 150 200 250 300 z 0.1 0.2 0.3 0.4 0.5 y12¤ 0 50 100 150 200 250 300 z 0.1 0.2 0.3 0.4 0.5 J1 p0=0 p0= 1 5y0 p0= 1 2y0 p0=0.65y0

Fig.11 Norm and current on the first site for increasing perturbations, other sites are phase-shifted versions of the first site. For small p₀ the eigenmode corresponding to l¤>0.17 is excited, as p₀ increases the oscillations become larger, with longer period. At p₀=0.65 y₀ the dynamics changes appearance, but the current is never negative and the TC never flips for any of these perturbations.

Increasing the perturbation slightly to p₀=0.675 y₀ changes the solution dramatically as seen in Fig.12: 0 20 40 60 80 100 120 140 z 0.1 0.2 0.3 0.4 0.5 0.6 y12¤ , 0 20 40 60 80 100 120 140 z 0.1 0.2 0.3 0.4 0.5 0.6 y22¤ , 0 20 40 60 80 100 120 140 z 0.1 0.2 0.3 0.4 0.5 0.6 y32¤ , 20 40 60 80 100 120 140 z -0.2 -0.1 0.1 0.2 J1 , 20 40 60 80 100 120 140 -1.0 -0.5 0.5 1.0 TC ÚmJm

Fig.12 At p0=0.675 y0 the solution changes. The current behaves like a vortex breather and the TC

is flipping with it. Phase space is reduced such that y_m¤2= y_m+3¤2. The Hamiltonian is about 1.5 times the unperturbed Hamiltonian.

This charge flipping breather solution persists with further increases in p0 but the dynamics is not

stable with larger p0. At p0= y0 the solution becomes chaotic after a while as can be seen in Fig.13:

0 20 40 60 80 100 120 140 z 0.2 0.4 0.6 0.8 1.0 1.2 1.4 y12¤ 0 20 40 60 80 100 120 140 z 0.1 0.2 0.3 0.4 0.5 0.6 y22¤ 0 20 40 60 80 100 120 140 z 0.1 0.2 0.3 0.4 0.5 0.6 y32¤

Fig.13 The dynamics are unstable and become chaotic at z>130.

The point at which chaos sets in is decreasing with increasing p0. Further numerical investigations

show that a charge flipping solution is only obtained by perturbing with the smallest non-zero eigen-value perturbation when y₀=0.4.

In the case with TC=2, solutions (Eq.38) are stable with respect to small perturbations for all ampli-tudes y₀ (by Eq.44). The smallest nonzero eigenvalue for y₀=1, l¤>0.38, corresponds to oscilla-tions with period >16.5. Perturbing with this eigenvector yields (Fig.14):

0 10 20 30 40 50 z 0.6 0.8 1.0 1.2 1.4 1.6 y12¤ , 0 10 20 30 40 50 z 1.6 1.7 1.8 1.9 2.0 2.1 2.2 J1 p0=0 p0= 1 10y0 p₀=4 10y0

Fig.14 The norm and current for increasing p0. The TC is unchanged by these perturbations. For

p₀=1

10y0 a pure eigenmode corresponding to l¤>0.38 is excited. Larger perturbations become

increasingly chaotic. : 0 10 20 30 40 50 z 0.5 1.0 1.5 2.0 2.5 3.0 y12¤ , 10 20 30 40 50 z -3 -2 -1 1 2 3 J1 , 10 20 30 40 50 -2 -1 1 2 TC >

Fig.15 The norm and current for p0= 6

10y0. Here the solution is chaotic enough for the current to

: 0 10 20 30 40 50 z 0.5 1.0 1.5 2.0 2.5 3.0 y12¤ , 10 20 30 40 50 z -3 -2 -1 1 2 3 J1 , 10 20 30 40 50 -2 -1 1 2 TC >

Fig.15 The norm and current for p₀= 6

10y0. Here the solution is chaotic enough for the current to

reverse and the TC flips around.

The chaotic behavior in Fig.15 is maintained for larger values of y₀ except it takes over for smaller perturbations. For instance, with y₀=4 the chaos sets in at about p₀= 3

10y0.

Conclusion and Discussion

The behavior of the current J_m, for perturbed solutions ym p

=Hy0+ pmL e

äHkz+mqnL _{in the 6-site system}

where pm correspond to eigenvectors with small eigenvalues does not produce a reversed current

and no changing topological charge if the initial amplitude y0 is small. This could be because the

nonlinear interaction goes like y₀3 and becomes negligible for small y₀. Instead these stable solu-tions gradually change their behavior towards the limiting solution as y₀/ p₀®0, typically the change occurs when the Hamiltonian has been more than doubled.

By increasing y₀ towards the limit of linear stability for TC=1 solutions, y₀<1/2, a large perturbation corresponding to the lowest non-zero eigenvalue produced a charge-flipping breather solution. Using this perturbation one could change the vorticity or topological charge in the six site system. However, the breather solutions are more easily controlled and stable in the 4-site system. Unfortu-nately no single charge flip was found, but it is possible that a clever superposition of perturbations could achieve a one time reversal of the current. Further investigations in this regime where y0 is

large, trying new ways to perturb the solution could give interesting results. Perhaps if one allowed for higher order nonlinearities in the model the dynamics would change and the nonlinear interac-tion could interact more with the current Jm and produce a single flip of the TC.

The behavior of the current Jm, for perturbed solutions ym=Hy0+ pmL eHkz+mqnL in the 6-site system

where p_m correspond to eigenvectors with small eigenvalues does not produce a reversed current and no changing topological charge if the initial amplitude y0 is small. This could be because the

nonlinear interaction goes like y₀3 and becomes negligible for small y₀. Instead these stable solu-tions gradually change their behavior towards the limiting solution as y₀/ p₀®0, typically the change occurs when the Hamiltonian has been more than doubled.

By increasing y₀ towards the limit of linear stability for TC=1 solutions, y₀<1/2, a large perturbation corresponding to the lowest non-zero eigenvalue produced a charge-flipping breather solution. Using this perturbation one could change the vorticity or topological charge in the six site system. However, the breather solutions are more easily controlled and stable in the 4-site system. Unfortu-nately no single charge flip was found, but it is possible that a clever superposition of perturbations could achieve a one time reversal of the current. Further investigations in this regime where y₀ is large, trying new ways to perturb the solution could give interesting results. Perhaps if one allowed for higher order nonlinearities in the model the dynamics would change and the nonlinear interac-tion could interact more with the current J_m and produce a single flip of the TC.

Acknowledgments

Magnus Johansson’s guidance has been invaluable.

References

[1] D.K. Cheng, “Field and Wave Electromagnetics”, 2:nd Ed. Addison-Wesley. 1989. [2] F. Lederer et al., “Discrete solitons in optics”, Phys. Rep. 463, 2008 1-126.

[3] A.S. Desyatnikov et al., ”All-optical discrete vortex Switch”, Phys. Rev. A, 83, 2011 063822. [4] M. Öster, “Stability and Mobility of Localized and Extended Excitations in Nonlinear Schrödinger Models”, Linköping Studies in Science and Technology. Dissertations No. 1072. 2007.

[5] J. C. Eilbeck and M. Johansson, “The Discrete Nonlinear Schrödinger Equation - 20 Years On”, Proc. of the 3rd Conf. Localization & Energy Transfer in Nonlinear Systems (June 17-21 2002, San Lorenzo de El Escorial Madrid) ed L Vázquez et al (World Scientific, New Jersey, 2003), pp 44-67 (arXiv:nlin.PS/0211049).

[6] H. Goldstein, C.Poole and J.Safko, “Classical Mechanics”, 3:rd Ed. Addison-Wesley. 2002. [7] S. M. Jensen, “The nonlinear coherent coupler”, IEEE J. Quantum Electron., QE-18, 1982, 1580-1583.

[8] M. Johansson et al., “Dynamics of breathers in discrete nonlinear Schrödinger models”, Physica D, 119 1998 115-124.

[9] L. Casetti and V.Penna, “Vortex Structures in a Chain of Coupled Bosonic Wells and the Mott Regime”,J.Low Temp. Phys, 126, 2002 455