Massimiliano Guasoni

Matteo Conforti

Costantino De Angelis

Light propagation in nonuniform plasmonic subwavelength waveguide arrays

We study light propagation in nanoscale periodic structures composed of dielectric and metal in the visible range. We demonstrate that diffraction curves of nonuniform waveguide arrays can be tailored by varying the geometric and dielectric features of the waveguides. The results obtained from a proper formulation of Coupled Mode Theory for non uniform arrays are validated through numerical solution of Maxwell equations in frequency domain.

I. INTRODUCTION

The miniaturization of photonic devices for confining and guiding electromagnetic energy down to nanometer scale is one of the biggest challenges for the information technology industries [START_REF] Maier | Plasmonics: Fundamentals and Applications[END_REF]. In the last years, photonic crystals technology allowed to gain one order of magnitude in the miniaturization of components such as waveguides and couplers with respect to conventional (i.e. based on total internal reflection) optics. However when the size of a conventional optical circuit is reduced to the nanoscale, the propagation of light is limited by diffraction. One way to overcome this limit is through surface plasmon polaritons [START_REF] Economu | Surface Plasmons in thin films[END_REF], which are evanescent waves trapped at the interface between a medium with positive real part of dielectric constant and one with negative real part of dielectric constant, such as metals in the visible range. Even though this phenomenon has been known for a long time, in the last years there is a renewed interest in this field, mainly motivated by the wide range of potential applications that sweep from the realization of biologic nanosensors [START_REF] Anker | Biosensing with plasmonic nanosensors[END_REF], to sub-wavelength imaging [START_REF] Fang | Sub-diffraction-limited optical imaging with a silver superlens[END_REF], to the merging of electronic circuits to photonic devices [START_REF] Ozbay | Plasmonics: Merging photonics and electronics at nanoscale dimensions[END_REF].

On the other hand, control of light propagation by means of periodic photonic structures is a fundamental issue that is attracting a lot of interest in the scientific community. In particular, arrays of evanescently coupled waveguides are unique structures that exhibit the peculiar properties of discrete systems. Indeed, light propagation in waveguide arrays is characterized by strong confinement of the field into the individual waveguides and the observable exotic phenomena are due to the weak coupling between the waveguides. As a result, modes of the whole structure can be approximated by a superposition of a discrete set of localized modes, thus light propagation can be considered truly discretized [START_REF] Christodoulides | Discrete self-focusing in nonlinear arrays of coupled waveguides[END_REF], [START_REF] Eisenberg | Discrete spatial optical solitons in waveguide arrays[END_REF]. Recently there is a great research effort in the field of discrete effects in plasmonic structures, and some peculiar outcome of discreteness were reported for metal-dielectric waveguide

The authors are with the Dipartimento di Elettronica per l'Automazione and Consorzio Interuniversitario per le Scienze Fisiche della Materia, Università di Brescia, via Branze 38, Brescia 25123, Italy. arrays, for example Bloch oscillations [START_REF] Lin | Spatial Bloch oscillations of plasmons in nanoscale metal waveguide arrays[END_REF], negative refraction [START_REF] Fan | All-Angle Broadband Negative Refraction of Metal Waveguide Arrays in the Visible Range: Theoretical Analysis and Numerical Demonstration[END_REF] and diffraction management [START_REF] Conforti | Subwavelength diffraction management[END_REF].

Nonuniform waveguide arrays have received increasing attention, since a more complex engineering of the periodic structure can provide further degrees of freedom. In this context, the first reported example concerns the usage of zigzag waveguide arrays (i.e., the cascade of arrays characterized by alternating tilt angles) in order to obtain diffraction management. Binary arrays composed of waveguides with alternating widths have been thoroughly studied [START_REF] Sukhorukov | Discrete gap solitons in modulated waveguide arrays[END_REF], [START_REF] Longhi | Multiband diffraction and refraction control in binary arrays of periodically curved waveguides[END_REF] since they exhibit interesting features, such as double refraction, due to their intrinsic two-band nature.

In the present work we study the behavior of nonuniform metal-dielectric waveguide arrays composed of waveguides with different dielectric cores, that determine strong variations of the coupling coefficients. Coupled mode theory (CMT) is extended in order to deal with plasmonic modes and varying coupling coefficients, and further improvements are proposed in order to take into account the different widths of the two array bands. Moreover, we demonstrate that the ability to control the magnitude of the coupling between the waveguides opens the way to the design of binary waveguide arrays with unusual properties, such as almost flat diffraction curves, that are required, for example, to achieve self collimation [START_REF] Kosaka | Self collimating phenomena in photonic crystals[END_REF]. Finite Element solution of Maxwell equations in nonuniform plasmonic waveguide arrays are reported in order to assess the validity of the analytical treatment.

The paper is organized as follows. After this Introduction, in sec. II we introduce the nonuniform arrays and report some results of the standard CMT. In sec. III we develop an extension of coupled mode theory for non-uniform waveguide arrays that exploits the unsymmetric coupler as the fundamental cell; this extension is motivated by the unacceptable low accuracy of the standard CMT for nonuniform plasmonic waveguide arrays. In sec. IV we calculate Bloch modes of the array and the energy distribution between different bands. Due to the high energy concentration characteristic of plasmonic waveguide we obtain simple approximate (yet accurate) expressions for the energy distribution between bands. In sec. V we exploit the theoretical results of previous sections to engineer diffraction in nonuniform waveguide arrays, reporting as a relevant example the design of arrays characterized by a flat diffraction curve. We end with the conclusion in sec. V.

