Constantine N Alexeyev

Alexander V Volyar

Maxim A Yavorsky

orbital angular momentum, and discrete current in a double-ring fiber array

Keywords: number(s): 42, 81, Qb

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

The study of propagation of light in fiber arrays dates back to a theoretical paper by Jones, in which he studied an infinite one-dimensional array of monomode evanescently coupled optical fibers [1]. Since then the physics of infinite fiber arrays has been studied in detail for the cases of one and two dimensions. Along with elaboration of the models of coupling [2] a large amount of attention has been paid to conceptual questions concerned with a novel type of diffraction-the so-called discrete diffraction-that takes place in waveguide arrays and opens new vistas in molding the flow of light [3]. Recently impressive results have been obtained while studying curved waveguide arrays [4] and effects of nonlinearity in arrays [5]. The physics of fiber arrays has proved to be so profound as to enable demonstration of a number of physical effects, initially discovered in other fields not immediately connected with waveguide optics [6]. The present interest in fiber arrays is also concerned with their applications for creation of beams with subwavelength resolution [7] and their use for manipulating quantum states of light [8].

Among fiber arrays one can single out an important class of one-dimensional finite arrays with periodic boundary conditions, which corresponds to the so-called circular fiber arrays. The question of light's propagation in such arrays with various geometries has been addressed in a number of papers, including both classical works [9,[START_REF] Snyder | Optical Waveguide Theory[END_REF] and more recent studies [START_REF] Hizanidis | [END_REF]. Practical interest in such arrays is concerned with their application for creation of cylindrical beams [12]. Recently, the question of circular fiber arrays has drawn attention in connection with orbital angular momentum (OAM) generation by such arrays [13] and its evolution in them [14]. Recently Desyatnikov et al. studied the behavior of discrete vortex solitons and vortex breathers in a circular array of nonlinear waveguides [15]. Besides, circular fiber arrays have the structure typical for some types of photonic crystal fibers (PCFs) [START_REF] Joannopoulos | Photonic Crystals: Moulding the Flow of Light[END_REF]. In this connection studying their properties may shed light on certain properties of light's propagation in PCFs. * c.alexeyev@yandex.ua However, among circular fiber arrays only one type of a single-ring array has been studied to date. The only modification of this model concerns circular arrays with one central fiber [START_REF] Snyder | Optical Waveguide Theory[END_REF]. Meanwhile, it is of interest to generalize the model of a circular array to the case of N concentric arrays.

In this paper we propose the method of solving the problem of modes of concentric evanescently coupled arrays, which consist of the same number of identical fibers. This method is illustrated on the example of a double-ring circular array. The first aim of our paper is to obtain the modes and the spectrum of their propagation constants in a double-ring array of identical monomode weakly guiding fibers. We study energy transfer between the rings and determine the conditions at which the energy is mostly confined within them. The second aim of this paper is to study the OAM of supermodes of such arrays and the influence of array's parameters on the value of specific OAM carried by the mode.

II. SUPERMODES OF A DOUBLE-RING ARRAY

Consider a double-ring fiber array, which consists of two concentric evanescently coupled circular arrays each composed of N identical monomode fibers. The centers of the fibers are positioned at the vertices of homothetical N -gones as is shown in Fig. 1. The coupling constant for the inner ring is I , coupling in the outer ring is denoted by J , and coupling between the nearest fibers, which belong to different rings, is K. This asymmetric choice of exchange constants allows one to describe all possible configurations of the system in study. There are two equivalent ways to study the modes of coupled arrays. The first one is connected with composing the set of coupled-mode equations in envelope amplitudes a i (z) of the ith fiber field [1]. The second approach, outlined in [START_REF] Snyder | Optical Waveguide Theory[END_REF] and recently developed in [17,18], implies constructing the perturbation matrix of the total Hamiltonian-like operator that describes propagation of light in the array. In the following we will limit our consideration to the frameworks of the second method.

According to [17], propagation of light in an array of coupled fibers is governed in the scalar approximation by the waveguide equation: where ∇ t = ( ∂ ∂x , ∂ ∂y ), k is the wave number in vacuum, β is the propagation constant, n cl is the refractive index of the medium, which the array is embedded in, and e t is the transverse part of the electric field introduced as E(x,y,z) = e(x,y)e i βz (see notations in [START_REF] Snyder | Optical Waveguide Theory[END_REF]). Perturbation V caused by the presence of the array can be written for identical step-index fibers as V = 2k 2 n 2 co i δ(ρr i ), where = (n 2 con 2 cl )/2n 2 co , ρ is the core's radius, and n co is its refractive index. It is assumed here that the ith fiber is associated with the local polar coordinate frame, in which r i coordinate is defined. For infinitely spaced fibers the solution of Eq. ( 1) can be represented as a linear superposition of scalar solutions ψ(r), localized at the ith fiber: |e t = i C i |i , where |i = col(0,0...ψ(r)...0,0) and ψ stands in the ith place. As the separation of the fibers diminishes the modes of the array should be determined through the standard procedure accepted in quantum mechanics [START_REF] Davydov | Quantum Mechanics[END_REF]: One should average Eq. ( 1) over the basis {|i }. Then the structure of modes is obtained through solving the eigenvector equation for the obtained matrix. As can be shown [17], for a one-dimensional infinite array of identical fibers the obtained matrix in the nearest-neighbor approximation has the following elements:

