{'fraction_non_alphanumeric': 0.034358096920189866, 'fraction_numerical': 0.010348824373551165, 'mean_word_length': 4.354957883848086, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Supermodes of Hexagonal Lattice Waveguide Arrays

We present a semi-analytical formulation for calculating the supermodes and corresponding Bloch factors of light in hexagonal lattice photonic crystal waveguide arrays. We then use this formulation to easily calculate dispersion curves and predict propagation in systems too large to calculate using standard numerical methods.


I. INTRODUCTION
Coupled optical waveguides are an important class of devices for the manipulation and processing of light and include, for example, the widely used directional coupler [1]. In waveguide arrays, which can be seen as generalizations of the directional coupler, the propagation of light is more complicated. The spreading out of light into neighboring waveguides in an array can be seen as a discrete version of diffraction. By varying the angle of the incoming light, Eisenberg et. al. showed that the magnitude and sign of this discrete diffraction could be manipulated in ways not possible in regular diffraction [2]. The modeling of directional couplers and complex waveguide arrays can be achieved with fully numerical methods, but the ability to predict unique phenomena such as discrete diffraction, is possible only by harnessing the efficiency and physical insight gained by modal methods. Although an exact modal method exists, Yeh showed that much analytical insight can be achieved in the tight binding limit [3]. In this limit the overarching coupled waveguide modes or supermodes of the system can be represented as linear combinations of the single waveguide modes with coupling between adjacent guides determined by the overlap between the modal fields and the neighboring waveguides. This method shows that for uniform waveguides the even superposition is always the fundamental coupled waveguide mode. This is consistent with the oscillation theorem, according to which, in a 1D geometry, the mode with the fewest number of nodes has the lowest energy [4]. Yeh's method allows us to define the symmetry of the individual supermodes and accurately predicts the evolution of light in the individual waveguides within the array.
The coupling of photonic crystal waveguides (PCWs) has been a prominent area of research because of their short coupling lengths and unique dispersion properties, making them suitable for use in ultra compact devices. In square lattice PCW arrays, Locatelli et. al. showed that, depending on the number of rows between the waveguides, the coupling coefficient can change sign, inverting the sign of the discrete diffraction coefficient throughout the array [5]. The ability to invert the sign of the diffrac-tion has applications in fields of negative index media and imaging. Locatelli's work arises from the work by de Sterke et. al. into the coupled modes of square lattice PCWs [6]. They found that, as the number of rows between the waveguides varies, the symmetry of the fundamental mode alternates between even and odd because the underlying Bloch mode undergoes a π-phase change every transverse period. Therefore the order of the modes can invert, while still conforming to the oscillation theorem. Using a novel method combining the scattering matrix and tight binding methods, Botten et. al. generalized this formulation to arrays with arbitrary numbers of waveguides [7]. They found that, as the number of rows between the waveguides is increased, the order of the supermodes completely inverts, which implies that the diffraction coefficient has reversed. In that work a single propagating plane wave approximation is made, mapping the problem onto that of the homogeneous waveguide array. Through node counting it is shown that the order of supermodes in the square lattice PCW array upholds the oscillation theorem.
Despite the significant advantages of hexagonal lattice PCW systems, leading to their prominence in experimental research, modal methods for hexagonal lattice PCW systems are lagging behind their square lattice counterparts. This is because of the analytic complexities which arise when modeling the hexagonal lattice symmetry, for example, the non orthogonality of the lattice vectors (see Fig. 1(a)).
In hexagonal lattice coupled PCWs, the geometry of the system can either be in-line or staggered depending on the number of rows between the waveguides (see Fig. 1(b)). Ha et al. have shown that this staggered geometry is particularly efficient in the coupling of slow light between the waveguides [8]. In the staggered geometry, the coupled PCW modes are degenerate at the Brillouin zone (BZ) edge, where the two approaching modes have equal and opposite group velocities. Sukhorukov et. al. showed that this is associated with the presence of vortex modes, satisfying the Bragg condition [9]. In previous work, we show that additional degeneracies can occur within the BZ, where the coupled PCW modes braid around each other. By analyzing the modes of the cladding, we find a region in the BZ where the two least evanescent Bloch modes decay at the same rate, the difference in the Bloch mode periodicity leads to the beating in field, causing braiding [10].
Generalizing this theory to accommodate an array of coupled PCWs is not straightforward. There are a number of extra obstacles that arise in modeling this system. First, the advantage of using modal method analysis for complex problems, such as an array, comes from exploiting the symmetries of the system, but in the array these symmetries are much harder to find. To achieve the required symmetry in a system of more than two waveguides, we need to represent the field in each waveguide by one quantity which must symmetric about the waveguide center. But in this formulation, we wish to model systems where waveguides are created by modifying the radius or refractive index of a row of inclusions, rather than being completely removed. This means that we cannot treat the waveguide as a uniform medium-instead we must work completely in Bloch mode bases, in which we cannot easily relate the field at the center of the waveguide to its edges. Secondly, when comparing vertically aligned waveguide centers in the staggered geometry, we find they are at inequivalent points of the unit cell and therefore are not related by Bloch's theorem. Thirdly, since it has been shown there are two equally important Bloch modes in the hexagonal lattice [10], the treatment must be generalized to the vector case, necessitating a treatment using matrices rather than scalars. Lastly, the basis vectors are not orthogonal which means that the Bloch factors in directions parallel and transverse to the waveguide do not decouple.
For this reason we report here a formulation which allows us to calculate the modes of a arbitrary number of coupled PCWs in a hexagonal lattice. We consider the tight-binding limit, and, in spite of the complex geometry, find that the resultant matrix has the same symmetry as that of Yeh. The key difference is that the coupling between adjacent waveguides is expressed in terms of reflection and transmission matrices [7], rather than over-  lap integrals. This allows for a simple, numerically efficient, tractable modeling tool which provides the physical understanding to exploit the properties of coupled hexagonal PCWs. We show, for example, how the braiding behavior, previously shown for two coupled PCWs, generalizes to the PCW array. The outline of this paper is as follows: Section II presents a detailed derivation of the formulation, showing the main result of this paper, before Section III shows results achieved using the new formulation, comparing them to established methods and extending it to consider problems too complex to be practically calculated using standard numerical methods. Section IV sums up the results of this paper.

II. MODAL FORMULATION
In this section we derive a formulation for calculating the supermodes of a hexagonal lattice array with an arbitrary number of PC waveguides. We begin by deriving a resonance condition that may be used to find the dispersion relation and modes of a single PCW. We then calculate the reflection and transmission properties of the barriers between waveguides in the PCW array. Finally we use these quantities, together with a tight binding approximation, to relate the fields in each waveguide and derive a resonance condition that may be used to find the dispersion relation and supermodes of the entire PCW array.
We work exclusively in the Bloch mode bases of the two PCs that comprise the system of coupled waveguides; the Bloch mode basis is the natural basis for a periodic medium. At a single frequency, the field in a PC may be written as a superposition of propagating and evanescent Bloch modes [11], therefore we represent field by vectors of Bloch mode amplitudes. In any given unit cell we may write a vector c − of downward propagating/decaying mode amplitudes, and a vector c + of upward propagating/decaying modes.
In each bulk PC, the Bloch amplitudes in different unit cells are related by Bloch's theorem. In our notation, the first lattice vector a 1 = (d x , 0) is parallel to the waveguide. The other lattice vector is a 2 = (−d x /2, −d y ) (see Fig. 1(a)). Bloch's theorem may be used to define the Bloch factor, µ Rs = e iks·a R , which is the ratio of the amplitudes of mode s at two points separated by the lattice vector a R . Here k s is the Bloch mode's wavevector; in our method, its x component k x is specified for a particular frequency, and its y component is found numerically using a Rayleigh multipole method [12].
Due to an impedance mismatch at the PC's edges, the Bloch modes may couple at an interface between PCs [13]. Therefore in our calculations we need to consider a vector of several Bloch modes in each PC. The number of modes required for an accurate calculation varies and is related to the number of propagating grating diffraction orders; for the examples studied here, two Bloch modes in each direction are sufficient, but we include five modes for additional accuracy. The unique mode-braiding behavior of hexagonal lattice PCWs arises from the existence of equally dominant evanescent Bloch modes [10], and therefore fundamentally requires multiple Bloch modes to be modeled.
We can consider a PC as a stack of gratings and apply the grating diffraction equation. At an interface between PCs where the periodicity in x is conserved, coupling can only occur between modes with the same k x This means that for an entire system of modes, µ 1s = e ikxdx , which represents translation by a 1 , is fixed. On the other hand, µ 2s is different for each mode. Given a vector c − of downward propagating/decaying mode amplitudes in a unit cell, the amplitudes in the unit cell positioned by a 2 with respect to the first is Λc − , where Λ = diag(µ 2s ). For a vector c + of upward propagating/decaying modes, the corresponding vector in the unit cell displaced by a 2 is Λ c + , where Λ = diag(µ 2s ). It may be shown using reciprocity that for hexagonal lattice PC Λ −1 = Λe ikxdx [11]. When we relate Bloch mode amplitude vectors c ± in different unit cells, we generally do so for cells displaced by integer multiples of a 2 .
Each PC has its own set of Bloch modes. Across an interface between PCs ( Fig. 2(a)), the Bloch modes of the two PCs couple. Analogously to thin film optics, the outgoing Bloch mode amplitudes c + 1 and c − 2 are related to the incoming Bloch mode amplitudes c − 1 and c + 2 by reflection and transmission matrices, which may be calculated most conveniently from numerically calculated PC impedances [13]: A. The single waveguide problem First, we derive the resonance condition for a single PC waveguide, which may be used to find the PCW's dispersion relation and modes. This provides an introduction to our nomenclature and methods, as well as deriving a result we use in Sec. II C to find PCW array supermodes in which the field in the individual PCWs is approximately even. The PC waveguide consists of a central waveguiding region (PC 2) sandwiched between two PC mirrors (PC 1). The resonance condition is derived as a matrix analog to that of an optical slab waveguide. We do not assume that the waveguiding region is a uniform dielectric: we allow the possibility that it consists of a row of holes-this allows an extra degree of freedom when tailoring the dispersion of the PCW.
We write down the relationship between the Bloch mode amplitudes on each side of the waveguide: where c − is the vector of downward mode amplitudes at the top of the waveguide, and c + is the vector of upward mode amplitudes at the bottom of the waveguide (Fig. 2b). From the up-down symmetry surrounding the interface, we know R 21 = R 21 . The asymmetric factor e ikxdx arises due to the nonorthogonality of the two lattice vectors. Its presence becomes particularly inconvenient in Sec. II C, which requires a high degree of symmetry. We remove this factor by definingc The term x ROW corresponds to the shift in x associated with a translation by (−y R /d y )a 2 , where y R is the y coordinate of the unit cell for which c ± is defined. Note that the c + and c − in Eqs.
(2) are defined on opposite sides of the waveguide region ( Fig. 2(b)), and so their x ROW differ by d x /2: this is what removes the symmetry-breaking quantity e ikxdx . The original quantities c ± are defined in unit cells separated by multiples of the diagonal lattice vector a 2 . For every second row of unit cells, multiplication by e −ikxx corresponds to a translation by an integer multiple of a 1 , and soc ± has the straightforward interpretation of being the amplitude vector of the Bloch modes in the unit cell with x = 0 in that row. For the other rows,c ± has no direct interpretation, butc ± e ikxdx/2 is the vector of Bloch mode amplitudes at x = d x /2. This new basis for Bloch mode amplitudes necessitates the definitioñ whereΛ replaces Λ in relating the downwards amplitudes of Bloch modes in different rows of a PC. The reflection and transmission matrices are unchanged. With our new quantities, we can rewrite Eqs. (2) as relating shifted Bloch mode amplitudes (Fig. 2(c)) By substituting Eq. (5b) into Eq. (5a), we find the resonance condition for a single PCW mode, So, the condition for a single waveguide mode is The eigenvectorc − associated with the zero eigenvalue gives the waveguide mode in terms of the Bloch modes of PC 2, the in-band PC. Using Eqs. (5) and (7), we derivẽ where Eq. (9a) holds for even PCW modes and Eq. (9b) holds for odd PCW modes. In an experiment it is easier to couple light into symmetric modes, so we will concentrate on the even PCW modes. The contribution of even modes to the field in a waveguide may be summarized as which may be shown to represent the amplitudes of the even superpositions of Bloch modes in a single waveguide (Eq. (12)). To find the waveguide mode's field H(x, y) inside a waveguide unit cell, we convert our calculated result,c + andc − , back into the basis of c + and c − . If we choose c − such that it is defined with respect to the origin, then by Eq. (3) there is no difference between the two bases, i.e., c − =c − and the amplitudes of the upward Bloch modes in this cell areΛ 2c + . Therefore the magnetic field H swg at position r in the unit cell is where Ψ + (r) and Ψ − (r) are respectively vectors of the fields of the forward and backward Bloch modes. A similar expression may be written for the electric field. It is straightforward to use the transmission matrix T 21 to calculate the Bloch mode amplitudes of PC 1 fromc + , and a similar procedure may be followed to calculate the field of the waveguide mode throughout PC 1.

B. The transmission and reflection matrices across the waveguide barrier
In this paper we consider PCW arrays consisting of layers of two different photonic crystals, PC 1 and PC 2, where PC 1 is the bulk PC, which is in bandgap, and PC 2 is the waveguide . The barrier consists of rows of PC 1, and the waveguide is one row thick. PC 1 and PC 2 may be dielectrics or PCs, and need not have the same background refractive index.
This section focuses on a part of the array consisting of just two neighboring waveguides and the barrier separating them (see Fig. 3a), to find the reflection and transmission matrices across the barrier. We neglect reflections off the outermost edge of each waveguide, and aim to relate the vectorsc − j ,c + j ,c − j+1 andc + j+1 to find the reflection and transmission matrices in each direction across the barrier, whereR,T andR ,T are reflection and transmission matrices from above and below the barrier respectively. If we setc + j+1 = 0 (see Fig. 3a), then the only Bloch mode incident on the barrier is from above, i.e.c − j , and our relations in Eqs. (13a) and (13b) becomec + j = Rc − j andc − j+1 = Tc − j . The relations betweenc + j andc − j+1 and the Bloch modes inside the barrier,c ± b are given bỹ From the up-down symmetry, we know T 12 = T 12 and R 12 = R 12 . We can rearrange these relations to find c − j+1 in terms ofΛ 2c − j : Similarly, by rearranging the relations to findc + j in terms ofΛ 2c − j , we derivẽ If we instead setc − j = 0, then Eqs. (13a) and (13b) becomec + j = T c + j+1 andc − j+1 = R c + j+1 . By following the same method as above, we deriveT =T andR = R.

C. PCW array dispersion relation and supermodes
In a waveguide array, interactions occur between every pair of waveguides. For large arrays of many waveguides, it is unfeasible to model all these interactions. Following Yeh [3], we ignore the interactions between non-adjacent waveguides, making a nearest neighbor approximation, which is valid in the tight binding regime. In doing so, we assume that the field decays substantially across the barriers between waveguides, i.e., we assume that |µ | 1 for every element µ ofΛ 1 . Formally, we neglect terms of order O(Λ 2 1 ). Doing so excludes the interactions between non-adjacent waveguides, as well as the Fabry-Perot-like terms in the interactions between adjacent waveguides.
Under this approximation, the transmission and reflection matrices in Eqs. (16) and (17) simplify tõ The reflection matrix for the waveguide barrier is approximated by the reflection off an infinitely thick barrier since any effect from the neighboring waveguide consists of terms of order O(Λ 2 1 ). We can now relate the Bloch mode amplitudes of each of the waveguides in an M -waveguide PCW array. The vectorsc + j andc − j that represent upward and downward Bloch mode amplitudes in waveguide j are as defined in the previous section. The topmost and bottommost waveguides, 1 and M , interact with only one neighbor each. The relationship betweenc − 1 andc + 1 is similar to that of a standalone PCW (Eq. (5)); the difference is that c + 1 includes a contributionTΛ 2c + 2 from waveguide 2: Similarly, in waveguide M ,c − M includes a contribution from waveguide M − 1: The other waveguides 1 < j < M , interact with both their neighbors (see Fig. 3b): In these equations, we represent the field in each waveguide by two different quantities: on the top edge as the amplitude vector of the forward modesc − , then on the bottom edge as the backward modesc + . In order to apply methods analogous to those used with coupled dielectric waveguides [3], we need to represent the field in each waveguide by one quantity. We definẽ e j only represents the contributions of even superpositions of Bloch modes in the waveguide, as noted in Sec. II A. From here on we represent the amplitudes in waveguide j byẽ j . We start by using Eq. (22) to write the 2M equations (19)-(21) as M equations relating thẽ e j :ẽ In order to replace the residualc ± j in the last terms of these equations, we recall our approximation discarding terms of order O(Λ 2 1 ). The field in each waveguide is an even single waveguide mode, plus interactions of order O(Λ 1 ) from adjacent waveguides: thus for each waveguide in the array, using Eq. (9a) we writẽ c + j =c − j +O(Λ 1 ). Note that thec ± j terms in Eq. (23) are all multiplied byT, which is of order O(Λ 1 ). Therefore for these terms we can ignore the O(Λ 1 ) contributions from other waveguides and writec ± j ẽ j . Applying this substitution to Eq. (23), we write the equations solely in terms ofẽ j , in the form of a tridiagonal block matrix equation, where A = R 21Λ2 and B = 1 2TΛ 2 . This is the resonance condition for a coupled waveguide array, in the same way that Eq. (6) is the resonance condition for a single waveguide. For each frequency, a range of k x can be scanned for the existence of supermodes by testing whether the condition may be satisfied by each individual k x .
The difference between Eq. (24), and the expression for conventional waveguide arrays [3], is that each element in Eq. (24) is a n × n matrix rather than a scalar. If n Bloch modes are considered in each waveguide, then the matrix in Eq. (24) is of dimension nM . For a large number of waveguides M , it can be computationally expensive to solve the eigensystem of this matrix. Therefore, in a block-matrix analog to the conventional treatment [3], we exploit the tridiagonal block symmetry of the matrix and block-diagonalize it analytically using a similarity transformation. The problem then reduces to solving the eigensystems of M matrices, each with dimension n. Analogous to Yeh's treatment [3], the elements of the symmetric and unitary similarity transform Φ that block-diagonalizes the tridiagonal matrix are given by where I has dimension n and To demonstrate how to use Eq. (24) to calculate a dispersion relation, we give the example of the two waveguide case, for which This can be block diagonalized using a similarity transformation to give where, from Eq. (26) Φ = 1 (29) Equation (29) decouples into two nontrivial solutions: (I − (A + B))x 1 = 0 withx 2 = 0 (even supermode), and (I − (A − B))x 2 = 0 withx 1 = 0 (odd supermode). The coupled waveguide array's dispersion relation may be constructed from the frequency-k x pairs at which at least one of these conditions is satisfied. So the condition for an even supermode is det(I − (A + B)) = 0, and for this modeẽ 1 =ẽ 2 = 1 2x 1 , wherex 1 is an eigenvector of I − (A + B) with an eigenvalue of zero. This is the Bloch mode analogue of the result achieved in our previous work [10].
Due to the block tridiagonal symmetry of the matrix in Eq. (24), similar conditions may be obtained even for very large systems of coupled waveguides. Generally, the condition for the sth PCW array supermode is with all otherx i = 0, where σ s = 2 cos(sπ/(M + 1)). [3]. The fields in each waveguide for supermode s can be found by where Φ s denotes block column s of Φ The procedure to find the dispersion relation of a coupled PCW array therefore involves calculating A and B for frequency-k x pair and checking whether det(I − (A − σ s B)) = 0 for any σ s . At a frequency and k x for which this determinant is zero, then the null vector of The resonance condition Eq. (30) can be written directly in terms of the original reflection, transmission and propagation matrices, using Eq. (4): (32) The first two terms in Eq. (32) give the resonance condition of the single waveguide in a hexagonal lattice. Therefore the third term, of order Λ 1 , is the effect of the neighboring waveguides in the tight-binding approximation. In the photonic band gap, all Bloch modes in the cladding are evanescent and have Bloch factors |µ 1 | < 1. As we increase the distance between waveguides, increases, and Λ 1 → 0. Therefore the resonance condition of the array approaches that of the single waveguide as we increase the distance between the waveguides, as expected.

D. PCW array modes and propagation
At the k x associated with supermode s, Eq. (30) may be used to calculatex s . The otherx i = 0, and these values are used implicitly in Eq. (31) to find the supermode's field amplitudesẽ j in each waveguide j.
In choosing to represent the field in each waveguide by the vectorẽ j , we assumed that the field is an even superposition of upward and downward Bloch modes in each waveguide. In Yeh's treatment of coupled waveguides [3], the field in the waveguide array is written as a superposition of single waveguide modes, with one amplitude per waveguide. We do the same by writingx s ẽ, whereẽ is the vector of Bloch mode amplitudes calculated for a single waveguide. This approximation is consistent with an eigenvalue perturbation theory approach: k x is found including terms of order O(Λ ), but the modal field is found ignoring these terms.
From Eq. (31), it follows that the field of supermode s in waveguide j isẽ We now generalize to consider a superposition of supermodes with amplitudes w s , and calculate the resulting field in each waveguide,ẽ j =ã jẽ , which we simply represent by the scalar amplitudeã j . Theñ i.e.ã =φw, (34b) whereã and w are vectors ofã j and w s respectively. Propagation in the PCW array is most simply represented in the basis of its supermodes, each of which simply acquires phase as it propagates. We use the propagation constant k x of each supermode s to write its Bloch factor µ 1s = exp(ik x d x ) Given an initial field w(0) = φ −1ã (0) at the start of the waveguide array, the amplitude in waveguide j after propagating p periods is i.e.ã(p) =φ L p w(0), where L = diag(µ 1s ). Writing the initial condition w(0) in terms of waveguide rather than supermode amplitudes, we see that the field at any integer number of periods p along the PCW array is where φ = φ −1 . Therefore, knowing only the dispersion relation, we can analytically calculate the amplitudesã(p) of each waveguide mode at an arbitrary number of periods p into the structure. The actual field in each waveguide may be calculated fromã(p), and the field H swg (r) of the single waveguide mode over the domain 0 < x < d. Recall that in Eq. (3), we transformed the Bloch mode amplitudes c into a mathematically convenient basis and wrote them asc. Unlike in Sec. II A, we cannot set the origin of the coordinate system such that a j =ã j for all waveguides j, therefore we must explicitly perform a transformation into a more natural basis. We map eachã j to the amplitude a j of the single mode waveguide in the unit cell p periods into the structure, by writing where κ j = exp(ik x (x R − pd)), with x R the coordinate of the unit cell in waveguide j (more specifically, the point in that unit cell marked with a cross in Fig. 1(a), p periods into the structure. The origin of coordinates may be chosen such that x R = pd for every unit cell p along waveguide j = 1: then κ j = 1 for all waveguides that are in line with waveguide 1, and κ j = exp(ik x d/2) for waveguides staggered by ∆x = d/2 with respect to waveguide 1. Here, k x is approximated to that for the single waveguide mode. Finally, we write the field in the waveguide array as a superposition of translated single waveguide modes H swg (r − R j ) and their amplitudesã(p): where R j = (x R , y R ) is the coordinates of the nearest unit cell to r inside waveguide j.
If we neglect the evanescent tails of single waveguide modes, then we can write the field in or near waveguide j of the array in terms of one waveguide mode only: If we confine our interest to the field inside the waveguide at a single point per unit cell, then r − R j is constant and only one value H swg (r − R j ) is used; the mode may be renormalized such that H swg (r − R j ) = 1. Then the field throughout the waveguide array is given by the elements a j of a(p), which can be calculated analytically.

III. RESULTS
We now apply the methods we have developed to calculate dispersion relations for arrays of coupled waveguides, as well as the propagation of light through them. Where the use of conventional numerical tools is practical, we compare the results of our method to these. For simplicity, in all our examples we consider arrays constructed from the same two PCs. Both PCs are regular hexagonal lattices of cylindrical inclusions with radius r = 0.3 d, in a dielectric background with n = 3. The difference between the PCs is that the barrier, PC 1, has air holes with n = 1, whereas the waveguide material, PC 2, has inclusions with n = 1.5. Our method also applies to more conventional W1 PC waveguides, and can easily be extended to dispersion engineered waveguides: we consider an example with infiltrated holes in PC 2 in order to demonstrate that our method works even when the waveguiding region is not a uniform dielectric.
We consider two arrays based on this waveguide. The first is a system of M = 4 waveguides, each separated by a barrier of = 4 rows of PC 1. The second array is a system of M = 31 waveguides, also separated by = 4 rows of PC 1. The second array involves a structure with a width of over 150 unit cells, which is not feasible to simulate using conventional numerical tools. The separation of = 4 means that the waveguides are staggered.
We consider 5 Bloch modes in each PC, which means that for each frequency-k x pair, testing for the existence of a supermode only involves calculating the determinants of M = 4 different 5 × 5 matrices (Eq. (30)). At the frequencies we consider, we could get results of acceptable accuracy using only 2 Bloch modes.
In Sec. III A, we use Eq. (30) to calculate the dispersion relations of these structures. Then, in Sec. III B we calculate the discrete diffraction pattern of the field as it propagates through the waveguide array.

A. PCW Array Dispersion relations
First, we calculate the dispersion curves of the supermodes of the 4 waveguide system described above, for light polarized with the H field out of the plane, at normalized frequencies d/λ ∈ (0.29, 0.30). The dispersion curves are shown in Fig. 4. At most frequencies there are four supermodes based on the even single waveguide mode, and their dispersion curves agree well with those calculated using the Fictitious Source Superposition (FSS) method, which is considered to be highly accurate for this kind of problem [14]. Having established the accuracy of our method with a relatively small array of 4 waveguides, we now show its power by finding the dispersion curves of the supermodes of a 31-waveguide array (Fig. 5). We see that for this increased number of waveguides and supermodes, the dispersion curves of the supermodes start to form a band, but the dispersion curve still exhibits the braiding behavior shown for the coupled PCW system [10].

B. Driven PCW Arrays
We now turn our attention to field propagation through the two PCW arrays. We seek to simulate a waveguide array where the field is incident from a single waveguide. In the M = 4 waveguide array, we set a(0) = (0, 0, 1, 0) T , corresponding to an initial condition where the entire field is in the 3rd waveguide. At a normalized frequency of d/λ = 0.293, there are four propagating supermodes, at k x d = −2.68, −2.70, −2.74, and −2.80. A single waveguide with these parameters has k x d = −2.72 We choose the negative values of k x as they are associated with positive group velocity (Fig. 4). We calculate their field at one point per unit cell over the first 100 periods of the waveguide array, i.e., we calculatẽ a(p) from p = 0 to p = 100 (Fig. 6). The approximation made when writingx s ẽ has normalized residual errors of (2.2, 0.3, 0.3, 2.1) × 10 −3 for the four supermodes.
These results are indistinguishable from the results of an FSS simulation, although neither method simulates the coupling at the interface between the input waveguide and the waveguide array. In order to include coupling effects, we use a well established Bloch modal supercell scattering method [11] to model propagation through a single waveguide into the waveguide array. We sample the field at the same points in each unit cell for each waveguide in the array. There is good agreement between the supercell method and our results (Fig. 6). Now we consider the system of M = 31 coupled waveguides. Again, we start with all the field in the central waveguide, so all elementsã j ofã(0) are 0, exceptã 16 = 1. We calculate the propagation of the field through 1100 periods of the waveguide array (Fig. 7). For the first 450 periods, a typical discrete diffraction pattern is observed. After around 500 periods, the diffracted light reaches the outer waveguides, and edge effects start to play a significant role in the field profile. Once the

IV. DISCUSSION AND CONCLUSION
The semi-analytical formulation developed here provides an accurate modal method to solve for the modes of coupled PCW waveguides. An advantage of this method is that the computational time increases linearly as the number of waveguides increases, in the sense that for M waveguides, M equations of the form of (32) need to be solved to obtain the PCW modes at each frequency. The advantage here is that each mode is identified by its symmetry and solved for independently. Therefore, the memory requirements remain constant as M increases. This contrasts with supercell methods, where as the number of waveguides increases, the computation domain grows, increasing both the required memory and computation time. Furthermore, in supercell methods the modes of different symmetry are not readily identifiable, making crossings in the dispersion curves, such as those shown in Figures 4 and 5 difficult to recognize.
Our method is, in essence, a tight-binding formulation and therefore it warrants comparison to a conventional tight-binding formulation based on overlap integrals [3]. In such conventional formalisms one solves for the isolated single PCW mode and assumes that the modes of the coupled PCWs can be written in terms of a superposition of isolated single PCW modes. The interaction between adjacent PCWs is found by computing overlap integrals. The formulation presented here is identical in terms of the assumptions it is based on, i.e. adjacent waveguides couple weakly, however in place of overlap integrals interaction between waveguides is modeled by a transmission matrix. There is one key advantage in using our modal method based calculation: we find the PCW modes by finding the propagation constants k x satisfying (32) at a given frequency. In contrast, an overlap integral based formulation for PCWs gives the shift in frequency of supermodes from the single PCW mode. The former is much more desirable for modeling propagation problems, as once an operation frequency is chosen the relevant modes can be obtained in a single set of calculations. At the level of Bloch modes and supermodes, the similarities between the acoustic and electromagnetic treatments mean that this modal method may also be applied to arrays of elastic waveguides [15,16].
In the work of de Sterke et al [6], it is shown that although the fundamental mode can be either an even or odd superposition of PCW modes, when taking a cross section along the PCWs, the fundamental mode always possesses one less node than the second mode. Here we find that, when taking a cross section along the coupled PCWs, node counting cannot be used to predict the order of modes in the hexagonal lattice. This is because the oscillation theorem only holds for one dimensional Sturm-Liouville problems. In the case of square lattice PCWs, the modes have a frequency below the first Wood anomaly [6]. This allows for the use of a scalar approach where coupled PCWs are mapped on to coupled uniform waveguides. The rich physics found in hexagonal lattice PCW arrays, is due to the fact that there are two dominant Bloch modes and a vectorial approach is needed. We have observed that near the center of the Brillouin zone, where one operates below the first Wood anomaly, the oscillation theorem holds in the hexagonal lattice. However, as the Bloch wavevector approaches the Brillouin zone edge, crossing the first Wood anomaly and especially upon entering the braided region, the modes become inherently vectorial in nature, ie. their Bloch mode expansion requires at least two Bloch modes, and thus cannot be mapped on to a scalar amplitude like in the square lattice. Therefore there is no general node counting algorithm governing the order of modes in hexagonal PCW arrays.