II. NON UNIFORM PLASMONIC ARRAYS

We consider a one dimensional (1D) array formed by the alternation of two dielectric layers (cores) divided by a metallic Fig. 1. The plasmonic array studied in this paper. Two cores with dielectric constant ε 1 and ε 2 (white and black respectively) are alternated and divided by a metal layer (grey) whose dielectric constant is εm. Width of cores is 2b, width of metal layers is 2a. The fundamental period Λ is 4a + 4b. Dashed rectangles surround 3 consecutive basic cells of the array, that according to the CMT developed in this paper are plasmonic couplers. one (cladding), as sketched in Fig. 1. In a CMT approach we can consider, as basic cell, the single isolated waveguide formed by a dielectric layer surrounded by metal (Fig. 2). Let's call ε 1 , ε 2 and ε m the relative dielectric constants of dielectrics and metal in the array, and 2a and 2b the widths of metallic and dielectric layers, respectively. Metal dielectric constant ε m is calculated by means of Drude's model, i.e.

ε m = 1 -ω 2 p /(ω 2 -iγω),
where ω p and γ are the plasma and collision frequency of the metal, respectively. Supposing that ε 2 > ε 1 (from here and for the rest of the article), we work at frequencies such that ω < ω p / √ 1 + ε 2 . Under this condition, the fundamental mode of the single waveguides is TM-even and can be expressed as a sum of decaying exponentials both in dielectric and metal; if 2b is sufficiently small the guides are monomodal, that is the necessary condition to obtain a bimodal array. Indeed the array supports two different modes since the alternated guides have different propagation constants β 1 and β 2 . According to CMT model, diffraction functions of the two array modes are [START_REF] Sukhorukov | Discrete gap solitons in modulated waveguide arrays[END_REF]:

k z(1,2) (k) = β 1 + β 2 2 ± β 1 -β 2 2 2 + 4C 2 cos k 2 ( 1
)
where k is the tilt angle normalized with respect to the period of the array (-π < k < π) and C is the coupling coefficient between modes of adjacent waveguides, whose value is [START_REF] Hardy | Coupled mode theory of parallel waveguides[END_REF]:

C = ωε 0 +∞ -∞ [ε arr (x) -ε g1 (x)] E x1 (x)E * x2 (x) dx + +ωε 0 +∞ -∞ ε g2 (x) [ε arr (x) -ε g1 (x)] ε arr (x) E z1 (x)E * z2 (x) dx (2)
In Eq.( 2) ε 0 is the dielectric constant of vacuum, E x1 , E x2 , E z1 and E z2 are the modal electric fields of the two isolated waveguides, ε g1 (x) and ε g2 (x) are the dielectric profiles of the two unperturbed waveguides and ε arr (x) is the dielectric profile of the whole array.

Looking at Eq.( 1) we see that by controlling ∆β = (β 1β 2 )/2 and C we can obtain diffraction curves with very different amplitudes, where for amplitude we mean k z(1,2) (π)k z(1,2) (0). For example when ∆β >> C, diffraction curves become nearly flat. The propagation constant β of a single guide with relative dielectric constant ε and width 2b at pulsations ω that satisfy the inequality ε/|ε m (ω)| < 0.5 is well approximated by Eq.(3) (see Appendix):

β ≈ √ ε w c 0 2 + w p |ε m |c 0 b (3)
where c 0 is light speed in vacuum. It's easy to note that β is very sensitive to ε, so that a sensible difference ∆β arises (β 2 /β 1 = ε 2 /ε 1 ) by coupling waveguides with cores with different dielectric constant ε 1 and ε 2 . In contrast, the dependence of β on the width b is usually weaker than the dependence on ε. For example, if we use silver (neglecting losses: w p = 13.6884e15 rad/s, γ = 0 rad/s), at any pulsation smaller than ω p / √ 1 + ε, β is almost independent on b, for b > 40nm : so ∆β will be quite small when coupling waveguides with the same dielectric but different width.

III. COUPLED MODE THEORY FOR PLASMONIC ARRAYS

It's well known that CMT gives good results when the modes of two coupled waveguides are, with good approximation, a linear combination of the modes of the single unperturbed waveguides. From Eq.(1) we derive that the second array mode has a diffraction curve k z2 (k) whose values can be smaller than ωε 2 /c 0 , implying that field becomes sinusoidal in the cores of the waveguides with dielectric ε 2 . In these conditions a linear combination of single waveguides modes (that are combination of exponentials in the cores) can not well approximate the real array mode, making the value of k z2 (k) calculated from Eq.(1) unreliable.

To overcome this problem we consider a refined CMT model where the single cell of the array consists of two adjacent waveguides (i.e. a waveguide coupler) with dielectrics ε 1 and ε 2 , respectively (Fig. 2). We expect this alternative formulation of CMT to work better, because it takes into account the effective fundamental basic cell of the non uniform array (that is the coupler). We will consider linear combination of the coupler modes, that can be sinusoidal in waveguides with dielectric ε 2 , approximating better the second array mode.

Let's consider bimodal couplers (any coupler has at least two TM-modes), and let's use a CMT model in which array modes are linear combinations of the two modes of all the couplers in the array, so that the total transverse electric field E x (x, z) can be written as:

E x (x, z) = n A n,a (z)E x(n,a) (x) + A n,b (z)E x(n,b) (x) + E res (x) (4)
where E x(n,a) (x) and E x(n,b) (x) are the two modes (denoted with a and b) of the n-th coupler in the array, A n,a (z) and A n,b (z) are their amplitudes and E res (x) is the residual field. Supposing that the residual field is negligible and following the treatment reported in [START_REF] Hardy | Coupled mode theory of parallel waveguides[END_REF], we find that amplitudes are connected by the relations below: [START_REF] Christodoulides | Discrete self-focusing in nonlinear arrays of coupled waveguides[END_REF] In Eq.( 5) and Eq.( 6) A n,i (i = a, b) is the derivative respect to z of the amplitude A n,i . Any term R i1,jl is the correlation between mode i in the coupler n and mode j in coupler m = n + l -1 (l = 0 with reference to the coupler at the left of the coupler n and l = 1 with reference to the coupler n and l = 2 with reference to the coupler at its right, see Fig. 1) , that is:

A n,a + R a1,a2 A n+1,a + R a1,a0 A n-1,a + +R a1,b2 A n+1,b + R a1,b0 A n-1,b = i(β a + k a1,a1 )A n,a + +iC a1,b1 A n,b + iC a1,a2 A n+1,a + +iC a1,a0 A n-1,a + iC a1,b2 A n+1,b + iC a1,b0 A n-1,b (5) A n,b + R b1,b2 A n+1,b + R b1,b0 A n-1,b + +R b1,a2 A n+1,a + R b1,a0 A n-1,a = i(β b + k b1,b1 )A n,b + +iC b1,a1 A n,a + iC b1,b2 A n+1,b + +iC b1,b0 A n-1,b + iC b1,a2 A n+1,a + iC b1,a0 A n-1,a
R i1,jl = +∞ -∞ E x(n,i) (x)H * y(m,j) (x)dx (7)
under the normalization condition R i1,i1 = 1, while R i1,j1 = 0 (i = j) because different modes in the same coupler are orthogonal. β a and β b are the propagation constants of modes a and b in any basic unperturbed coupler, while any term C i1,jl is the coupling coefficient between mode i in the coupler n and mode j in coupler m = n + l -1 . This coefficient can be written as:

C i1,jl = β i R i1,jl + k i1,jl , (8)
where the term k i1,jl is:

k i1,jl = ωε 0 +∞ -∞ [ε arr (x) -ε cn (x)] E x(n,i) (x)E * x(m,j) (x)+ + ε cm (x) [ε arr (x) -ε cn (x)] ε arr (x) E z(n,i) (x)E * z(m,j) (x) dx, ( 9
)
where ε cn (x) and ε cm (x) are the transverse dielectric profiles of adjacent and unperturbed couplers n and m, while ε arr (x) is the transverse profile of the whole array. Using as solutions A n,a = Ae ikn+ikzz and A n,b = Be ikn+ikzz (-π < k < π) we can rewrite Eq.( 5) and Eq.( 6) in matrix form:

R A B k z = C A B ( 10
)
Where:

R = P 11 P 12 P 21 P 22 , C = Q 11 Q 12 Q 21 Q 22 , ( 11
)
with

Q 11 = β a + k a1,a1 + C a1,a2 e ik + C a1,a0 e -ik , Q 12 = C a1,b1 + C a1,b2 e ik + C a1,b0 e -ik , Q 21 = C b1,a1 + C b1,a2 e ik + C b1,a0 e -ik , Q 22 = β b + k b1,b1 + C b1,b2 e ik + C b1,b0 e -ik , P 11 = 1 + R a1,a2 e ik + R a1,a0 e -ik , P 12 = R a1,b2 e ik + R a1,b0 e -ik , P 21 = R b1,a2 e ik + R b1,a0 e -ik , P 22 = 1 + R b1,b2 e ik + R b1,b0 e -ik .
Diffraction curves k z (k) are the eigenvalues of R-1 C; being coefficients R ix,jy << 1 in Eq.( 11), we can well approximate R-1 ≈ (2I -R). In the product (2I -R) C we can neglect all terms different from C i1,jl , β a , β b , β a R i1,jl and β b R i1,jl because they are much smaller. In this way, matrix (2I -R) C becomes:

(2I -R) C = M 11 M 12 M 21 M 22 (12)
where

M 11 = β a + k a1,a1 + k a1,a2 e ik + k a1,a0 e -ik M 12 = k a1,b1 + (k a1,b2 + 2∆βR a1,b2 )e ik + (k a1,b0 + 2∆βR a1,b0 )e -ik M 21 = k b1,a1 + (k b1,a2 -2∆βR b1,a2 )e ik + (k b1,a0 -2∆βR b1,a0 )e -ik M 22 = β b + k b1,b1 + k b1,b2 e ik + k b1,b0 e -ik and ∆β = (β a -β b )/2.
If we neglect the residual field E res (x), the system conserves the energy, implying that R-1 C must be hermitian. This property imposes the following equalities:

k a1,a2 = k a1,a0 = C a , k b1,b2 = k b1,b0 = C b , k a1,b2 + 2∆βR a1,b2 = k b1,a0 -2∆βR b1,a0 = C ab , k a1,b0 + 2∆βR a1,b0 = k b1,a2 -2∆βR b1,a2 = C ba , k a1,b1 = k b1,a1 = C.
Moreover, calling β a + k a1,a1 = βa and β b + k b1,b1 = βb , we can rewrite [START_REF] Longhi | Multiband diffraction and refraction control in binary arrays of periodically curved waveguides[END_REF] in this way:

(2I -R) C = βa + 2C a cos(k) C + C ab e ik + C ba e -ik C + C ba e ik + C ab e -ik βb + 2C b cos(k) (13)
Diffraction curves are the eigenvalues of ( 13):

k z(1,2) = βa + βb 2 + (C a + C b ) cos(k) ± βa -βb 2 + (C a -C b ) cos(k) 2 + |C + C ab e ik + C ba e -ik | 2 (14)
Even when difference (ε 2 -ε 1 ) is small, the two modes of any coupler concentrate energy in a different way: the first confine the most part of energy in the core with dielectric constant ε 2 , while the second in the core with dielectric ε 1 , so that the term ( βaβb ) increases and coefficients C, C ab and C ba becomes negligible respect to [( βaβb )/2 + (C a -C b ) cos(k)] 2 , implying:

k z1 = βa + 2C a cos(k) k z2 = βb + 2C b cos(k) (15)
That is a very simple solution where the two diffraction curves depend only on the coupling between same modes in two adjacent couplers. Their amplitudes are respectively 4C a and 4C b , and they can take very different values, in contrast to what results from Eq.( 1). This result can be found also following [START_REF] Guasoni | Peculiar properties of photonic crystal binary waveguide arrays[END_REF], where interaction with second neighboring guide is taken in account.

IV. MODES AND ENERGY CONCENTRATION IN PLASMONIC

ARRAYS

Now we want to show that the coupling coefficients depend quite exclusively on how field of the unperturbed basic coupler concentrates at the metal-dielectric interfaces. This fact allow us to deeply understand the link between non uniformity of the array (ε 2 = ε 1 ) and its diffraction curves. Let's call g 1 , g 2 and g m the regions with dielectric ε 1 , ε 2 and ε m in the coupler (Fig. 2), and let's start with observing that the first coupler mode is very similar to the mode of the single isolated guide with core ε 2 , excluding the region g 1 . We can then approximate the first coupler mode in all regions except g 1 with the mode of the single guide with core ε 2 . An analogous argument holds true for the second coupler mode.

Coupling coefficient C a = k a1,a0 can be very well approximated considering only the fields overlapping in region g 2 of coupler (n -1) ( see Fig. 3) : 3. The array is shown together with transverse electric field E x of modes a of couplers n-1 and n (dotted and solid lines respectively), that are E x(n-1,a) and E x(n,a) according to the explanation of section IV. Regions g 1 and g 2 of each coupler are shown ((g 1 , n -1), (g 2 , n -1) for coupler n -1; (g 1 , n), (g 2 , n) for coupler n). In the region (g 1 , n -1) the modes E x(n-1,a) and E x(n,a) are much smaller than in region (g 2 , n -1), so that their overlapping in (g 1 , n -1) can be neglected in the calculation of k a1,a0 ( Eq.( 16)). Large discontinuities of the electric field are due to the big dielectric constat change between metal and dielectric. where we distinguish the overlapping between the x and z components of electric fields (C a and C a , respectively). Subscript (g 2 , n-1) refers to the region g 2 of coupler (n-1). Let's consider now an array with dielectrics ε 1 and ε 2 , and let's study non uniformity effects by varying ε 2 keeping ε 1 fixed. In (g 2 , n -1) the fields E x(n,a) (x, ε 2 ) and E z(n,a) (x, ε 2 ) behave nearly such as the mode of the single guide with core ε 2 , that decays in metal with constant T m (see Appendix). Being T m practically independent on ε 2 , also the field E x(n,a) (x, ε 2 ) is quite independent on ε 2 in (g 2 , n-1) except for a constant, so that we can write

C a ≈ ωε 0 (ε 2 -ε m ) g2,n-1 E x(n,a) (x)E x(n-1,a) (x) * C a + ωε 0 (ε 2 -ε m ) g2,n-1 E z(n,a) (x)E z(n-1,a) (x) * C a (16)
x[nm] E x (g 1 ,n-1) (g 2 ,n-1) (g 1 ,n) (g 2 ,n) ε 1 ε m ε 2 ε m ε 1 ε 2 ε m ε m ε m Fig.
E x(n,a) (x, ε 2 ) ≈ E x(n,a) (x, ε 1 )h e
x(L,a) (ε 2 ), where h e x(L,a) (ε 2 ) is the value of the transverse electric field of the coupler mode a at the external left interface and, by definition, h e x(L,a) (ε 1 ) = 1 (see Fig. 4). For the same fact,

E z(n,a) (x, ε 2 ) ≈ E z(n,a) (x, ε 1 )h e z(L,a) (ε 2 )
, where h e z(L,a) (ε 2 ) is the value of the longitudinal electric field of coupler mode a at the external left interface. Similar arguments holds true for fields E x(n-1,a) (x, ε 2 ) and E z(n-1,a) (x, ε 2 ): in (g 2 , n -1) they are sum of evanescent waves decaying with constant T d (see Appendix) that doesn't vary too much as function of ε 2 (when the array is bimodal); so using

E x(n-1,a) (x, ε 2 ) ≈ E x(n-1,a) (x, ε 1 )h i x(R,a) (ε 2 ) and E z(n-1,a) (x, ε 2 ) ≈ E z(n-1,a) (x, ε 1 )h i z(R,a) (ε 2 )
we approximate very well Eq.( 16), being h i x(R,a) (ε 2 ) and h i z(R,a) (ε 2 ) the values of the transverse and longitudinal electric field of the coupler mode a at the internal right interface. On the basis of the arguments reported above, we can easily deduct the following relations:

C a (ε 2 ) C a (ε 1 ) ≈ ε m -ε 2 ε m -ε 1 h e x(L,a) (ε 2 )h i x(R,a) (ε 2 ) C a (ε 2 ) C a (ε 1 ) ≈ ε m -ε 2 ε m -ε 1 h e z(L,a) (ε 2 )h i z(R,a) (ε 2 ) (17)
Reasoning in the same way for C b = k (b1,b2) we can write:

C b (ε 2 ) C b (ε 1 ) ≈ h i x(L,b) (ε 2 )h e x(R,b) (ε 2 ) C b (ε 2 ) C b (ε 1 ) ≈ h i z(L,b) (ε 2 )h e z(R,b) (ε 2 ) (18)
Eq.( 17) and Eq.( 18) show how variations of couplings as function of non uniformity are related only on fields concentration at the interfaces of the unperturbed basic coupler. We note that C a (ε 2 ) and C a (ε 2 ) are always decreasing function of ε 2 , while C b (ε 2 ) and C b (ε 2 ) exhibit a minimum when the condition k z2 (k) ≈ ωε 2 /c 0 is satisfied, corresponding to the threshold between a decaying or oscillating field in the core of the waveguide g 2 . This fact is not intuitive, because we could expect that the more the difference (ε 2 -ε 1 ) is increased the more the modes of adjacent couplers are decoupled.

It is worth showing now how to calculate the Bloch modes of the array from CMT. From Eq.( 4) we see that the transverse electric field propagating in the array (neglecting E res (x)) can be written as:

E x (x, z) = n [A 1 (k)e ikz1z + A 2 (k)e ikz2z ]e ikn E x(n,a) (x) + n [B 1 (k)e ik z1 z + B 2 (k)e ik z2 z ]e ikn E x(n,b) (x) = = n B 1 (k)[r 1 (k)E x(n,a) (x) + E x(n,b) (x)]e ikn e ik z1 z + n B 2 (k)[r 2 (k)E x(n,a) (x) + E x(n,b) (x)]e ikn e ikz2z (19)
where the couples (A 1 (k), B 1 (k)) and (A 2 (k), B 2 (k)) are the eigenvectors of Eq.( 13) corresponding to the eigenvalues k z1 (k) and k z2 (k) respectively, and r 1,2 (k) = A 1,2 (k)/B 1,2 (k). From Eq.( 19) we deduct that the transverse electric field of two Bloch modes m E(1,k) (x) and m E(2,k) (x) are:

m E(1,k) (x) = n [r 1 (k)E x(n,a) (x) + E x(n,b) (x)]e ikn m E(2,k) (x) = n [r 2 (k)E x(n,a) (x) + E x(n,b) (x)]e ikn (20)
Mode m E(1,k) (x) can be written as (see [START_REF] Guasoni | Peculiar properties of photonic crystal binary waveguide arrays[END_REF]):

m E(1,k) (x) = = r 1 (k)E x(0,a) (x) + E x(0,b) (x) * δ Λ (x)e ikx/Λ = = u E(1,k) (x) * δ Λ (x) e ikx/Λ , ( 21
)
where

u E(1,k) (x) = r 1 (k)E x(0,a) (x) + E x(0,b) (x) e -ikx/Λ , ( 22
) and E x(0,a) (x) and E x(0,b) (x) are the two modes of the first unperturbed basic coupler in the array (modes of other couplers are translation of these two), Λ is the array period, δ Λ (x) is the Dirac comb that makes u E(1,k) (x) periodic and * is the convolution operator.

In the same way we can deduct that mode m E(2,k) (x) is:

m E(2,k) (x) = = r 2 (k)E x(0,a) (x) + E x(0,b) (x) * δ Λ (x)e ikx/Λ = = u E(2,k) (x) * δ Λ (x) e ikx/Λ , ( 23
)
where

u E(2,k) (x) = r 2 (k)E x(0,a) (x) + E x(0,b) (x) e -ikx/Λ , ( 24
) The corresponding transverse magnetic fields of Bloch modes are:

m H(i,k) (x) = ωε 0 ε arr (x) k zi (k) m E(i,k) (x) i = 1, 2 (25)
Let's take now a plane wave E x (x) = e ikxx at the array input, then it can be expressed (see [START_REF] Guasoni | Peculiar properties of photonic crystal binary waveguide arrays[END_REF]) as linear combination of array modes:

E x (x) = a 1 m E(1,k) (x) + a 2 m E(2,k) (x) + E res (x) (26) H y (x) = a 1 m H(1,k) (x) + a 2 m H(2,k) (x) + H res (x) (27)
where E res (x) and H res (x) refer to the input component projected upon radiative modes: we are supposing that only two modes propagate in the array, that are m E(1,k) (x) and m E(2,k) (x) (m H(1,k) (x) and m H(2,k) (x) are the correspondent magnetic field modes), with k x = (k + 2πm)/Λ (m integer) [START_REF] Guasoni | Peculiar properties of photonic crystal binary waveguide arrays[END_REF]. We consider (-π/Λ < k x < π/Λ) in order to achieve a good coupling over the propagating array modes. The total input power density is:

+∞ -∞ E x (x)H * y (x) 2 dx = = En 1 |a 1 | 2 +∞ -∞ m E(1,k) (x)m * H(1,k) (x) 2 dx + + En 2 |a 2 | 2 +∞ -∞ m E(2,k) (x)m * H(2,k) (x) 2 dx + + +∞ -∞ E res (x)H * res (x) 2 dx ( 28
)
where En 1 and En 2 represent power densities carried by the first and second array mode respectively. Coefficients a 1 and a 2 can be deducted by the orthogonality relation between modes [START_REF] Hardy | Coupled mode theory of parallel waveguides[END_REF]:

a i = +∞ -∞ E x (x)m * H(i,k) (x)dx +∞ -∞ m E(i,k) (x)m * H(i,k) (x) i = 1, 2 (29)
From Eq.( 25) we see that:

m H(1,k) (x) = ωε 0 ε arr (x) k z1 u E(1,k) (x) * δ Λ (x) e ikx/Λ ≈ ≈ ωε 0 ε cou (x) β a u E(1,k) (x) * δ Λ (x) e ikx/Λ = = u H(1,k) (x) * δ Λ (x) e ikx/Λ ( 30
)
where

u H(1,k) (x) = r 1 (k)H y(0,a) (x) + H y(0,b) (x) e -ikx/Λ , ( 31
) and ε cou (x) is the relative dielectric profile of the unperturbed basic coupler. The approximation done is justified by the fact that k z1 (k) ≈ β a and the energy of the overlapping fields between two adjacent coupler modes is much smaller then the energy of the modes themselves. With the same treatment for m H(2,k) (x) we deduce that:

m H(2,k) (x) = u H(2,k) (x) * δ Λ (x) e ikx/Λ (32)
In Eq.( 30) and Eq.(32) H y(0,a) (x) and H y(0,b) (x) are the transverse magnetic fields of the modes of the unperturbed coupler. The equations above allow us to deduce the power carried by the two array modes written in Eq.( 28) from the analysis of the unperturbed coupler. Indeed we can write:

a 1 = +∞ -∞ e ik x x m * H(1,k) (x)dx +∞ -∞ m E(1,k) (x)m * H(1,k) (x)dx = = +Λ/2 -Λ/2 u * H(1,k) (x)dx +Λ/2 -Λ/2 u E(1,k) (x)u * H(1,k) (x)dx (33)
where the integration can be done over one period. Using orthogonality condition between E x(0,a) , H y(0,b) and 33) can be written as 2(|r 1 (k)| 2 + 1) under the normalization condition that:

+Λ/2 -Λ/2 E x(0,i) (x)H * y(0,i) (x) 2 dx = 1 i = a, b (34)
In order to calculate numerator in Eq.( 33) we can approximate H y(0,a) and H y(0,b) by considering them constant over the two semiperiods of the coupler. So we set H y(0,a) (x) = h a1 when (-Λ/2 < x < 0) and H y(0,a) (x) = h a2 when (0 < x < +Λ/2). Constants h a1 and h a2 are the mean value of H y(0,a) (x) in the regions g 1 and g 2 of the coupler respectively (see Fig. 5). In the same way we set H y(0,b) (x) = h b1 when (-Λ/2 < x < 0) and H y(0,a) (x) = h b2 when (0 < x < +Λ/2). With these approximations numerator becomes:

+Λ/2 -Λ/2 u * H(1,k) (x)dx = = (-iΛ/k)(r 1 (k)h a1 + h b1 ) * (1 -e -ik/2 ) + (-iΛ/k)(r 1 (k)h a2 + h b2 ) * (e ik/2 -1) (35)
The treatment to calculate a 2 is exactly the same as above, replacing r 1 (k) with r 2 (k). Then we can now estimate the power density ratio En 1 /En 2 :

En 1 En 2 ≈ N 1 N 2 D 2 D 1 , ( 36
)
where

N i = |(r i (k)h a1 + h b1 )(1 -e ik/2 ) + (r i (k)h a2 + h b2 )(e -ik/2 -1)| 2 , and D i = |r i (k)| 2 + 1.
As shown in the next section, this approximation is very good and we can deduce that the two modes carry nearly 50% of power whatever the tilt angle k is. Indeed even with small non uniformity h a2 >> h a1 , h b1 >> h b2 , |r 1 (k)| >> 1 and |r 2 (k)| << 1, because terms of the main diagonal in matrix of Eq.( 13) are much grater than the other two. Then as first approximation we can write En 1 /En 2 ≈ |h a2 | 2 /|h b1 | 2 . Being the first and the second coupler modes concentrated in g 2 and g 1 respectively, we can approximate h a1 = h b2 = 0; imposing normalization condition (Eq.( 34)) we obtain 3)). Ratio ε 2 /ε 1 must be limited in order to maintain bimodality, so we consider usually (1 < ε 2 /ε 1 < 1.5), that justify why the two modes carry nearly the same amount of power independently on k. In order to obtain a good coupling between the input and the array modes, transverse electric and magnetic fields of the coupler modes have to be the more constant as possible in the core where they concentrate the most of their energy. Then it's important to evaluate ratio h c /h i between value h c of the mode at the centre of the core and value h i at the interfaces (see Fig. 4). The characteristics of mode a in g 2 can be deducted by the single guide with dielectric ε 2 , so in case of mode a ratio h c /h i can be approximated as 2/(e

|h a2 | 2 β a /ε 2 = |h b1 | 2 β b /ε 1 , then En 1 En 2 ≈ |h a2 | 2 |h b1 | 2 ≈ ε 2 ε 1 , (37)
T d (ε 2 )b + e -T d (ε 2 )b ), where T d is function of ε 2 (see Appendix).
The same treatment can be done for mode b of the coupler in guide g 1 , replacing T d (ε 2 ) with T d (ε 1 ).

V. DESIGN OF A FLAT-DIFFRACTION ARRAY IN THE VISIBLE BAND