∇ 2 t + k 2 n 2 cl + V e t = β2 e t , ( 1
H ik = β2 δ ik + C(δ i,k+1 + δ i+1,k ), ( 2
)
where C is a coupling constant expressed through the standard overlap integral. Imposing a periodic boundary condition |i = |i + N one can describe with matrix elements (2) a circular array of N fibers. The first term on the right of Eq. ( 2) is usually omitted since it does not affect the mode structure.

Here we also omit the dependence of constant β on the radial number m.

A. Spectrum and supermodes

For a double-ring array of 2N fibers the dimension of the H matrix should be 2N. It is easily verified that for the numeration of fibers used in Fig. 1 the corresponding perturbation matrix H 2N , which determines the structure of modes, is given by the following expression:

H 2N = I M N K1 N K1 N J M N , ( 3
)
where 1 N is the unity matrix of rank N and (M N ) ik = δ i,k+1 + δ i+1,k . The problem of finding the supermodes of the doublering array is equivalent to the problem of diagonalization of matrix H 2N . Fortunately, in this case the latter problem can be easily solved. To this end one should recall that the matrix M N is diagonalized by a lattice Fourier transform expressed through the diagonalizing unitary transformation:

MN = U N M N U + N , ( 4
)
where the elements of rank-N matrix

U N are (U N ) kl = 1 √ N exp(iϕ k l) and ϕ k = 2πk N . If one defines a rank-2N matrix C 2N = ( U N 0 N 0 N U N )
, where a rank-N matrix 0 N has zero elements, then

C 2N H 2N C + 2N = I D N K1 N K1 N J D N . ( 5
)
Here a diagonal matrix D N is composed of the spectrum of a circular array of N fibers, so that

(D N ) nk = 2 cos ϕ n δ nk . ( 6
)
Note that here there is no summation over n.

To diagonalize the matrix (5) it is useful to notice that it consists of diagonal matrices, which commute with each other. Moreover, a function of a diagonal matrix D is also a diagonal matrix, whose elements are the same functions of the elements of D:

(f (D)) ik = f (D ik ). (7)
These properties enable one to operate with D matrices as with C numbers. Then the problem of diagonalization of the matrix Eq. ( 5) reduces to the one of diagonalization of the matrix α ≡ ( a b b c ) composed of real numbers. The latter is diagonalized by the transformation

α = P αP + = 1 2 (a + c) + R 0 0 -R , ( 8
)
where

P = 1 2 √ R (2R -a + c) 2b 2R -a + c -2b 2R + a -c (9)
and

R = 0.25(a -c) 2 + b 2 .
One has a formal correspondence: I D N → a, J D N → c, K1 N → b, which enables us to construct the matrix P 2N of the transformation, which diagonalizes the matrix Eq. (5). Then the resulting diagonal matrix will have the structure dictated by Eq. ( 8):

H2N = P 2N C 2N H 2N C + 2N P + 2N = ⎛ ⎝ I +J 2 D N + I -J 2 2 D 2 N + K 2 1 N 0 N 0 N I +J 2 D N - I -J 2 2 D 2 N + K 2 1 N ⎞ ⎠ . ( 10
)
Having taken the functions of diagonal commuting matrices with the use of Eq. ( 7) one can obtain the spectrum of propagation constants of supermodes of a double-ring array:

β 2 m = β2 + (I + J ) cos ϕ m ± (I -J ) 2 cos 2 ϕ m + K 2 . ( 11
)
The first term here originates from the first term on the right of Eq. (2). It is interesting to study the limiting case of symmetric exchange: J = I = K, which corresponds to two parallel chains of fibers. The spectrum correction in Eq. ( 11) then turns into

β 2 m -β2 = J (2 cos ϕ m ± 1). ( 12
)
The structure of eigenmodes is obtained through the standard procedure. The state of the double-ring array is described by a 2N -dimensional vector X 2N = col(q 1 ...q N , r 1 ...r N ) ≡ col( Q N , R N ), where q i describes the amplitude of the ith fiber in the inner ring, and r i stands for the amplitude of the ith outer ring fiber. As is evident, the eigenmodes of the system are determined by the elements of the vector P 2N C 2N X 2N . It is somewhat convenient to express these eigenmodes in terms of eigenmodes of noninteracting inner and outer rings X m and Y m correspondingly [14]:

X m = 1 √ N N-1 k=0 exp(iϕ m k)q k , ( 13
)
Y m = 1 √ N N-1 k=0 exp(iϕ m k)r k , where m = 0,1 • • • N -1.
Then after a little algebra the expressions for the modes can be brought to the following form:

Z (+) m = cos α m X m + sin α m Y m , ( 14
) Z (-) m = -sin α m X m + cos α m Y m , where tan 2α m = K 4(I -J ) sec ϕ m , 0 < α m < π
2 and no summation over m is implied. The spectra of "+" and "-"labeled supermodes Z (±) m are obtained, correspondingly, from + and branches of Eq. [START_REF] Hizanidis | [END_REF]:

β m = β + (I + J ) 2 β cos ϕ m ± I -J 2 β 2 cos 2 ϕ m + (K/2 β) 2 , ( 15
)
where we have allowed for β I,J,K. The modes in Eq. ( 14) take on the simplest form in the case of a symmetric exchange J = I = K. In this case they are given by symmetric and antisymmetric combinations of single-ring supermodes:

Z (±) m = 1 √ 2 (X m ± Y m ).
Evidently, this form conveys the idea of in-phase and out-of-phase propagation of individual supermodes.

B. Multiring arrays

The scheme suggested can be easily generalized to the case of multiring arrays. Indeed, consider an array of three concentric circular arrays. The coupling constants within the circular arrays are I , J , and G, between the inner and the middle arrays K 1 and between the outer and the middle arrays, K 2 . Using the scheme of numeration analogous to the one given in Fig. 1 one can easily obtain the structure of perturbation matrix for such an array:

H 3N = ⎛ ⎜ ⎝ I M N K 1 1 N 0 N K 1 1 N J M N K 2 1 N 0 N K 2 1 N GM N ⎞ ⎟ ⎠ . ( 16
)
Using the diagonalizing matrix C 3N = diag(U N ,U N ,U N ), where "diag" stands for a diagonal matrix, one can bring H 3N to the form

C 3N H 3N C + 3N = ⎛ ⎜ ⎝ I D N K 1 1 N 0 N K 1 1 N J D N K 2 1 N 0 N K 2 1 N GD N ⎞ ⎟ ⎠ . ( 17
)
Analogously, the problem of diagonalization of matrix Eq. ( 17) reduces to diagonalization of a rank-3 matrix h = (

I K 1 0 K 1 J K 2 0 K 2 G
). The spectrum of H 3N can be obtained from the spectrum of h in the analogous manner. If the ith spectrum branch of h has the form F i (I,J,G,K 1 ,K 2 ), i = 1, 2, 3, it corresponds to the set of N eigenvalues:

β 2 mi = β2 + F i (2I cos ϕ m ,2J cos ϕ m ,2G cos ϕ m ,K 1 ,K 2 ). ( 18
)
The eigenvectors of H 3N are also obtained through the eigenvectors of h.

If x i = col(a i ,b i ,c i ) is the eigenvector of h then the eigenvector mi of H 3N reads as mi = a i X m + b i Y m + c i Z m , ( 19
)
where we use designations of Eq. ( 14) and

Z m = 1 √ N N-1
k=0 exp(iϕ m k)s k , s k being the amplitude of the kth fiber in the outer ring. As is seen, the problem of eigenmodes is reduced to the standard problem of eigenvectors of a finite-dimension matrix. The example provided describes a regular procedure for finding eigenmodes of an N -ring circular fiber array.

III. ENERGY TRANSFER IN A DOUBLE-RING ARRAY

The problem of energy transfer management in circular fiber arrays has always been in the focus of attention of communication optics, especially in connection with the search for the ways of all-optical switching and routing [START_REF] Hizanidis | [END_REF][START_REF] Hudgings | [END_REF]. In a general context this can be seen as the steps towards implementation of the so-called guidonics [21]. In that regard it is desirable to study the process of energy evolution in double-ring arrays on the basis of a rigorous theory. The interest to this question is also justified due to the fact that certain types of PCFs [START_REF] Bjarklev | Photonic Crystal Fibres[END_REF] feature the presence of pronounced ring structures [START_REF] Hu | [END_REF]24]. In this section we study numerically the spreading of energy in the array in study. This process is strongly determined by geometrical parameters of the array, which influence the values of the coupling constants I,J,K. To obtain the exact expressions for these constants let us assume that the field of the fundamental mode located on a fiber is represented in the Gaussian approximation as [START_REF] Snyder | Optical Waveguide Theory[END_REF]

