Table Of ContentStationary scattering from a nonlinear network
Sven Gnutzmann1, Uzy Smilansky2,3, and Stanislav Derevyanko4
1School of Mathematical Sciences, University of Nottingham, Nottingham NG7 2RD, UK
2 School of Mathematics, University of Cardiff, Cardiff CF24 , UK
3 Dep. Physics of Complex Systems, Weizmann Institute of Science, Rehovot, Israel
4 Nonlinearity and Complexity Research Group, Aston University, Birmingham, UK
Transmission through a complex network of nonlinear one-dimensional leads is discussed by ex-
tendingthestationaryscatteringtheoryonquantumgraphstothenonlinearregime. Weshowthat
the existence of cycles inside the graph leads to a large number of sharp resonances that domi-
natescattering. The latterresonances arethen shown to beextremely sensitivetothenonlinearity
and display multi-stability and hysteresis. This work provides a framework for the study of light
propagation in complex optical networks.
1
1
0 Thestudyofquantumgraphshasgaineditspopularity 1 g = 0 2
2
in recent years [1] not only because graphs emulate suc- ψ + ψ 6 ψtrans
n cessfully complex mesoscopic and optical networks, but in refl
a also because they manage to reproduce universal prop- 4
J
erties (such as level statistics, transmission fluctuations
4
andothers)observedingenericquantumchaoticsystems.
Here we generalize quantum graph theory to the nonlin-
]
D ear domain. The theory will be applied in particular to 3
C show the effect of nonlinearity on transmission through
networks of nonlinear fibers. Our model may also be FIG.1. (coloronline)AgraphwithV=4vertices,B=6bonds
.
n used as a simple yet non-trivialmodel where the univer- and N=2 leads. The incoming, reflected and transmitted
li sal properties derived from detailed numerical computa- waves are shown on therespective leads.
n
tions of Bose-Einstein condensates in non-regular traps
[
[2–7] could be further investigated.
have the coordinate x [0, ) and x =0 is at the ver-
2 Scattering is studied as a stationary process. The l ∈ ∞ l
v tex where the lead is attached. The wave function on
main finding is that the sharp resonances which dom-
2 the graph is a bounded piecewise continuous and differ-
inate scattering in networks with complex connectivity
8 entiable function on the edges written collectively as
lead to a dramatic amplification of the nonlinear effects:
3
2 while non-resonant scattering hardly deviates from the Ψ(x)= ψ (x ) B+N (1)
. predictionsofthelineartheory,tuningtheparametersto { e e }e=1
1
a nearbyresonance(withoutchangingthe incomingfield where ψ (x ) is the wave functions on the edge e. While
1 e e
0 intensity) brings the system into the nonlinear regime the model can describe far more general settings we re-
1 which is signaled by multi-stability and hysteresis. For strict ourselves in this exploratory work to the discus-
: this reason we revisit the theory of scattering in the lin- sion of stationary scattering. The wave function on edge
v
i ear regime and demonstrate that sharp resonances with e satisfies the stationary nonlinear Schr¨odinger equation
X large amplification of the incoming wave inside the sys- (NLSE)
r temareveryfrequentforgraphscomparedtoothercom-
a plex(chaotic)scatteringsystems. Theoriginofthiseffect d2ψe +g ψ 2ψ =Eψ . (2)
can be related to the topology of the graph(existence of − dx2 e| e| e e
e
cycles)andleads to apower-lawdistribution forthe am-
plification. Here, ge is real nonlinear coupling parameter which we
assume constant on each edge (but it may take different
valuesondifferentedges). E istakenpositivethroughout
this manuscript,E =k2 andk reducesto thewavenum-
I. THE NONLINEAR SCHRO¨DINGER ber (propagation constant in fiber optics) in the linear
EQUATION ON GRAPHS
case. Setting g = 0 on all edges will reduce our model
e
to a standard(linear) quantum graph. Note that for ap-
Consider a general metric graph which consists of V plications in fiber optics, g ψ 2 E and the nonlinear
| | ≪
vertices connected by B internal bonds and N leads to term is a small perturbation of an otherwise linear wave
infinity as illustrated in Fig. 1. The bonds and leads equation. Then the stationary equation above describes
will be collectively referred to as edges. The bonds are the spatial evolution of the amplitude of a continuous
of finite lengths L , b = 1, ,B and are endowed with wave beam rather than a wave envelope (for which one
b
···
coordinates x [0,L ] (with a definite choice for the wouldhaveanon-stationarynonlinearSchr¨odingerequa-
b b
∈
direction in which x increases). The semi-infinite leads tion [8]).
b
2
A. The wave function on a single bond task. Since the purpose here is to display general fea-
tures of wave propagation through a nonlinear network,
The solutions of (2) on a single bond can be obtained we chose a “minimal” set of local matching conditions
[9] by writing ψ(x) = r(x)eiθ(x) (omitting the index e commonly used in the linear case: we require continuity
for the moment). Then, (2) is equivalent to two coupled and that the sum over all outgoing derivatives of wave
ordinary differential equations functions on edges adjacent to j to be proportional to
the commonvalue φ of the wave-functionat the vertex,
j
d d
=0 and =0 (3a)
ψ (0)=φ (e=1,...,v ) (7a)
dxH dxL e j j
where and
=1 dr 2+ r2 dθ 2+ Er2 gr4 (3b) vj dψ
e
H 2(cid:18)dx(cid:19) 2 (cid:18)dx(cid:19) 2 − 4 =λ φ . (7b)
dx (cid:12) j j
dθ Xe=1 e(cid:12)xe=0
=r2 . (3c) (cid:12)
L dx The constants λ (cid:12)are arbitrary real parameters. In the
j
classical particle picture, the imaginary part of this con-
These equations formally describe a classical particle in
a central potential V(r) = (E/2)r2 (g/4)r4 on the 2D dition ensures conservation of angular momenta (6) and
− the real part can be expressed via
plane where x now takes the role of time. The Hamil-
tonian energy and the angular momentum are con-
H L vj
stants of motion. Note that the angularmomentum (3c)
p (0)=λ r , (8)
e j j
reduces to the flux
Xe=1
∗
Im ψ dψ/dx (4) where p is the radial momentum associated with the
L≡ e
particle on the edge e, and r is the radial coordinate
j
carriedby the wavefunction. Forgivenvaluesfor and at the vertex. When v = 2 the matching condition is
H j
thefullsolutionisobtainedintheformoftheintegrals equivalent to replacing the vertex by a δ potential with
L
strength λ .
j
x= 2 2V 2/r2 −1/2dr (5a) In linear quantum graphs theory, it was found useful
Z (cid:0) H− −L (cid:1) to express the matching conditions in terms of a vertex
θ = r−2 2 2V 2/r2 −1/2dr (5b) scattering matrix σ(j) which connects the coefficients of
LZ H− −L incoming and outgoing waves [1]. Though in the nonlin-
(cid:0) (cid:1)
earsettingsthelackofasuperpositionprincipleprohibits
which can be reduced to elliptic integrals [9].
a decomposition into incoming and outgoing waves, the
concept can be taken over formally. Defining
B. Matching conditions at the vertices 1 dψ(0)
ain,j = kψ (0)+i (9a)
e 2k (cid:18) e dx (cid:19)
Two physical requirements guide our choice of the e xe=0
1 dψ(0)
matching conditions at the vertices: continuity and the aout,j = kψ (0) i (9b)
e 2k (cid:18) e − dx (cid:19)
conservation of flux. Consider a vertex j with vj ad- e xe=0
jacent edges and set x = 0 as the vertex coordinate
e and collecting them in vectors ain/out,(j) =
on all the edges e emanating from j. Using the clas-
T
sical point-particle analogue, continuity implies that at ain/out,j,...,ain/out,j the matching conditions
1 v
“time” x = 0, all radii r and angles θ assume the (cid:16) (cid:17)
e e e become
same values. Flux conservation is equivalent to angular
momentum conservation, aout,(j) =σ(j)ain,(j) (10a)
vj
with
=0 . (6)
e
L
Xe=1
There is a large family of mathematically acceptable σe(je)′ =v1 (cid:18)1+e−2iarctanvλjjk(cid:19)−δee′ (10b)
j
matchingconditionswhichsatisfythelatterrequirements
– e.g. all matching conditions that define a self-adjoint The flux conservation follows from the uni-
linear Schr¨odingeroperatoronthegraph(see[10])satisfy tarity of the vertex scattering matrix and
flux conservation. The matching conditions appropriate Le = Im ψe∗(0)dψe(0)/dxe = k |aoeut,j|2−|aien,j|2 so
for any particular experimental setting should in prin- that v = 0 becomes ai(cid:0)n,(j) 2 = aout,(cid:1)(j) 2
e=1Le j| | j| |
ciple be derived ab initio, which is clearly a non-trivial - impPlying flux conservatioPn. In the seqPuel we will
3
assume that the vertex matching conditions (7a) and Though the general solution of the NLSE on each edge
(7b) are satisfied on all vertices. and the matching conditions are all known explicitly,
To finish the discussion of the vertex matching con- it is generally not possible to solve the correspond-
ditions we note that matching conditions for the time ing set of equations and obtain the scattering function
dependent NLSE on star graphs where discussed previ- s(k,ain,leads)inclosedform. Therefore,weshallcontinue
ously in [11, 12]. The more relevant to the present work thediscussioninthefollowingsectionbypresentinganu-
is Ref.[11] where the authors treated rigorously the case merical solution of a relevant example.
v = 2 (see also [2]). For weak nonlinearity the vertex
j
scattering matrix (10b) follows from their derivation.
II. RESONANT SCATTERING FROM A
NONLINEAR NETWORK
C. Scattering from a nonlinear network
A. A simple examplary model
In the linear case the transport through a quantum
graph with N leads can be described by an N N uni- We study the graph shown in Fig. 1. The six bond
×
taryscatteringmatrixS(k)whichconnectsincomingand lengths were chosen by a random number generator in
outgoing amplitudes on the leads the interval 0 < L < 1 and rationally independent
b
within the numerical accuracy. We have kept the same
aout,leads =S(k)ain,leads. (11)
set of lengths for all numerics that is presented in this
The scattering matrix S(k) can be expressed explicitly manuscript (see [14] for the actual values). While we do
in terms of the vertex scattering matrices σ(j), the bond not show results for different (random) choices we have
lengths and the wave number k [13] in the form checkedthat they leadto qualitatively equivalentresults
(the statistical properties that we will mention are also
1 quantitatively equivalent).
S(k)=ρ+τ T(k)τ . (12)
out in
1 T(k)σ Two linear leads (L for ’left’ and R for ’right’) are
int
−
attachedatthevertices1and2respectively,withg =
L,R
Here T(k) is a diagonal 2B 2B matrix with diagonal
entrieseikLb thatgivethepha×sedifferenceofaplanewave 0a.inA2 isntacitdioennatrfyrowmavtehewlietfht lEea=d iks2paarntdialilnytetnrasintsymIiitnte=d
at the two ends of the bond b (each bond appears twice | L|
to the right lead and partially reflected
due to the two possible directions of a plane wave). The
matricesρ,σint,τin,andτout arebuiltupfromthevertex ψ (x )=ψ +ψ =ain e−ikxL +r(k,ain) eikxL
scattering matrices σ(j). I.e. ρ is an N N matrix that L L in refl L L
× (cid:0) (17(cid:1)a)
contains all direct scattering amplitudes (if all leads are
attached to different vertices this is a diagonal matrix); ψR(xR)=ψtrans =aiLn t(k,aiLn) eikxR. (17b)
the2B 2Bmatrixσ containsallscatteringamplitudes
int
from o×ne (directed) bond to another inside the graph; Global gauge invariance implies that the phase of aiLn
ethvaenttcuoanlltyaiτninsacnadtteτroiuntgaraem2pBlit×udNesanfrdomN×th2eBlemadastriincteos Tcahnebreeflcehctoisoennaanrdbittrraanrsilmy,isssoionwecoteaffikecieaniLnts=r(k√,IaiiLnn)≥an0d.
the bonds, and from the bonds out to the leads. They t(k,aiLn)willbecomputedasfunctionsofkandIin. Inthe
can be combined to one unitary matrix linearcase(ge =0forallbonds)thereflectionandtrans-
missionamplitudesarethematrixelementsofa2 2scat-
×
Σ= ρ τout . (13) tering matrix S(k). In the nonlinear case they form the
(cid:18)τin σint(cid:19) two components of the scattering (vector valued) func-
tion s=ain(r(k),t(k))T. Flux conservation implies
UnitarityofΣimpliestheunitarityofthescatteringma- L
trix S(k) The unitarity of S(k) ensures global flux con-
r(k,ain)2+ t(k,ain)2 =1 . (18)
servation | L | | L |
aout,leads 2 = ain,leads 2. (14) In the present setting the importance of the nonlinear
| | | | effects is controlledby I . The linear theory is obtained
in
Moreover, in a linear system the scattering matrix S(k) inthelimitI 0. However,thestrengthofthenonlin-
in
is independent of the incoming amplitudes ain,leads. earity will not b→e uniform as a function of k because the
In the nonlinear case transport is described by a N- intensity inside the graph structure may vary strongly,
component nonlinear scattering function especially near resonances as we will show below.
Using the formalism described above, the solution of
aout,leads =s(k,ain,leads) . (15)
the scattering problem reduces to a finite set of non-
linear equations in a high-dimensional space which re-
Flux conservation on each vertex implies that the scat-
quires rather substantial computer resources. The nu-
tering function conserves the norm
merical complexity can be alleviated further by consid-
s(k,ain,leads)2 = ain,leads 2. (16) ering the special case with just one nonlinear bond˜b, so
| | | |
4
that g = δ . In the present simulation we choose for the bond (1,2) in the linear limit. While the in-
b ± b,˜b
˜b = (1,2) but we have also checked that other choices tensity on the bond fluctuates around the incoming in-
tensity, I , there are also large peaks at narrow res-
yield quantitatively similar behavior. in
onances. For example, near the marked sharp reso-
nance in Fig.2 the intensity on the bond is two orders
of magnitude ( 320 times) higher than the intensity
B. Linear scattering: resonances and amplification ≈
of the incoming beam. Over a larger spectral interval
(0 < k < 20000) we found several other resonances
The key to the understanding of the amplification of with amplification factors α > 105 and a distribution
the nonlinear effects resides with the scattering in the P(α) = K−1 Kδ(α α(k))dk with a power law decay,
linear regime. The upper panel of Fig. 2 shows the (lin- P(α) α−s wR0ith s −2.85. The algebraic decay of this
ear) reflection probability r(k,0)2 as a function of the ∼ ≈
| |
wave number k for an interval of moderate wave num-
1
bers where the typical length of a bond is about 14-
16 wave lengths. Several resonances are clearly visible.
Quantumgraphswith incommensurate bond lengths are 10−4
a paradigm of quantum chaotic scattering [13] where
)
the statistics of resonances (i.e. of their location and α
(
their widths) is very close to the universal predictions P 10−8
of random-matrix theory. This implies that the width
of the resonances is distributed over a broad range of
values as illustrated by Fig. 2. The width of a single 10−12
resonance is inversely proportional to the decay time of 10−2 1 102 104
α
thecorrespondingresonantstate. Narrowresonancesare
associatedwithwaveswhicharetrappedinthestructure
FIG. 3. (color online) Double logarithmic plot of the numer-
forlongtime,whichisexpressedinthestationaryformal-
ically obtained distribution oftheamplification factor onthe
ismbyrelativelylargevaluesofthe wavefunctiononthe
nonlinear bond for thegraph depicted in figure1.
bonds. ThelowerpanelinFig.2showsthe amplification
distribution is a special feature of networks which would
1 usually notbe expected in other chaotic scatteringmod-
elssuchasscatteringthroughachaoticquantumdot. To
explainthedifferencebetweenanetworkandamoregen-
2
)| eral chaotic scattering system let us go back to Eq. (12)
k,0 0.5 which describes scattering through the network. Simi-
(
r
| larequations havebeenusedto modelchaoticscattering
(e.g. by an average over the matrix Σ from Eq. (13)).
0 In either case the factor (1 T(k)σint)−1 in the second
−
term is responsible for the resonances and one may ex-
100
pect large amplification factors whenever T(k)σ has
int
10 an eigenvalue near unity (an exact eigenvalue unity in-
dicates the existence of a bound state which is confined
α
1 to the bonds). In any generic model of chaotic scatter-
ingunimodulareigenvaluesofthesubunitarymatrixσ
int
0.1 (which acts on a vector of 2B coefficients, one for each
directedbond)arestronglysuppressed. Forotherknown
180 185 190 195 200
examples of resonant scattering in one-dimensional non-
k
linearSchr¨odingersystems [5–7],the equivalentofσ is
int
FIG. 2. (color online) Upper panel: reflection probability one number with modulus smaller than 1. In all these
|r(k,0)|2 forthegraphdepictedinfigure1inthelinearlimit. modelshighamplificationfactorsareeithercut-offorex-
The arrow marks thenarrow resonance which is discussed in tremely rare. However for networks with the standard
the text. Lower panel: intensity amplification factor α(k) vertex matching conditions (7a) and (7b) (and λ = 0 –
in logarithmic scale in the linear limit aiLn → 0 on the bond see below for λ = 0) the situation is drasticallyjdiffer-
j
connecting vertices 1 and 2 in figure1. 6
entas everycycle createdfromthe bonds ofthe network
supports an eigenvalue unity of σ (moreover cycles of
int
factor even length support eigenvalues minus unity). In fact,
let us consider a cycle that consists of three bonds b ,
1
Lb ψ (x )2dx b , and b , then corresponding eigenvector a of σ with
α(k)= 0 | b b | b (19) 2 3 int
R L I uniteigenvaluevanishesonalldirectedbondsthatdonot
b in
5
belong to the cycle, and has values 1 on the directed 1 α
±
bonds thatbelongtothe cycle(the twosignscorrespond 300
to two different directions to go through the cycle). For
200
rationally dependent bond lengths this implies that one
0.75
may chose k such that 100
2)|
eikLb1 =eikLb2 =eikLb3 =1 (20) inrk,a(L 0.5 197.92 k 197.922
|
which shows the existence of embedded bound states in
the continuum of scattering states. The construction is
equivalent to that of perfectly scarred states (see [16]) 0.25
– these scarred states vanish exactly on the vertex and
197.918 197.92 197.922
for each bond on the cycle the bond length is an integer k
multiple of the wave length 2π/k. For incommensurate
bond lengths, there are no perfectly scarred states on FIG. 4. (color online) Reflection probability near thenarrow
resonancemarkedinFig.2. Thecentral(blue)curveistheres-
the cycle as the condition (20) can neverbe met exactly.
onance in the linear limit. The (green/red) curves left/right
Howeverthemappingk (eikLb1,eikLb2,eikLb3)isaner-
→ from thecentral one correspond to 8 increasing values of the
godicflowona3-dimensionaltorus,onethusfindsvalues incomingintensity|aiLn|2(inequalstepsfrom0.0006to0.0048)
forthewavenumberk whichapproximatecondition(20)
in the attractive/repulsive case (g = 1/ g = −1). The inset
to arbitrary precision. In the exemplary model we have shows the corresponding amplification factor.
usedforourcalculationthegraphstructurecontainstwo
independentcyclesthatbothcontainthenonlinearbond.
The two cycles can be chosen such that they consist of
Fig. 4 which resolves the narrow resonance in the reflec-
three bonds.
tion probability (and the amplification factor) for vari-
In the above discussion we have assumed the the ver-
ous values of the incoming intensity. In the attractive
tex potentials λ vanish. The vertex scattering matrices
j (repulsive) case the resonance moves to the left (right)
(10b) show however that for sufficiently high wave num-
as I increases. However some parts of the curve move
in
ber k these potentials are not relevant. Note also, that
fasterthanotherswhicheventuallyleadstoamultivalued
vertex matching conditions which are entirely different
dependency above the critical value I . This implies
crit
from (7a) and (7b) do not necessarilyhave a similar dis-
multistability – an experiment would show hysteresis.
tribution of narrow resonances.
In both the attractive and the repulsive cases the crit-
icalincomingintensitywheremultistabilitysetinisnear
I 0.002. For this incoming intensity the strength of
crit
C. Implications for nonlinear scattering: ≈
nonlinearity inside the graphis typically (i.e. away from
multistability the resonance) on the order of ν 5 10−8 – at the
typ
resonanceit is howeverν =1.5 1∼0−5×. These findings
res
The strength of nonlinearity on the nonlinear bond are similar to previous work on n×onlinear resonant scat-
may be measured by the effective parameter teringfromquantumdots[3,4]andfromone-dimensional
structures [5–7]. Our model generalisesthe latter results
g g Le
ν = | e| ψ 2 | e| ψ (x )2dx . (21) by allowing for additional topological complexity.
e e e
E | | ≡ EL Z | |
e 0
For a fixed incoming intensity I , and with g = 1 the
in e
effective nonlinearity will proportional to the| a|mplifica- D. Application in nonlinear fiber optics
tion factor
The numerically found power-law distribution of the
1
ν(E)= α(E)I . (22) amplification factor α can be expected to be a generic
in
E
feature (at least for similar types of vertex matching
Eveniftheincomingintensityistoolowtoinducenotice- conditions). This implies that by tuning the parameters
ablenonlineareffectsoffresonance,atnarrowresonances of an experiment to a sufficiently narrow resonance one
the high fields on the bonds are expected to behave in may find arbitrarily high amplification factors and thus
a nonlinear way. This qualitative picture is supported the nonlinearity effects of multistability and hysteresis
by the numerical simulations. For incoming intensities may be observed at considerably lower incoming veloci-
up to I 0.005 the nonlinearity is either not relevant tiestheninourexample. Nonethelessletustranslatethe
in
≈
at all or can be taken into account as a perturbation for importantparametersofourmodeltoafiber-networkex-
almosttheentirek-spectrum. However,nearthe marked perimental setup.
resonancetheamplificationbytwoordersofmagnitudeis For a CW optical beam propagating in a single mode
sufficienttogiverisetostrongnonlineareffectsthatcan- telecom fiber with the linear refraction index n = 1.5,
0
not be described perturbatively. This is shown in figure nonlinear Kerr coefficient n = 2.4 10−16cm2/W, and
2
×
6
effective mode area S = 50µm2 operating at the tele- onemayequivalentlychangethelengthsoftheedgesina
eff
com wavelength λ = 2πc/ω = 1.55µ the strength of the controlledwaybyslowlyvaryingthetemperature. Letus
nonlinearity(onthenonlinearbond)canbeestimatedas also mention that multistability in scattering from non-
[8] linear cristals has been observed experimentally [15].
8n ω2
0
ν = n P . (23)
c2β2S 2 ave
eff
III. CONCLUSION
Here, β is the propagationconstantandP is the bond
ave
averaged power ψ 2 of the beam in watts. For the pur-
| | To conclude, the theory presented here shows how
poses of the current estimate one can neglect mode dis-
the interplay between complex topology and nonlinear-
persion of the fiber and assume β n (2π/λ). Then
≈ 0 ity gives rise to a pronounced amplification of nonlinear
from(23)itfollowsthatinordertoachievetheresonance
levelofnonlinearityν =1.5 10−5 fromthenumerical effects. We would like to stress that while our model
res × is highly idealized, the strong amplification of intensity
example considered above at the telecom wavelength λ,
near narrow resonances is a universal effect that can be
the required power level must be as high as P 6kW
ave ∼ expectedinanylinearcomplexnetwork. Anyco-existing
which is above the thresholds for stimulated Raman and
nonlinearity that may be negligible off resonance will be
Brillouin scattering for the fiber length of a few meters.
drasticallyamplifiedat a resonance. We believe that the
Thiscanbeoffsetinseveralways: onecanconsidermode
latter effect can be observed in actual experiments with
dispersion and operate at lower frequencies so that the
interconnectedopticalfiberseventhoughthemodelitself
effective mode index βc/ω is lower; or, reduce the effec-
mayneedfurtheradjustment tofit the detailsofsuchan
tive mode area by a factor up to 10 by changing the size
experiment. Moreover we believe that NLSE on metric
of the fiber core or else one could use highly nonlinear
graphs as presented here will be a very useful paradigm
non-silica fibers where the value of the nonlinear coeffi-
system where the interplay between topology and non-
cient n can be enhanced up to two orders of magnitude
2 linearity can be studied qualitatively.
[8]. Thus, using the estimates based on the present sim-
ACKNOWLEDGMENTS
ulation,onecanexpectthatnonlineareffectswillappear
at power levels of 1 10W.
÷
In an experiment varying the wave number in a con- We wouldliketo thank Y.Silberberg,N. Davidson,P.
trolled way may not be feasible – so let us mention that Schlagheck, and T. Kottos for fruitful discussions.
[1] T.Kottos, U.Smilansky,Ann. Phys.274, 76(1999); S. [10] V. Kostrykin,R. Schrader,J. Phys.A 32, 595 (1999).
Gnutzmann,U.Smilansky,Adv.inPhys.55,527(2006). [11] J.Holmer,J.MarzuolaandM.Zworski,Commun.Math.
[2] P.Leboeuf, N. Pavloff, Phys. Rev.A 64, 033602 (2001). Phys. 274, 187 (2007)
[3] T. Paul, K. Richter,P.Schlagheck, Phys.Rev.Lett. 94, [12] Z. Sobirov, D. Matrasulov, K. Sabirov, S. Sawada, K.
020404 (2005). Nakamura, Phys. Rev.E 81, 066602 (2010).
[4] T.Paul,PLeboeuf,N.Pavloff,K.Richter,P.Schlagheck [13] T.Kottos,U.Smilansky,Phys.Rev.Lett.85,968(2000).
Phys.Rev.A 72, 063621 (2005). [14] Choiceofbondlengthsforthenumericspresentedinthis
[5] K. Rapedius, D. Witthaut, H.J. Korsch, Phys. Rev. A manuscript: L12 ≈ 0.8412, L13 ≈ 0.4429, L14 ≈ 0.4142,
73, 033608 (2006). L23 ≈ 0.7483, L24 ≈ 0.5137, L34 ≈ 0.5103 (up to order-
[6] K. Rapedius, H.J. Korsch, Phys. Rev. A 77, 063610 ing these arethefirst six numbersgiven byoctave using
(2008). the commands rand(’seed’,13); rand(1,6)).
[7] K.Rapedius,H.J.Korsch,J.Phys.A41,355001(2008). [15] P.E. Barclay, K. Srinivasan, O. Painter, Optics Express,
[8] G.P.Agrawal,TheNonlinearFiberOptics,4-thed,(Aca- 13, 801 (2005).
demic Press, 2007). [16] H. Schanz, T. Kottos, Phys. Rev. Lett. 90, 234101
[9] L.D. Carr, C.W. Clark, and W.P. Reinhardt Phys. Rev. (2003).
A 62, 063610 (2000); Phys. Rev.A 62, 063611 (2000).