In this section we design a non uniform plasmonic array in order to have flat diffraction curves in all the visible range (450 nm-750 nm) and in order to test all predictions done in the previous sections. Thanks to a flat diffraction curve, is possible to prevent either the beam divergence or diffraction broadening, enabling flexible design of light path in plasmon integrated optics [START_REF] Kosaka | Self collimating phenomena in photonic crystals[END_REF].

Dielectric and metal layers of the array are 120 nm and 40 nm wide respectively; ε 1 = 1.0 (air) and metal layers are silver (ω p = 13.61e15 and γ = 0, because we can neglect losses for propagation distances that we consider). In Fig. 6 the amplitudes of the two diffraction curves related to the two modes of the array at 450 nm and 750 nm are shown as function of ε 2 . We consider only values of ε 2 that allow the array to be bimodal. Difference (k z1 (π) -k z1 (0)) is the amplitude of the first diffraction curve (band 1),(k z2 (π)k z2 (0)) is the amplitude of the second (band 2). As predicted, the amplitude of the first curve always decreases as function of ε 2 , while the second exhibits a minimum of its absolute value. A good choice to obtain a flat second curve in all the visible range is to set ε 2 = 1.5, for which we have nearly the lowest amplitude of band 2 in all the range considered.

The chosen parameters guarantee that h c /h i > 0.8 for both the coupler modes in all the visible range. In Fig. 7 we show the comparison between the amplitude of diffraction curves in the case of a uniform array (ε 1 = 1.0, ε 2 = 1.0 or ε 1 = 1.5, ε 2 = 1.5) and a non uniform array (ε 1 = 1.0, ε 2 = 1.5) as In Fig. 8, in the case of ε 1 = 1.0, ε 2 = 1.5, we report a comparison between the numerically calculated amplitudes and those predicted by CMT using Eq. [START_REF] Guasoni | Peculiar properties of photonic crystal binary waveguide arrays[END_REF]. Relative errors are very small, less then 3% in the worst case. In the case of band 2, numerical and predicted curves are undistinguishable. As expected the amplitudes of both the diffraction curves reduces in modulus by increasing the wavelength. In fact the modulus of the dielectric constant of metal increase with wavelength, enabling a stronger confinement into the waveguides.

In Fig. 9 we report the good agreement between coupling coefficients C a and C b calculated with Eq.( 8) and those approximated by Eq.( 17) and Eq.(18). It is worth noting that the second band coupling constant C b has a non-monotonic behavior, exhibiting a minimum of the absolute value at around ε 2 = 2, in accordance with Fig. 6. In Fig. 10 we show numerically calculated diffraction curves at 450 nm and 750 nm, and those predicted by CMT in Eq.( 15): also in this case, the agreement is very good. Interestingly enough the curvature changes sightly from 450nm to 750nm case, indicating the broadband feature of the diffraction engineering.

In Fig. 11 we report a comparison between the numerically calculated power density ratio En 1 /(En 1 + En 2 ) and the CMT prediction from Eq.(36). Numerical values derive from the projection of an input field (E x (x), H y (x)), where E x (x) = e ik x x , over the exact (numerically calculated) Bloch Modes; predicted values are obtained as explained in the previous section. Also in this case relative errors are very small and we can see that the two modes carry nearly the same amount of power independently on tilt k.

In Fig. 12 and Fig. 13 we show the outcome of some finiteelement frequency-domain simulations of Maxwell equation in uniform and non uniform arrays. In order to test the previous approximations we inject into the arrays a narrow beam, corresponding to a broad spatial spectrum. If the array bands are flat, any plane wave component at the array input should be refracted and diffracted very little, entailing propagation without spatial broadening of the input beam. In Fig. 12 the input wavelength is set to 750 nm, the propagation length is 5 µm along z direction, and three different cases are compared. On the left propagation in a uniform array with ε 1 = ε 2 = 1.0 is simulated, and we notice that diffraction effects are much more accentuated than in the central and right figure, where non uniform case (ε 1 = 1.0, ε 2 = 1.5) is considered. In these two last cases, array bands are very flat and for a propagation distance of 5 µm they allow the field to propagate without distortion. In the central figure the core with dielectric ε 1 is excited whereas in the right one the core with dielectric ε 2 is excited. In the first case power is quite totally coupled on band 2, while in second case is quite totally coupled on band 1, so we could notice difference in propagation because the two bands have different amplitudes. Using a wavelength of 750 nm this difference is negligible for a propagation of 5 µm, but becomes noticeable when using a wavelength of 450 nm, as we can see in Fig. 13.

As far as losses are concerned, we verified that the propagation in the diffraction-engineered devices is not influenced at all by including a lossy model for the metal [ experimental values [START_REF] Johnson | Optical constants of noble metals[END_REF] gives ε(600nm) = -16 -0.44i], being the only effect a reduction in the transmitted power. Moreover the decay length of the fundamental mode of the waveguides with ε 2 = 1.5 dielectric core is, for example, L D (600nm) = [2Im(k z )] -1 = 11.95µm, much longer than the device length, indicating that all the relevant dynamics can take place without being suppressed by absorption.

We conclude this section reminding that, as said in section II, diffraction curves of the array are not sensitive on non uniformity consisting in the different width of adjacent cores: in Fig. 14 an example is shown in which cores have different widths 2b=120 nm (fixed) and 2c that sweeps from 120 nm to 180 nm. On the contrary, dielectric is constant in all cores and is set to 1.0. As we can see, diffraction curves change not too much even when 2b and 2c are very different.