= E exp - r 2 2r 2 0 , ( 20
)
where E is some constant and r 0 is the effective transverse size of this mode. Here r is supposed to be measured from the center of the fiber. As can be easily shown using the well-known results [START_REF]Selected Papers on Coupled-Mode Theory in Guided-Wave Optics[END_REF], for two monomode fibers the coupling constant C can be expressed through the following overlap integral:

C = 4k 2 n 2 co exp - L 2 2r 2 0 1 0 x exp(-x 2 )J 0 i L r 0 x dx, ( 21
)
where L is the distance between the centers of the fibers and J 0 is the Bessel function. To calculate the exchange constants I , J , and K with the use of Eq. ( 21) one has to set in it the parameter L to 2r 1 sin(π/N), 2r 2 sin(π/N), and r 2r 1 , correspondingly, where r 1,2 are the radii of the rings (see Fig. 1). The plots of the coupling constants are shown in Fig. 2.

To study the evolution of an initial energy distribution created at the input end of the double-ring array one should notice that in the main approximation the energy is supposed to be localized at the sites of the lattice. Its value is determined by the squared modulus of the field amplitude at the fiber. It is therefore convenient to characterize the energy distribution pattern with the set of such squared amplitudes. Let us describe the amplitudes at z = 0 through two N -dimensional vectors

Q N0 = col(q 10 • • • q N0 ) and R N0 = col(r 10 • • • r N0 )
, where q i0 and r i0 stand for the amplitudes of the ith fiber in the inner and outer rings, correspondingly. Using Eqs. ( 14) and ( 15) one can easily obtain for the resulting field (z) the following expression:

(z) = N-1 k=0 (A k q k + B k r k ), (22)
where the field amplitudes A k , B k at the kth fiber of the inner and the outer rings, correspondingly, are given by

A k = 1 N N-1 l,m=0 exp i ϕ l k -ϕ m l + I + J 2 β z cos ϕ l × q m0 cos R l z 2 β + i cos 2α l sin R l z 2 β + ir m0 sin 2α l sin R l z 2 β , ( 23
)
B k = 1 N N-1 l,m=0 exp i ϕ l k -ϕ m l + I + J 2 β z cos ϕ l × iq m0 sin 2α l sin R l z 2 β + r m0 cos R l z 2 β -i cos 2α l sin R l z 2 β , ( 24
)
and R l = (I -J ) 2 cos 2 ϕ l + K 2 . Figure 3 presents examples of intensity distribution evolution if one excites with equal weights the group of neighboring fibers of the inner ring at the input end of the array. As is seen, the overall evolution of intensity features a rather complicated behavior, showing intricate diffraction patterns [Fig. 3(a)]. As the radius of the inner ring r 1 increases at a constant spacing between the rings the picture of evolution somewhat simplifies, provided the scale along the longitudinal axis is the same [Fig. 3(b)]. Finally, the patterns turn into stratified columns [Fig. 3(c)]. It should be stressed that in those figures the intensity distribution in the outer ring is plotted in relative units since in absolute units the maximal intensity in the outer ring turns out to be much less than the one in the inner ring for Figs. 3(a) and 3(b). To demonstrate this fact it is sufficient to study the total energy W 2 = N k=1 |B k | 2 stored in the second ring. Averaged energy W2 ≡ W 2 z , where broken brackets stand for averaging over z, can be brought to the form

W 2 z = 1 2N N-1 l=0 |V l | 2 sin 2α l - 1 4 sin 4α l (V * l U l + c.c.) + 1 2 (3 + cos 4α l )|U l | 2 , ( 25
)
where

V l = N-1 m=0 exp(-ilϕ m )q m0 and U l = N-1
m=0 exp(-ilϕ m )r m0 are decomposition coefficients of the input field distribution over supermodes X l and Y l of the inner and outer rings. Numerical simulations show that as r 1 increases at a constant spacing between the rings the average energy in the outer ring (if one excites the inner ring fibers) also increases (see Fig. 4). Finally, it becomes equal to the average energy of the inner ring (2.5 in a.u.).

To explain such features of energy transfer let us analyze the influence of geometrical parameters of the array on the spectra of supermodes. Using Eqs. (15) and ( 21) one can plot the spectral curves as functions of r 1 (see Fig. 5). There are two main features of this family of curves, which for convenience we plotted in separate figures. First, there are numerous intersections between the curves at smaller values of r 1 , presented in Fig. 5(a) (had we plotted all the curves in the same figure, there would have been many more intersection points). Such events are associated with the so-called accidental degeneracy. In general, this is a desirable property, which helps one to reduce intermodal dispersion. Though the main efforts are made to reduce the chromatic dispersion in PCFs [START_REF] Reeves | [END_REF], the FIG. 4. (Color online) The total averaged over z energy W2 (in arbitrary units) stored in the second ring vs reduced radius of the first ring r 1 /r 0 . It is assumed that the fibers 9-13 of the inner ring are excited with equal weights and the spacing between the rings is constant. The other parameters of the array are the same as in Fig. 3. Note that the initial value of energy (in a.u.) in the inner array is 5 (five fibers excited with unity weights). question of intermodal dispersion management also evokes a certain amount of interest in researchers [27]. In double-ring arrays the possibility to manage the spectral characteristics is primarily connected with the ability to control the coupling constants I,J,K by varying the geometrical parameters of the array, and in this way allows more versatility as compared to single-ring arrays. Nevertheless, no evident peculiarity of energy transfer is connected with such a phenomenon. It is the other feature of spectral curves which is responsible for the observed effects.

As is seen from Figs. 5(a) and 5(b), as r 1 increases all the spectral curves tend to two horizontal asymptotics, which leads to an N -fold degeneracy of the spectrum. Analytically, this can be explained as follows. At r 1 → ∞ (at a constant inter-ring distance) the difference between I,J tends to zero. Simultaneously, the spectrum Eq. ( 15) goes over into

β m = β + I β cos ϕ m ± K 2 β .
Since also I → 0 there appears degeneracy in m. This degeneracy leads to suppression of energy spreading over the ring. However, between the rings the energy exchange remains possible. The main factor that facilitates such inter-ring exchange is the coincidence of modal coefficients cos α m , sin α m in Eq. ( 14), which at I → J tend to 1/ √ 2. The resulting modes are given by symmetric and antisymmetric combinations of single-ring modes [28], which propagate with phase velocities that differ by K/ β. In this way the energy oscillates between the rings almost without spreading over them, which is clearly seen in Fig. 3(c). The darker stripes indicate at the total transfer of energy to the other ring. The spatial scale of such oscillations depends on the inter-ring coupling and has the order of β/K.

In general, energy transfer between the rings is determined by the value of modal coefficients. For example, if one excites in the inner ring its mode X m , the amplitude of the mode given after Eq. ( 14), one can deduce that a complete tunneling of energy into the outer ring takes place whenever tan 2α m is infinite. This takes place at I → J , which explains the almost complete energy transfer for large r 1 . However, there exists the possibility of complete mode tunneling even if the difference between I and J is comparable with the value of K. Since tan 2α m is infinite also at cos 2ϕ m = 0, the X modes with m = N 4 , 3N 4 would possess this property of complete tunneling at arbitrary relations between the coupling constants. The energy stored in such modes is not confined to the inner ring and can be completely transferred to the outer ring, and vice versa. This is possible, however, only if N is a multiple of 4.