VI. CONCLUSIONS

In this paper we studied non uniform plasmonic arrays. We have found that diffraction curves of the array modes are very sensitive to non uniformity consisting in the alternation of two different dielectrics ε 1 and ε 2 in the cores of the array, while not to non uniformity consisting in the alternation of different widths of cores.

In order to overcome problems arising from the CMT model whose basic cell is the single waveguide, we proposed a CMT model where the basic cell is the unperturbed coupler. We noticed that even for little non uniformity the shape of diffraction curves is sinusoidal and depends only on the couplings between equal modes in adjacent couplers; these couplings can assume different values, entailing different amplitudes for the two diffraction curves. We have seen that the trend of coupling coefficients as function of non uniformity (that is ε 2 -ε 1 ) depends quite exclusively on how modes of unperturbed basic coupler concentrate energy at its interfaces. Moreover, by using CMT, we have obtained a good approximation of Bloch modes of the array and of the amount of energies carried by the array modes as function of the input tilt.

In the last section, basing on all results obtained, we proposed an example of a flat-diffraction array in the visible band, showing the correctness of all predictions done in the previous sections.

VII. APPENDIX

In a single plasmonic guide formed by a dielectric layer 2b wide, with relative dielectric constant ε and surrounded by metal with relative dielectric constant ε m , the fundamental mode is TM-even one. Imposing the continuity of the tangential electric and magnetic fields, whose propagation constant is β, we find that:

Fig. 2 .

 2 Fig. 2. In the upper and in the center figure two single plasmonic guides are shown. The cores (width 2b) have dielectric constant ε 1 and ε 2 respectively, and are surrounded by metal with dielectric constant ε m . In the lower figure there is the coupler used as basic cell in the CMT developed in this paper. Two cores (width 2b) with dielectric constant ε 1 and ε 2 , and divided by a metal layer (width 2a, dielectric constant εm), are surrounded by metal.

Fig. 4 .

 4 Fig. 4. Solid line is the field Ex of mode a in a coupler with ε 1 =1.0, ε 2 =1.5, width of central metallic layer (ε m ) 40 nm and widths of dielectric layers 120 nm. Wavelength is 750 nm. Black dots denoted with A, B, C, D represent values of field at the interfaces of a coupler used in Eq.(17) and Eq.(18). They are respectively h e x(L,a) , h i x(L,a) , h i x(R,a) , h e x(R,a) . Double dashed arrows indicates h c and h i respectively, that are values of mode at the center and at the interface of core in guide g 2 of the coupler.

Fig. 5 .

 5 Fig.5. Approximation of coupler modes (transverse magnetic field). Exact modes (solid lines) are approximated with a constant over the two basic semiperiods of the array (dot lines). Constants h a1 , h a2 , h b1 , h b1 represent the mean value of modes over the correspondent semiperiod.a)mode 1; b)mode 2

Fig. 6 .

 6 Fig. 6. Amplitude of the diffraction curves as function of ε 2 , when ε 1 = 1.0. Solid lines for band 1, dotted lines for band 2. a) 750 nm; b) 450 nm

Fig. 7 .Fig. 8 .

 78 Fig. 7. Amplitude of the diffraction curves as function of wavelength. Solid line for ε 1 = 1.5, ε 2 = 1.5; dotted line for ε 1 = 1.0, ε 2 = 1.0; dashes thin line for ε 1 = 1.0, ε 2 = 1.5. a) band 1; b) band 2

Fig. 9 .Fig. 10 .Fig. 11 .

 91011 Fig. 9. Coupling coefficients Ca and C b as function of ε 2 , when ε 1 = 1.0. Solid lines represent values calculated with classic CMT model (Eq.(8)); dotted lines represent values C a (ε 2 ) + C a (ε 2 ) and C b (ε 2 ) + C b (ε 2 ), where C a (ε 2 ), C a (ε 2 ), C b (ε 2 ) and C b (ε 2 ) are approximated by Eq.(17) and Eq.(18) after calculation of C a (ε 1 ), C a (ε 1 ), C b (ε 1 ) and C b (ε 1 ). Wavelength is 750 nm. a) Ca; b) C b

Fig. 12 .

 12 Fig. 12. Propagation in the array with 2a=40 nm, 2b=120 nm. Single guide excitation. Wavelength λ is 750 nm. a) ε 1 = 1.0, ε 2 = 1.0; b) ε 1 = 1.0, ε 2 = 1.5, core ε 2 excited; c) ε 1 = 1.0, ε 2 = 1.5, core ε 1 excited.

Fig. 13 . 1 ]Fig. 14 .

 13114 Fig.13. Same as in Fig.12, but wavelength λ is 450 nm.

  T d b T m b = ε |ε m | tanh(T d b) (38)where T d and T m are the decaying constant of the field in dielectric and metal respectively, withT d = β 2 -(ω/c 0 ) 2 ε and T m = β 2 -(ω/c 0 ) 2 ε m . Substituting T d b = x and k 2 = (w/c 0 ) 2 (|ε m | + ε)b we obtain: x √ x 2 + k 2 = ε |ε m | tanh(x)(39)When ε/|ε m | < 0.5 we can approximate left and right hand side of Eq.(39) with (1/k)x and (ε/|ε m |)(1/x) respectively, and k with w p b/c 0 . In this way we found that x = εk/|ε m |, then we obtain approximated solutions for the decaying and propagation constants:T d ≈ εw p |ε m |bc 0