Among other phenomena, which take place in a circular arrays of fibers, it is worth studying the effect of the presence of the second ring on such phenomena as discrete diffraction [4] and formation of linear azimuthons [14,29]. The most interesting is the case of substantial energy exchange between the rings since the opposite case is reduced to the evolution in a single-ring array. Numerical simulations show that the presence of the second ring does not alter the picture of diffraction if one excites only one fiber in the inner ring. As is seen in Fig. 6, the picture of energy spreading remains the same save for the presence of dark horizontal stripes, which can be analogously explained as the energy transfer to the outer ring. In the same way the second ring affects the movement of an azimuthon. Figure 7 shows the spreading of energy over the array if one excites a group of adjacent fibers in the inner ring. The initial amplitude is assumed to be equal, whereas the initial phase difference between the neighboring fibers of the excited domain is ϕ = 2πl/N, where l is an integer. In this way the excited domain represents a part of the phase distribution of discrete optical vortices (DVs). Evident overall rightward movement of the excited domain should be treated as an azimuthal evolution of the excited light spot now recognized as a linear azimuthon. The other ring is responsible only for the appearance of dark horizontal stripes.

There is one more interesting effect that accompanies energy spreading in double-ring arrays-the effect of energy focusing to the inner ring. To illustrate this effect consider the case of strong energy exchange between the rings and study the evolution of the total energies the rings decreases; the larger energy becomes stored in the inner ring. This is quite natural since the total energy tunneling into the other ring is possible only in some limiting cases: The energy, which is initially concentrated in the inner ring, only partially tunnels into the outer ring. The reverse is also true, which complies with the equivalency of the rings. However, if one simultaneously excites both rings (q m0 = 12 k=8 δ m,k and p m0 = 12 k=8 δ m,k ), there will take place a final concentration of energy in the fibers of the inner ring [see Fig. 8(b)]. At some z up to 90% of the initial energy is concentrated in the inner ring-an effect not observed upon a separate excitation of the rings. This effect is best demonstrated by the plot of average energy W2 where a pronounced dip is present [Fig. 8(c)].

W 1 = N k=1 |A k | 2 and W 2 = N k=1 |B k | 2 of

IV. ORBITAL ANGULAR MOMENTUM OF SUPERMODES

The interest in the question of OAM of light-the classical topic in electrodynamics-has been renewed since the studies of OAM of Laguerre-Gaussian beams [30]. Since then it has been studied in an almost boundless body of research (see, for example, [31] and references therein). Recently the attention of theoreticians has been attracted by the question of OAM of beam [START_REF] Ya | [END_REF] and fiber arrays [14,33]. In this section we generalize the approach suggested in [14] to the case of OAM of a double-ring array.

Following [START_REF] Berry | Proc. SPIE[END_REF], let us calculate the ratio of time-averaged OAM linear density M z to averaged linear energy density W with the help of the relation

M z W = 1 ω ψ| lz |ψ ψ|ψ , ( 26
)
where ω is radiation's frequency, lz = -i ∂ ∂ϕ and the scalar product is defined as a|b = S a * bdS, S being the total cross section. In using the scalar variant of a more general relation, in which lz is substituted by lz + σz [START_REF] Berry | Proc. SPIE[END_REF], we assume that the field of the array has linear polarization and thus can be described by a scalar function ψ.

The field of the supermode Z (+) m in the cross section can be represented in the global polar coordinates (see Fig. 1 for designations) as (

Z (+) m = E exp -r 2
) 27
Here we have taken into consideration that the distance r jn between a point (r,ϕ) and the centre of the j th fiber in the nth ring is r 2 jn = r 2 + r 2 n -2rr n cos(ϕϕ j ) [14]. Using the technique described in [14], one can obtain after lengthy but straightforward calculations,

Z (+) m | lz |Z (+) m = πN 2 E 2 N-1 j =0 sin ϕ j sin(mϕ j ) × cos 2 (α m )r 2 1 exp - r 2 1 r 2 0 sin 2 ϕ j 2 + sin 2 (α m )r 2 2 exp - r 2 2 r 2 0 sin 2 ϕ j 2 + sin(2α m )r 1 r 2 exp - (r 1 -r 2 ) 2 4r 2 0 × exp - r 1 r 2 r 2 0 sin 2 ϕ j 2 . ( 28
)
For Z (-) m mode one can obtain the analogous expression by making in Eq. ( 28) a formal change: cos α → sin α and sin α →cos α. For linear energy density of the supermode one gets in the same way:

Z (+) m |Z (+) m = πE 2 r 2 0 N N-1 j =0 cos(mϕ j ) cos 2 (α m ) exp - r 2 1 r 2 0 sin 2 ϕ j 2 + sin 2 (α m ) exp - r 2 2 r 2 0 sin 2 ϕ j 2 + sin(2α m ) × exp - (r 1 -r 2 ) 2 4r 2 0 exp - r 1 r 2 r 2 0 sin 2 ϕ j 2 , ( 29
)
and, analogously, for Z (-) m supermode. Figure 9 shows the dependencies of specific OAM defined by Eq. ( 26) for Z (+) m supermode as functions of reduced distance d/r 0 between the rings [Fig. 9(a)] and the number m of the supermode [Fig. 9(b)]. One should notice that starting from some number the OAM of the mode decreases whereas the topological charge of the DV still increases. This illustrates a general statement that the OAM and the topological charge of the optical field are not rigidly correlated. Comparing Eq. ( 28) with the corresponding result of Ref. [14] [Eq. [START_REF] Bjarklev | Photonic Crystal Fibres[END_REF]] one can identify the first and the second terms in curly brackets as contributions of inner and outer rings of the array [see also Eq. ( 13)], correspondingly. The third term should be identified as some "interference" OAM.

The last comment we would like to make concerns conservation of angular momentum. As has been proved in [14], the OAM of a mode superposition in a circular array is z independent. In a double-ring array, however, this is not so. Though this fact can be rigorously proved for any superposition of modes using the methods suggested in [14], we will demonstrate it on the example of the array with a symmetric exchange I = J = K. In this case the modes equation ( 14) assume the simplest form:

Z (±) m = 1 √ 2 (X m ± Y m
) with propagation constants β ± . If at z = 0 one excites only the mode X m of the inner ring, the field in the array will then evolve as

∝ X m cos(δβz) + iY m sin(δβz), (30)
where δβ = (β +β -)/2. As is evident, at certain z the power in the array will be concentrated on the outer ring; that is, the amplitudes q i of the fibers in the inner ring will be zero. In this case one can calculate the OAM using the expression for a circular array [14]:

M z ∝ r 2 i N-1 j =0 sin ϕ j sin(mϕ j ) exp - r 2 i r 2 0 sin 2 ϕ j 2 , ( 31
)
r i being the radius of the ring. As follows from this relation, M z at z = 0 and z n = πn/δβ do not coincide since for these cross sections the radii r i of single-ring arrays, which the energy is concentrated in, are r 1 and r 2 , respectively. Therefore, the OAM in a double-ring array is not conserved. One should emphasize, though, that this nonconservation does not relate to a time domain and has no direct connection to the Noether theorem. However, there is a subtlety to this otherwise evident result. As follows from Eq. ( 29) (see also Eq. ( 21) of Ref. [14]), the energy of a circular array also depends on its radius:

W ∝ N-1 j =0 cos(mϕ j ) exp - r 2 i r 2 0 sin 2 ϕ j 2 . ( 32
)
Though in zero approximation in r 2 i /r 2 0 its value is radius independent, the first-order corrections yield a slow dependence of energy on the array's radius. If this causes no paradoxes in the case of single-ring arrays, for multiring arrays to avoid this one should retain only the major term in the sum Eq. ( 32) to ensure conservation of energy. The nature of this paradox, probably, can be attributed to certain inconsistency of next-nearest-neighbor approximation used for obtaining the modes of circular arrays. Whereas Eq. ( 32) implies delocalization of energy over the whole array, this approximation disregards such effects in principle. To obtain self-consistent results one should retain in Eq. ( 32) only the first term, which implies localization of energy on a site. In other words, to avoid such controversies one should obtain the modes of an array through the procedure that takes account of all the couplings, which are present in the system. One possible way to elaborate the model could be using the expression for long-range coupling suggested in [START_REF] Martinez | [END_REF]. Otherwise, using incorrect expressions for modes may lead to unphysical results of the above mentioned type. It should be emphasized, though, that the source of mistakes is not in the scheme of calculation of energy and OAM, but in imperfect expressions for modes.

Concerning nonconservation of OAM, it is difficult to make a conclusion on the validity of this result. Since its value changes insignificantly we would rather attribute it to imperfection of the nearest-neighbor approximation.

V. LATTICE CURRENT AND "ORBITAL ANGULAR MOMENTUM"

The spreading of light over an optical system can be described by a number of appropriate optical currents [36]. In waveguide arrays evolution of delocalized energy can be effectively described through the same mathematical constructions as is done in the empty space, which has been successfully demonstrated by our previous paper [14]. It is desirable, however, to find out alternative ways of energy transfer description, which would allow for the fact that the energy of the array bears the features of pointlike localization being concentrated at the vertices of the fiber lattice. In such a situation one can search for lattice adaptations of common mathematical objects used for the description of standard continuous fields.

The simplest choice here seems to be adaptation of the notion of a continuous current j of a scalar complexvalued function ψ introduced in quantum mechanics as j ∝ i 2 ψ * ∇ψ + c.c. [START_REF] Davydov | Quantum Mechanics[END_REF]. A lattice one-dimensional analog of this current can be constructed using the following correspondence: ψ → ψ k and ∇ψ → ψ kψ k-1 , where k is the number of the vertex. These replacements lead to the following expression of a discrete current j k :

j k = i 2 ψ * k ψ k-1 + c.c. ( 33
)
Here ψ k is a complex amplitude at the kth vertex and the factor i/2 stands for certain chosen normalization. However simple, Eq. ( 33) can be very useful for understanding the details of energy transfer in discrete systems. First of all, description in terms of the discrete current is by orders much more technically simple than the one that employs the standard Poynting vector formalism. It can be effectively applied for understanding intricate questions concerned with energy transfer in fiber arrays. As an example let us consider its application to the question of DVs [15,37] in a circular fiber array. As has been shown [14], a mode

X m = 1 √ N N-1
k=0 exp(ikϕ m )q k of a circular fiber array represents a DV with topological charge m if N/2 > m, and m-N if N/2 < m. In the particular case N = 2m the field describes the set of edge dislocations and no DV is formed. The amplitudes for the mth mode read as ψ k = exp(ikϕ m ). Using Eq. ( 33) one can readily arrive at the result that for such a mode the discrete current is constant and is equal to

j (m) = sin ϕ m . ( 34
)
This result indicates the peculiarity of DVs connected with impossibility to create a DV with the charge more than N/2 at a given N . Indeed, if m surpasses [N/2] the value of j (m) changes its sign. This fact corresponds to the appearance of a DV with a negative topological charge due to a "stroboscopic" effect [14]. For even N at m = N/2 the discrete current j (m) is zero, which correlates with the appearance of degenerate edge dislocations of the phase front. In this way the introduced discrete current may be useful as an additional characteristic of energy transfer in the array. Moreover, it has a direct connection with the genuine OAM of a supermode. Indeed, as can be shown from Eq. ( 31) the OAM of the mth supermode of a single-ring array in the main approximation is equal to

M z ∝ r 2 1 sin 2π N sin ϕ m exp - r 2 1 r 2 0 sin 2 π N ∝ sin ϕ m . ( 35
)
Comparing Eqs. [START_REF] Berry | Proc. SPIE[END_REF] and ( 35) one can deduce a formal proportionality between OAM of the supermode and the discrete current.

It turns out that the thus introduced lattice current makes it possible to construct an integral of motion for an array of coupled fibers. Heuristic arguments for suggesting the form of this conserved quantity are based on a formal analogy between lattice current and the transverse component of the Poynting vector. Since for a single-ring array the product of the lattice current j k and the ring's radius r is analogous to OAM density (its z component), the total "lattice OAM" M can be obtained through summation of such point densities over the array. Since the radius is the same for each vertex, the resulting "lattice OAM" can be introduced as:

M = k j k . ( 36
)
As can be easily verified, for a single-ring array this sum is z independent for any condition of the array's excitation; that is, such "pseudo-OAM" represents a conserved quantity, or the integral of motion if one invokes the analogy between time t and longitudinal coordinate z in optical fibers.

For a double-ring array this formal analogy between OAM and M no longer holds since the conserved quantity is M = k j (1) k + j (2) k , (37) where j (i) k is the lattice current in the ith ring. Generally speaking, no information about the array's geometry is encoded in Eqs. (36) and (37). As will be shown further, the conserved quantity M rather relays to some set of differential equations, which describes the propagation of coupled fiber modes, than to any specific configuration of the array. Through Eqs. ( 22)-(24) using Eq. ( 33) one can show that M is determined by the boundary complex amplitudes in the following way:

M = N-1 l,m,m =0 sin ϕ l exp {il [ϕ m -ϕ m ]}(q * m0 q m 0 + r * m0 r m 0 ). (38)
This expression can be brought to a more lucid form, using mode decomposition coefficients introduced in Eq. ( 25) :

M = N-1 l=0 sin ϕ l [|V l | 2 + |U l | 2 ]. (39)
Comparing this result with Eq. ( 34) one can conclude that the sum of partial discrete currents of DVs induced in the rings is conserved. This conservation law resembles the one revealed for OAM of superposition of supermodes in a single-ring array [14]. The existence of such a specific invariant is not occasional and is closely connected with the basic properties of this system. To understand the reason for the appearance of this property it is useful to notice that the system in study can be equivalently described by the coupled-modes approach. Within its frameworks the field amplitudes Āk , Bk are the slowly varying envelopes and satisfy the following set of equations: Evidently, the property Eq. (43) takes place not only for the model of next-neighbors interaction but for a much wider class of interaction models. In this way the lattice "angular momentum" is the integral of motion for the system of coupled concentric circular fiber arrays.

It should be emphasized that this novel integral is a less fundamental quantity as, for example, the energy of the system. It can be shown that M does not conserve if the numbers of fibers in the outer and inner rings are different. In this case the property Eq. ( 41) of the matrix that describes the coupling between the rings does not hold. Also M is not the integral of motion in the presence of nonlinearity: Adding the terms proportional to |A k | 2 A k and |B k | 2 B k to the right-hand sides of Eq. (41) would destroy this symmetry. Recently nonconservation of topological charge in circular arrays of nonlinear fibers, which naturally implies nonconservation of the corresponding OAM, has been shown in Ref. [15]. Though there is no direct connection with OAM and the introduced lattice pseudo-OAM, our results somewhat correlate with the effect of DV's sign switching established in Ref. [15].

VI. CONCLUSIONS

In this paper we have studied the energy transfer and orbital angular momentum in a double-ring array of monomode fibers. On the basis of the obtained rigorous expressions for supermodes and the spectra of their propagation constants we have investigated the energy transfer between the rings of the array. We have determined the geometrical parameters of the array, at which the energy is mostly confined within the layers. We also showed how the developed method can be applied to the case of multiring arrays. We have calculated the orbital angular momentum carried by supermodes and showed that its value can be controlled by the changing array's geometrical parameters. We have introduced the discrete lattice current and shown that an integral of motion for circular arrays-the lattice pseudo-OAM-is connected with this quantity. We have also made the connection of the total discrete current with orbital angular momentum and shown that this integral can be useful while studying angular momentum of discrete optical vortices.

) 1 FIG. 1 .

 11 FIG. 1. (Color online) Geometry of a double-ring fiber array and the scheme of fibers' numeration.

FIG. 2 .

 2 FIG. 2. Coupling constants vs reduced radius of the first ring r 1 /r 0 . Array's parameters: N = 21, distance between the rings d = 3r 0 ; = 10 -3 , r 0 = 10λ He-Ne , n co = 1.48.

FIG. 3 .

 3 FIG. 3. Evolution with z of the pattern of light's intensity (|A k | 2 and |B k | 2 ) on the fibers of the rings; k is the number of the fiber. In these examples it is supposed that one excites with equal weights the fibers 9-13 of the inner ring (q m0 = 12 k=8 δ m,k ). The radius of the inner ring is (a) r 1 = 8r 0 ; (b) r 1 = 10r 0 ; (c) r 1 = 12r 0 . The other parameters of the array are N = 21, d = 3r 0 , = 10 -3 , r 0 = 10λ He-Ne , n co = 1.48. Panel (c) presents the case of strong energy exchange between the rings.

5 FIG. 5 .

 55 FIG. 5. (Color online) Corrections β (±) k = β (±) k -β to the scalar propagation constant vs reduced radius of the first ring r 1 /r 0 ; N = 21, d = 3r 0 , = 10 -3 , r 0 = 10λ He-Ne , n co = 1.48. For convenience of presentation the spectral curves are given in different figures [(a) and (b)]. Note that with the increase of r 1 /r 0 all the curves tend to horizontal asymptotics, which are the same for both figures (a) and (b).

  Y m in the outer ring can be shown to equal |sin 2α m sin βz 2 |, where β is the difference in propagation constants of supermodes Z (+) m , Z (-) m , composed of the modes X m and Y m . Maximal energy accumulated in the outer ring is determined by |sin 2α m | = |tan 2α m | √ 1+tan 2 2α m . Using the explicit form of tan 2α m

FIG. 6 .

 6 FIG. 6. Pattern of light's intensity (|A k | 2 and |B k | 2 ) in the rings as a function of z upon excitation of the 11th fiber of the inner ring (q m0 = δ m,11 );N = 21, r 1 = 12r 0 , d = 3r 0 , = 10 -3 , r 0 = 10λ He-Ne , n co = 1.48.

  FIG. 7. Pattern of light's intensity in the rings as a function of z upon equal weights excitation of the fibers 9-13 of the inner ring. The phase difference between the neighboring excited fibers at the input is 2π N so that q m0 = exp[i 2π N (m -8)] 12 k=8 δ m,k . Array's parameters: N = 21, r 1 = 12r 0 , d = 3r 0 , = 10 -3 , r 0 = 10λ He-Ne , n co = 1.48.

FIG. 8 .

 8 FIG. 8. (Color online) Energy focusing in the inner ring. The plots of the total energies W 1 , W 2 (in arbitrary units) of inner and outer rings vs z upon excitation of the inner ring only: q m0 = 12 k=8 δ m,k (a) and simultaneous excitation of the both rings: q m0 = 12 k=8 δ m,k ,p m0 = 12 k=8 δ m,k (b). (c) Plot of averaged over z energy W2 (in a.u.) stored in the second ring vs reduced radius of the first ring for the same excitation as in (b). Note that for (a), W 1 (z = 0) = 5, W 2 (z = 0) = 0; for (b) and (c), W 1 (z = 0) = W 2 (z = 0) = 5.

× 1 2r 2 0+

 12 cos α m exp 2rr 1 cos(ϕϕ j )r 2 sin α m exp 2rr 2 cos(ϕϕ j ) -

FIG. 9 .

 9 FIG. 9. Specific orbital angular momentum (L z = ω Mz W ) of the supermode Z (+) m vs (a) reduced distance d/r 0 between the rings; (b) number m of the supermode (at constant geometrical parameters of the array). Note that starting from some number the angular momentum of the mode decreases whereas the topological charge of the discrete optical vortex still increases. Array's parameters: N = 21, r 1 = 12r 0 , d = 3r 0 , = 10 -3 , r 0 = 10λ He-Ne , n co = 1.48.

FIG. 10

 10 FIG. 10. (Color online) Oscillations of the total lattice current of the inner ring reduced to the total energy of the ring (in arbitrary units). The fibers 9-13 of the inner ring are excited with equal weights. The phase difference between the neighboring excited fibers at the input is 2π N ; q m0 = exp[i 2π N (m -8)] 12 k=8 δ m,k . Array's parameters: N = 21, r 1 = 12r 0 , d = 3r 0 , = 10 -3 , r 0 = 10λ He-Ne , n co = 1.48.

  Āk + I km Ām + K Bk , (41)i d Bk dz = β Bk + J km Bm + K Āk ,where β is the propagation constant for a single fiber, and I km and J km are the exchange matrices that describe coupling within the rings. For the next-neighbors approximation they are proportional to δ k,l+1 + δ k+1,l . Using Eqs. (41) and their complex conjugates one can demonstrate that the lattice "angular momentum" is conserved: matrices possess the property I lm = I l+1,m+1 .(43)

ACKNOWLEDGMENTS

The authors acknowledge the helpful suggestion of the referee to make analytical averaging of the total energy over z.