{'fraction_non_alphanumeric': 0.043391458175385395, 'fraction_numerical': 0.014278493808440738, 'mean_word_length': 4.688277150866097, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Paraxial Hamiltonian for photons in two-dimensional photonic crystal microstructures

New solid-state physics based approach is developed for analysis of the paraxial light propagation in two-dimensional (2D) photonic lattices of coupled dielectric waveguides or microcavities. In particular, using Maxwell's equations, a non-Hermitian Hamiltonian eigenproblem with respect to the spinor wave function of a photon is obtained for energy-dissipating photonic microstructures. The Hamiltonian is suitable for almost the entire subclass of 2D structures encompassing arrays of semiconductor microcavities and microstructured photonic crystal fibers, characterized by light propagating mostly normal to the periodic lattice plane. Methods of numerical solution are discussed and the formalism is applied to a square array of coupled semiconductor microcavities, revealing reach possibilities for tailoring photonic band structure both in the photon energy and photon lifetime energy broadening domains. In particular, a feasibility to open a double photonic crystal band gap simultaneously in the energy and lifetime energy broadening domains is demonstrated.


INTRODUCTION
Photonic crystal structures offer unique possibilities for controlling light-matter interactions 1,2,3,4 and tailoring light propagation 5 by introducing lattice defects 6,7 or photonic crystals heterostructure barriers. 8,9 Among proposed so far photonic structure configurations, twodimensional (2D) lattices are highly attractive in technological aspect and have a realistic potential of finding applications in novel optoelectronic devices and integrated photonic circuits.
A particular subclass of 2D structures consists of arrays of coupled optical waveguides or microcavities. In these structures, a photonic mode propagates mostly along the waveguiding direction or cavity axis [vertical direction in Fig.1 (a)], such that only a small lateral component k ⊥ of wave vector k undergoes Bragg reflections. Such paraxial photonic structures employ lattices of period significantly exceeding the optical wavelength. Photonic crystal fibers 10 are an example of implementing this concept. Matrices of phase-coupled vertical cavity surface emitting lasers (VCSELs) are another example of such 2D structures. 11,12 These structures allow photonic wave functions behaviour to be directly examined in periodic and quasi-periodic lattices 13,14 and offer a possibility to implement photonic crystal heterostructures made of several photonic crystal materials with dissimilar band gaps. 9,15 In addition, they might be loaded with optical gain or loss, which are shown 16,17,18 to impact band edge energies at high-symmetry points of reciprocal lattice.
Many of these concepts related to photonic crystal structures originate from analogy to semiconductors (see for example, Refs. 19 and 20 in the case of author). However, despite an apparent simplicity of peculiar 2D photonic crystal microstructures presented above, a simple approach utilizing a standard Hamiltonian framework of the solid-state physics has not yet been reported so far for photonic crystals.
In this paper, a single-photon non-Hermitian Hamiltonian is obtained for 2D photonic lattices of parallel dielectric waveguides or microcavities. It is derived by introducing a paraxial gauge transformation, which converts Maxwell's equations into a non-Hermitian Hamiltonian eigenproblem with respect to a biorthonormal set of spinor wave functions for photons. (See Eq. (60), the main result of this paper.) Properties of this Hamiltonian in case of lattices with inversion symmetry (e.g., widely used square and triangular lattices) are discussed and a method of numerical solution based on biorthonormal plane wave expansion is outlined. The numerical solutions are obtained for the structures with square symmetry of the lattice, predicting a peculiar interplay of the photon energy and photon lifetime energy broadening bands. These results envisage new possibilities for application of dissipative photonic crystals heterostructures benefiting from the features of 2N -dimensional confinement of photonic envelop wave functions in Ndimensional photonic lattices.
For the first time, a truncated form of this Hamiltonian (adapted for the case of VCSEL arrays) has been briefly introduced in Ref. 21, without providing a discussion about its validity or proving orthogonality of its solutions. However, the symmetry and polarization structure of the main lasing modes predicted by the model has been confirmed in Ref. 21 by experimental measurements in square arrays of VCSELs. (See also Refs. 22,23.) A k·p approximation of this Hamiltonian restricted to the eight bands 20 has been solved analytically, yielding the dispersion characteristics of the low-order bands in square lattices. 24 A parabolic approximation obtained for the optical loss dispersion allows a single-band effective mass Hamiltonian to be introduced for modeling properties of the main lasing modes in VCSEL arrays. Its finite difference element representation with the numerical grid step being the array pitch has the form of Coupled Mode Theory equations, offering an efficient tool for calculating photonic envelope functions in complex quasiperiodic lattices. 9 This method, which is detailed in the later Ref. 25, has been used to interpret the experimental measurements of confined optical envelop functions in VCSEL-based photonic crystal heterostructures. 9,15 In Ref. 20, this Hamiltonian is used to analyze the Coriolis-Zeeman effect for photons in periodic lattices of microcavities, envisaging possibility of interactions between photons and gravitational field. The basic steps of analysis in Ref. 20 might be used as guidelines for following the treatment of photonic lattices in the present paper. However, note that in Ref. 20, the dissipative effects are not taken into account.
The paper is organized as follows. In Sec. II, using an equivalent cavity-unfolded representation of microcavities, the Hamiltonian for lattices of microcavities defined by mirror reflectivity patterning is obtained. In Sec. III, the Hamiltonian is extended to the case of microcavities or parallel dielectric waveguides defined by periodic variations of refractive index. Sec. IV brings together the non-Hermitian Hamiltonian terms obtained in previous sections, discuss the validity of such Hamiltonian approximation, outlines the method of solution based on biorthonormal plane wave expansion and reports numerical results of band structure computations in squarelattice photonic crystal microstructures.

II. REFLECTIVITY-PATTERNED MICROCAVITIES
The Hamiltonian for lattices of microcavities with periodic variations of mirror reflectivity is detailed here on example of mirror-patterned VCSEL arrays (Fig. 1). In Sec. II A, a typical structure of VCSEL array with reflectivity-patterned distributed Bragg reflector (DBR) is introduced. In Sec. II B, the constitutive equations in cavity-unfolded representation are obtained. Sec. II C reports the model Hamiltonian.

A. VCSEL array photonic crystal
Arrays of vertical cavity surface emitting lasers are an example of photonic crystal made of evanescently coupled microcavities. These particular photonic structures utilize a usual VCSEL wafer incorporating a one-wavelength cavity with an optical gain medium (e.g. quantum wells) sandwiched between two DBRs [ Fig.1 (a)]. The λ-cavity might be regarded as a defect in the thick λ/2-periodic DBR stack. In a solitary VCSEL, it defines the lasing mode wavelength, which fits the cavity roundtrip self- repetition conditions. In VCSEL arrays, the λ-cavity fixes the main (longitudinal) wave vector component [along the z-axis in Fig.1(a)].
The photonic lattice of VCSEL array is defined by reflectivity patterning, R(x, y), of the top DBR. 11 The crystal unit cell consists of a highly reflecting VCSEL pixel surrounded by a grid of lower reflectivity. In the bottom-emitting structures, this is accomplished by depositing a periodic pattern of metallic overlays on the top DBR, e.g., using gold pixels and chromium grid for pattern definition. 11 In the top-emitting structures, the pixel positions and the grid are defined by air openings in the metallic film (e.g., gold film) deposited on the top DBR. 26 The pattern is characterized in terms of a reflectivity contrast and lattice cell fill factor F F (the area ratio of the pixel and of the lattice cell). The resulting composite DBR is of high reflectivity at the pixel positions (R 2 ∼99.9%÷99.99%) and of low reflectivity contrast between the pixels and the grid (2δR∼1%). Since the roundtrip optical gain is low and uniformly distributed across the VCSEL structure, such shallow reflectivity pattern suffices to define position of the lasing microcavities.
For a typical InGaAs/GaAs/AlGaAs VCSEL array emitting at 960 nm wavelength, the lattice pitch Λ is about 5 µm, such that only a small transversal component k ⊥ ∼π/Λ of the propagation vector k (k ⊥ /k ∼ 0.03) undergoes Bragg reflections in the plane of periodic lattice. The angular spectrum of the most important loworder cavity modes is thus of narrow width ξ ∼ k ⊥ /k and centered about the main propagation direction (the zaxis direction). In what follows, the cavity modes are an-alyzed in the paraxial approximation, within an accuracy of second-order terms ξ 2 ∼10 −3 . The impact of reflectivity pattern R(x, y) is taken into account as a second-order perturbation. (For a typical VCSEL array, the amplitude reflectivity variations assume δR < ∼ 0.5 · 10 −3 .) The optical gain is uniformly distributed across the cavity 27 and has no influence on the cavity mode structure. Uniformly distributed optical losses, such as the material losses (∼ 0.1% in VCSEL arrays) and the output coupling losses at the bottom DBR (∼ 0.16%), do not influence a curvature of the dispersion curves for the cavity modes as well. Their additive contribution to the optical mode frequencies and losses is ignored in the subsequent analysis and can be easily taken into account by shifting the frequency and loss spectra of optical modes.
An optical mode of VCSEL array can be regarded as a superposition of standing waves coupled via Bragg scattering effects at the patterned mirror. The partial amplitudes of standing-wave harmonics are defined thus by the boundary conditions at the cavity mirrors and not as a result of propagation within the structure, like in usual photonic crystal materials exhibiting periodic variations of the dielectric constant. To apply the usual method of orthogonal plane wave (OPW) expansion, 28 an equivalent cavity-unfolded representation is needed.
Multiple reflections at the cavity mirrors effectively translate the entire cavity into a structure that is periodic along the cavity axis (z-axis). The unfolded PhC is three-dimensional (3D) and can be analyzed in terms of propagating plane waves. Optical modes of a VCSEL array are thus represented by electromagnetic Bloch waves propagating in the equivalent 3D photonic crystal.
A typical VCSEL wafer incorporates a few tens of various dielectric material layers. To take into account a detailed composition of the wafer structure, the cavity has to be unfolded using the outermost layer interfaces. Thus reflections at the metallic overlay pattern in the top DBR and at the wafer substrate in the bottom Bragg reflector have to be used for unfolding the VCSEL cavity [ Fig.1(a)].
At each reflection from cavity output coupling interfaces, the field amplitude reduces because of the waves leaving the cavity. In unfolded photonic crystal, the impact of subsequent reflections is reproduced by a periodic stack of energy dissipating layers (g-layers), each representing the effect of a single mirror reflection [ Fig.1(b)]. The zero propagation length between the incident and reflected waves at the output coupling interfaces of the cavity assumes that the equivalent g-layers are infinitesimally thin. Introduction of such thin absorbing layers allows both the energy dissipation (within the layers) and the continuity of electromagnetic field components (at each layer interface) to be taken by the model into account.
In between the g-layers, the unfolded photonic crystal has the same sequence of dielectric layers as encountered in the VCSEL microcavity. However, such a detailed description of wafer structure composition is needed for analyzing technical issues, like, e.g., optimizing a spectral overlap between the gain peak and the cavity modes in function of the active region temperature. 27 This description leads to intricate numerical model demanding lengthy numerical simulations. This paper is focused on general features of the light propagation behaviour in photonic lattices. Therefore, a simplified model system is used 21,22 [ Fig.1(c)]. It consists of the λ-cavity of VCSEL array in which the DBRs are replaced with ideal mirrors of effective DBR reflectivities. In this way, a Fabry-Pérot cavity is formed in which the reflectivity of the upper mirror is modulated in two directions parallel to its plane. Such Fabry-Pérot resonator replicates the boundary conditions for the field in VC-SEL cavity, reproducing thus the effect of fast oscillations within the periodic DBR stacks. Introduction of this model system is consistent with the Kapitza method 29 to analyze oscillations in a dynamic system. In this particular case, it consists in separating fast longitudinal oscillations of the electromagnetic field (along the cavity axis) and its slow lateral oscillations (in the x-y plane).
The model system treated here consists thus of a Fabry-Pérot cavity of the length L=λ/n with patterned reflectivity R(x, y) of one mirror [the upper mirror in Fig.1.(c)]. The second cavity mirror is of uniform reflectivity and does not impact the cavity mode structure. In what follows, it is assumed to be perfectly reflecting. The dielectric material in the cavity is of uniform refractive index n and impedance Z. The equivalent cavity-unfolded representation is a simple 2L-periodic (along the z axis) structure shown in Fig.1(d).
Note that the periods (along the z axis) of the cavityunfolded VCSEL array structure [ Fig.1(b)] and its simplified model [ Fig.1(d)] are different. As shown below, the absorbing g-layers have the same impedance as the neighboring dielectric materials such that Bragg reflections in the z-axis direction cannot occur at the g-layers. Therefore, the period of the cavity-unfolded structure has no impact on photonic bands. Furthermore, the longitudinal wave vector component of a mode is set by the self-repetition condition at the cavity roundtrip. However, in the cavity-unfolded representation, it is not defined until such additional condition is superimposed, yielding k z = 2π/L.

B. Material equations for a cavity-unfolded structure
There are many different conventions for the coordinate systems used to consider reflections, and the phases of reflection coefficients are dependent on the coordinate system. Here, the following convention is used for the complex vector amplitudes of the incident and reflected waves: whereσ is the reflection operator, R s,p are the amplitude reflection coefficients, indexes s and p stand for polarization perpendicular to the plane of incidence and in the plane of incidence, respectively. The convention (1) yields equal (in modulus and phase) reflection coefficients R s and R p at normal incidence, R s = R p (= R). Furthermore, in the range of propagation angles ξ=k ⊥ /k z considered here, the angular dispersion of the DBR reflectivity is negligible, such that the approximation R s =R p is valid at oblique incidence as well. Indeed, for ξ=0.03, the reflectivity of a typical DBR decreases on δR∼10 −5 . (The reflectivity amplitude at normal incidence is R∼1.) The notion of absorbing g-layers in the cavity-unfolded representation of a structure is linked to the duality of symmetry operations. 30 In the active notation, the reflection operatorσ in (1) assumes transformation of the vectors E(x, y, z) and H(x, y, z). In the passive notation this symmetry operation is applied to the coordinate axes. The effect of Rσ on the field can be then regarded as a result of propagation through a layer with equivalent transmittance T =R. To match the phases of reflected (at the cavity mirror) and transmitted (by the equivalent glayer) waves, the layer thickness has to be infinitesimally small. To impact the amplitude and phase of the transmitted waves in this limit, the material parameters of the layer should vary as ∝ h −1 , with h being the layer thickness.
In the cavity-unfolded representation, electromagnetic waves propagate in the periodic stack of absorbing glayers [ Fig.1(d)]. Material equations for this structure are obtained here by considering the amplitude of a plane electromagnetic wave in Fabry-Pérot (FP) cavity with mirrors of uniform reflectivity. The cavity mirror located at z=0 is of reflectivity |R|<1. As in Sec. II A, the second cavity mirror (at z=−L) is a perfect reflector. In the cavity-unfolded representation, the amplitude of the wave traveling towards the positive z-axis direction reduces by a factor of R at each g-layer: where ω is the frequency of electromagnetic wave, E 0 and H 0 are the amplitudes of electric and magnetic fields at the coordinate origin z=0, θ(z) is the unit step function θ(z)= {0 (z<0), 1/2 (z=0), 1 (z>0)}, and j is an integer number (j=0, ±1, . . .) enumerating g-layers (at z=0, ±2L, . . .) associated with the subsequent reflections at the cavity mirror. The electromagnetic wave (2) is written in the form of the Bloch wave composed of plane wave envelope function and periodic crystal part (the second exponential term is of the period 2L). The effective propagation vector k of the wave (2) is complex: where ωn/c and τ are the wave number and unit vector along the wave propagation direction as seen from the cavity, andẑ is the z-axis unit vector. Eqs.(2)- (3) show the evolution of wave amplitude and phase with successive roundtrips in the cavity. The imaginary part of k accounts for the decay of electromagnetic field in the cavity due to output coupling loss. The real part of k allows the cavity mode frequencies to be evaluated from a usual condition Re (2k z L) =2πq with q being an integer. In the cavity [ Fig.1(c)], and hence in between the glayers of unfolded photonic crystal [ Fig.1(d)], the electric and magnetic field components (2) assume the relationships A correspondence between the electromagnetic energy flux in the cavity and in the equivalent unfolded structure can be established assuming that the impedance of g-layer material is the same as in the dielectric material of the cavity. Otherwise, an impedance mismatch and Bragg scattering effects at periodic stack of g-layers will result in a retro-reflected wave propagating in the backward z-axis direction, having no counterpart in the FP cavity. The relationships (4) thus hold through the entire equivalent photonic crystal, including the absorbing g-layers as well.
within the accuracy of the ξ 2 -order terms, one obtains the constitutive relationships for the periodic structure representing cavity-unfolded Fabry-Perot resonator where ε=n/Z and µ=nZ are the dielectric constant and magnetic permeability in the cavity, z −∞ δ(ξ)dξ=θ(z). Eqs. (6)-(7) assume the impedance matching through the entire structure, which thus exhibits no artificial photonic band gaps in the z-axis direction.
Eqs. (6)-(7) are similar to the constitutive equations in a nonpermanent gravitational field induced by rotating coordinate frame. 31,32,33 Indeed, in Eq. (1), the reflection operatorσ at the cavity mirror includes the coordinate rotation (σ =ÎĈ 2 (ẑ) whereÎ represents the coordinate inversion andĈ 2 (ẑ) represents rotation by π about the cavity axis). The non-Galilean space-time metric 34 (g = 0) within the layers represents the effect of mirror reflections and can thus be attributed to the coordinate rotations used to unfold the cavity.
An apparent discontinuity of the electromagnetic field at g-layers [see Eqs. (2) and (7)] is in fact a result of infinitesimally small g-layer thickness. Using Eq. (7) and approximation δ(z)= lim h→0 in a layer of small finite thickness h (e.g., located at the coordinate origin). It can be seen that the tangential fields E and H as well as the normal components of D and B are continuous at each interface of the g-layer, in agreement with the boundary conditions at moving (rotating) interfaces. 35,36,37 Within the layer, the tangential components of D and B vary as ∝ h −1 , yielding a singularity in the limit h→0.
The effective refractive index in the layer varies with the wave propagation direction 20,33,38 τ : For a plane wave propagating in the positive z-axis direction (ẑ·τ >0), the phase accrual in the layer reads ω c n eff τ dr, yielding a complex number where all other terms vanish in the limit h→0. At each g-layer, the wave amplitude thus reduces by a factor of e iφ = R, in agreement with Eq. (2). The wave propagation in the positive z-axis direction of the structure (6)-(7) corresponds to multiple roundtrips in the FP cavity. However, this periodic structure is nonreciprocal. Thus, the amplitude of a wave propagating in the backward direction (τẑ < 0) increases by a factor of R −1 (where |R|<1) at each g-layer. This picture corresponds to an external electromagnetic field exciting oscillations in the FP cavity. The opposite zaxis directions in the cavity-unfolded structure (6)- (7) are thus related by the time reversal operation.
In the case of mirrors with uniform reflectivity, the vector g in (7) shows the same transformation properties as vector g = 1 c Ω×r composed of space-time components of metric tensor in the case of rotations. Namely, g is anti-invariant (g → −g) under the coordinate inversion (P ) and time reversal (T ) operations. It can be seen that Eqs.(6)-(7) do not show any particular symmetry under the P T -transformation (coordinate inversion followed by time reversal operation). Therefore, the energy spectrum assumes complex eigenvalues, as opposed to the real spectrum in the case of pseudo-Hermitian systems that are invariant under the P T -transform. 39,40 However, the P T transformation is important for the analysis presented here. Thus, in Sec. IV, square-lattice structures are treated using non-Hermitian Hamiltonian and biorthonormal set of wave function partners. The invariance of square lattice under (2D) coordinate inversion allows the partners of biorthonormal set to be defined as Ψ(t, x) andΨ(t, x) = Ψ * (−t, −x). (In the case of dissipative effects, the results of substitution t → −t and complex conjugation are different, 41 such that a wave function Ψ(t, x) and its T -transform Ψ * (−t, x) do not coincide.)

C. Hamiltonian
The material equations (6)-(7) define a cavity-unfolded periodic structure, which is equivalent to a Fabry-Pérot resonator with mirrors of uniform reflectivity. The vector g in (7) does not vary with the position in the x-y plane (parallel to the cavity mirrors). Different from FP cavity, lattices of coupled microcavities, such as VCSEL arrays discussed here, utilize mirrors with reflectivity patterning R(x, y). The cavity-unfolded representation of such mirror-patterned structures is assumed to obey the same relationships (6), however, with periodically varying (in the x-y plane) space-time coupling where the z-period of unfolded structure is 2L, and R(x, y) is the (amplitude) reflectivity of the cavity mirror (e.g., effective reflectivity of the top DBR in VC-SEL arrays). Note that in the model defined by Eqs. (6) and (11), like in the case of other photonic crystal structure models, the field radiation effects at the boundary discontinuities 42 are neglected. It can be seen that an electromagnetic Bloch wave propagating in photonic crystal with effective noninertiality (11) is of the same form as indicated in Eq. (2) but the field amplitudes E 0 and H 0 vary periodically in the x-y plane. In (2), the frequency ω (real number) is the independent parameter of motion, while the Bloch vector k (complex) is defined from the dispersion equation (3) and takes the effective propagation loss in the structure into account [last term in the right-hand side of Eq. (3)].
In fact, the independent variable parameterizing the dispersion curve can be chosen either as the wave number or as the frequency of electromagnetic wave. In what follows, to obtain the standard Hamiltonian form of equations, the vector k is used as a real-valued independent parameter. Respectively, the quantityhω qk is a quantized observable assuming complex values (q is the band index). The real part of observablehω qk is the photon energy and the imaginary part yields the lifetime broadening due to the optical loss.
Maxwell' equations in a photonic crystal with periodically varying noninertiality (11) are solved here by separating the fast and slow field oscillations in, respectively, longitudinal (z-axis) and lateral (x-y plane) directions: where the index q enumerates the photonic bands, the scalar function f qk (z) and the set of two vector functions e qk (r ⊥ ) and h qk (r ⊥ ) are, respectively, the fast and slow oscillating Bloch wave components. In (12), the field amplitudes e qk and h qk are normalized in vacuum. As shown below, such separation of variables is valid in conditions of the paraxial approximation and low contrast of photonic crystal lattice. As in (2), the longitudinal part of the wave (12) assumes a general form where the first and the second terms are, respectively, the plane wave envelope function and periodic Bloch function part. The complex parameter φ qk defines the magnitude of periodic variations in the Bloch function amplitude and phase. This parameter is of the order of ln R ∼ ξ 2 . Therefore, within the accuracy of the ξ 2 -order terms, the wave (13) assumes the expansion: where The weak periodic modulation η qk (z) shows step-like variations at the positions of g-layers and thus represents the effect of abrupt phase-amplitude variations at each reflection of the mirror. Note that the z-period average values of the 2L-periodic function η qk (odd function) and its derivative ∂η qk /∂z (even function) are null ( 1 2L L −L η qk (z)dz= η qk 2L ≃0 and ∂η qk /∂z 2L ≃0). By virtue of these properties, the longitudinal and lateral Bloch function components of electromagnetic wave (12) can be analyzed separately.
Substituting (12)- (15) in Maxwell' equations for the curl of E and H, and taking the z-period average, one obtains the equations for the slow component propagating in the lateral (x-y plane) direction: Then, by taking the difference between the z-averaged equation in (16) and corresponding Maxwell' equation, one gets the equations for the fast oscillations in the zaxis direction Eqs. (16) can be converted into a Hamiltonian eigenproblem, provided that a photonic wave function ψ qk (r ⊥ ) is introduced via gauge transformation of the fields e qk (r ⊥ ) and h qk (r ⊥ ).
Wave states of a photon are usually expressed in terms of a three-component spinor wave function 43 (a photon spin is one). The spinor indexes represent polarization state, while the spatial distribution of spinor components accounts for the angular momentum. Thus, formally, three spinor components (with indexes s z =0, ±1) have to be taken into account in Eqs. (16)- (17).
However, a state with s z =0 cannot be realized in free space 43 and only two spinor components (s z = ± 1) are independent. Thus, by virtue of the transversality of electromagnetic wave in vacuum, only two field components might have independent spatial distributions. Obviously, the same holds for homogeneous dielectric media or periodic photonic crystal structures. In the last case, the wave transversality condition is replaced by the coupling of the field components via Bragg scattering effects.
Therefore, a photonic state is defined here in terms of two independent spinor components forming thus a twocomponent vector wave function ψ qk (x, y). For paraxial wave (12), it is convenient to choose the components ψ qk (x, y) of wave function ψ qk (x, y) in the lateral (x-y plane) direction. A perturbation analysis (not shown here) reveals that within the accuracy of ξ 2 -order terms, the wave function components can be introduced in Eqs.(16)-(17) via operator relationships where (ẑψ qk ) = 0,Ê is the tensor operator and twice repeated Greek indexes indicate summation over the x, y and z components. An explicit form of expressions (18)- (19) reads Eqs. (18)- (19) are used here as a definition of the gauge transformation introducing a wave function ψ qk . This gauge transformation is validated in Eqs. (25) and (37). In the gauge (18)- (19), the electromagnetic Bloch wave (12) reads qk (x,y), (21) where e αβγ is the completely antisymmetric unit tensor of Levi-Civita. The functions ψ qk andẑ×ψ qk defines the electric and magnetic fields of the main polarization component. The first-order terms (∝ 1 kz ∂ ∂x ∼ξ) and the second-order terms (∝ 1 k 2 z ∂ 2 ∂x 2 ∼ξ 2 ) contribute in longitudinal and cross polarization components of the wave. Eq. (21) allows a complex polarization structure of inhomogeneous wave to be taken by the model into account and is in agreement with the results obtained for Gaussian beams. 44 It must be noticed that ψ qk (andẑ×ψ qk ) differs from the lateral (x-y plane) component of the field e qk (h qk ). However, the squared modulus |ψ qk (r ⊥ )| 2 characterizes the energy flux in the z-axis direction. Thus, taking a z-period average [ Fig.1(d)] of the Poynting vector In the paraxial approximation considered here, the energy flux in the z-axis direction is thus defined by the main polarization component of the wave. The longitudinal and cross polarization components contribute to the energy flux in the lateral direction.
As an example, consider a wave function ψ k = ψ 0 e ikxx+ikyy with ψ 0 being a constant vector in the x-y plane. It defines a plane wave propagating in the direction of k = (k x , k y , k z ). The energy flux associated with the wave is of the density S z = c 8π |ψ 0 | 2 along the z axis. The approximation (21) takes into account the transversality condition for electric and magnetic components of a plane wave. Thus, for the case of p-polarized (TM) wave (k ⊥ ×ψ 0 =0), the electric field has nonzero components e 0z =−k ⊥ /k z ψ 0 and e 0⊥ =(1 − 1 4 k 2 ⊥ /k 2 z )ψ 0 , yielding the field amplitude e 0 =(1+ 1 4 k 2 ⊥ /k 2 z )ψ 0 . The magnetic field oscillates in the xy-plane with the amplitude h 0 =(1+ 1 4 k 2 ⊥ /k 2 z )ψ 0 . The energy flux in the direction of k is thus S = S 0z (1+ 1 2 k 2 ⊥ /k 2 z ), in agreement with the vector calculus utilizing directional angle ξ=k ⊥ /k z of the wave. By substituting e → −h and h → e, a similar agreement can be readily proved for the case of s-polarized (TE) wave (k ⊥ ψ 0 =0).
Within the accuracy of the ξ 2 -order terms, the inversion of Eq. (20) is straightforward, yielding the expressions for the wave function ψ qk in terms of the fields e qk andẑ × h qk : These relationships can be represented in the tensor operator form that reads where 45 e 3αβ e 3γβ =δ αγ − δ α3 δ γ3 , and the inverse operator Thus,Ê is not the unitary operator andÊ −1 =Ê + .
The paraxial gauge transformation (21) (18) in (16) and applying, respectively, the operators e 3αβÊ −1 βγ andÊ −1 βα , in the first and second equations [in (16)] one obtains (25) where m 0 =nhk z /c is the effective mass andp ⊥ =−ih∇ ⊥ is the momentum operator in the lateral (xy-plane) direction. The paraxial gauge transformation (21) representing an arbitrary photonic state by a function ψ qk (r ⊥ ) is thus validated by the fact that Maxwell' equations for E and H take the same form in this paraxial gauge. The first term in the HamiltonianĤ [right-hand side of Eq. (25)] is associated with the paraxial propagation along the cavity z axis and accounts for the dispersion of the longitudinal wave vector component in the dielectric material of the cavity (m 0 c 2 /n 2 =hk z c/n). The in-plane kinetic energy (second term inĤ) and effective potential (third term) take into account the dispersion and Bragg scattering effects in the lateral direction due to periodic reflectivity pattern R(x, y) of the cavity mirror. The parameter m 0 can thus be interpreted as the lateral effective mass of a photon in an empty lattice (in the case of R=1).
The non-stationary Schrödinger equation for ψ qk (t, r ⊥ ) follows from (25) by substitution i∂/∂t → ω. In the case of R=1, within the accuracy of the time variable, the Schrödinger equation (25) is analytically similar to the scalar paraxial wave equation. However, unlike a scalar field amplitude in paraxial wave equation, the spinor function ψ qk in Eq. (25) cannot be associated directly with any of the six components of electromagnetic field. The same remark applies to comparison between the Eq.(25) and the scalar 2D Helmholtz equation for microcavities. 46 The periodic potential U (x, y)=i ch 2nL ln R(x, y) implies that ψ qk is a Bloch wave 22 composed of plane wave envelope function and periodic crystal part where spinorial and angular parts of periodic function u qk (r ⊥ ) allow the impact of photonic lattice symmetry to be taken into account in the analysis of electromagnetic field behaviour under the lattice rotations. Since the polarization anisotropy of mirror reflectivity R(x, y) is much smaller than the order of effects accounted for in the paraxial approximation [see Eq.
(1)], the spin-orbit coupling term is neglected in Eq. (25). All eigen states of the Hamiltonian (25) are thus doubly degenerate by polarization. As shown in Ref. 20, this degeneracy can be removed by a symmetry breaking effects in nonperma-nent gravitational field, when the photonic crystal rotates along the cavity z axis. Finally, the HamiltonianĤ in Eq. (25) is non-Hermitian (Ĥ + =Ĥ). The general properties of non-Hermitian Hamiltonians were extensively studied in the past, yielding the conditions for a discrete real-valued spectrum of eigensolutions. 39,40 In Ref. 47, the interested reader can find a comparison between the eigenproblems of Hermitian and non-Hermitian Hamiltonians. In the case discussed here [Eq. (25)], the most important results of these studies are related to the spectrum ofĤ and orthogonality condition of its eigenfunctions.
The Hamiltonian (25) does not exhibit a pseudo-Hermiticity, as opposed to a class of P T -symmetric Hamiltonians, which are invariant under the time reversal followed by coordinate inversion. Therefore, the spectrum ofĤ assumes non-paired complex eigenvalueshω qk . Their complex conjugateshω * qk can only be obtained in the spectrum of the adjoint operatorĤ + . The real part of eigenvalues Re(hω qk ) is the photon energy and the imaginary part is the energy broadening due to a finite lifetime of photons in the cavity, Im(hω qk ) =h/2τ qk .
The orthogonality of eigensolutions can be established using the biorthonormal set of functions 48 consisting of the concomitant partners ψ qk andψ qk of associated with the eigenproblemsĤ andĤ respectively. Taking the difference between the matrix elements of Eq. (27) and the Hermitian adjoint of Eq.(28), one obtains This relationship evidences 49 the orthogonality of functions associated with different eigenvalues. In the case of Hamiltonian (25), the spinorial structure of wave functions ψ qk andψ qk as well as the plane wave envelopes and periodic Bloch function parts 50 have to be taken by the orthogonality relationship into account, yielding Here, the integration runs over the entire x-y plane and qk (x, y) is the scalar product evaluated at a point r ⊥ =(x, y). For photonic states at the same point of the 2D Brillouin zone (k ′ ⊥ =k ⊥ ), the integral over the entire crystal (30) can be reduced to a lattice-cell integral 50 where Ω ⊥ is the lattice cell area. For a state |ψ qk , the probability distribution function for the coordinates is thus given by the productψ * qk ·ψ qk , as opposed to the usual expression |ψ qk | 2 .
Indeed the probability to find a photon at a point In the Schrödinger picture, this expression allows the stationary probability distribution to be associated with an eigen state. By contrast, the distribution |ψ qk | 2 varies in time as ∝ exp (−t/τ qk ). This behaviour is expected by the photon lifetime considerations, since |ψ qk | 2 defines the the density of the energy flux along the cavity z-axis. Note that in the case of Hermitian Hamiltonian [Im(U )=0], both distributions are stationary and indistinguishable.
Substituting a Bloch wave (26) for the ψ qk and using the expressionψ qk = e ik * ⊥ r ⊥ũ qk (r ⊥ ) for its concomitant partner, we obtain the orthogonality relationship for the periodic Bloch functions 51 The orthogonality of partners ψ qk andψ q ′ k ′ leads to the orthogonality relationship for the lateral components of the electric and magnetic fields (18). Since the gauge transformation operatorÊ [Eq. (19)] is independent of the material equations, the electromagnetic field associated with the partnerψ qk is defined by the relationships The scalar product in the integrand of Eq. (30) can be then expressed in terms of the electric and magnetic fields Substituting Eq. (24) forÊ −1 and using the relationship where twice repeated index α indicates summation over the x, y and z vector field components. In the right-hand side of Eq. (35), the integrands in the second and third terms are of the form (∇ ⊥ ·F) with F being a vector field given by one of the expressions in curly brackets. As in Eqs.(30)-(31), the plane wave envelope function and periodic Bloch part 50 have to be accounted for in the electric and magnetic field components contributing to F. It follows that F can be represented as a product F=e i(k ⊥ −k ′ ⊥ )r ⊥ f (r ⊥ ) with f (r ⊥ ) being a cell-periodic function.
Integration over the x-y plane can be then accomplished by using the standard methods of the crystal field theory, 50 yielding Stokes' theorem allows the last integral over the lattice cell to be transformed into the contour integral along the cell boundaries yielding where dS ⊥ and dl ⊥ are, respectively, the elements of the Wigner-Seitz cell and its boundary, and ∂{cell} is the counterclockwise oriented contour in the x-y plane.
The symmetry of the Wigner-Seitz cell assumes that for each point at the boundary, one can put in correspondence another point at the opposite boundary such that the two points are related by a primitive lattice translation. The contributions from such points in the integral ∂{cell}ẑ × f · dl ⊥ cancel out each other since at these points, the periodic function f (r ⊥ ) takes the same values, while the contour elements dl ⊥ are oriented in opposite directions.
It follows that the second and third integrals in the right-hand side of Eq.(35) vanish. Therefore, by virtue of (30), the orthogonality relationship for the slowly varying components (18) of the field reads where the integration runs over the x-y plane. Thus, within the accuracy of the paraxial approximation, the orthogonality relationship between the wave states from different bands [with indexes q and q ′ in Eq.(36)] is not influenced by the fast longitudinal (along z axis) oscillations (14) of the electromagnetic field. Using (12), separating the integrals for longitudinal envelope functions and periodic parts [as in Eq. (32)], and noting that η qk 2L =0 with accuracy ∼ξ 2 , one obtains the orthogonality relationship for electromagnetic waves that reads where integration is over the entire cavity-unfolded structure. Note that relationship (37) is in agreement with the waveguiding theories. 41,52 With the help of orthogonality relationship (31), the solution of equations for the fast longitudinal oscillations of the field [Eq. (17)] is straightforward. As in the case of Eq. (25), it follows by applying the gauge transformation (20) and the operators −ẑ×Ê −1 andÊ −1 in the first and second equations in (17). Within the accuracy of ξ 2order terms, both equations yield the same form of the master equation for the periodic part η qk (z) of the fast longitudinal wave component where φ qk ∼ξ 2 (or less) is the diagonal matrix element qk . The off-diagonal elements have smaller magnitudes and, within the accuracy ∼ξ 2 , their contribution cannot be accounted for in the right-hand side of Eq. (38). In this case, the left-hand side of Eq. (38) is also null ( ψ (α) q ′ k |ψ (α) qk ∂η qk /∂z≡0 for q ′ = q). The integral of Eq.(38) is the expression given in (15) but now with the parameter φ qk expressed as the diagonal matrix element of effective crystal potential. This verifies the separation of variables in the paraxial gauge transformation (21).
The matrix elements φ qk provide a way for more accurate formulation of the paraxial approximation conditions in the Hamiltonian (25). Thus, the effective crystal potential in (25) has been assumed so far to be of the relative order ξ 2 , that is ch n | ln R| 2L < ∼ ξ 2 m0c 2 n 2 . Using the roundtrip self-repetition condition in a λ-cavity k z =2π/L, the condition of low-contrast patterning can be represented in the form However, as evidenced by Eq. (38), it is more important that the effective potential variations are small in average, at the scale of the crystal lattice cell. For a particular photonic band, the last condition can be expressed by evaluating the diagonal matrix element in (39), yielding where the matrix element is the same as |φ qk | in Eq. (38). The condition (40) verifies the paraxial approximation (25) within a particular set of photonic bands. For structures with simple lattice cell topology (like in Fig. 1), this condition is less restrictive. For such structures, further simplification is possible in the most important case of low-order photonic bands, which are typically characterized by relatively smooth wave functions. Introducing the reflectivity contrast parameter δR as a measure of the maximum variations in ln R(x, y), the intraband matrix elements (40) can be then estimated as | ψ (α) qk | ln R|ψ (α) qk |∼F F |δR| for the lattices with small fill factor F F < 1 2 (e.g., for the lattices with small-size pixels in the case of structures shown in Fig.1) . In the opposite case of large fill factor F F > 1 2 (large-sized pixels in Fig.1), the matrix elements are ∼(1 − F F )|δR|. It follows that the paraxial Hamiltonian (25) is valid for the structures with reflectivity patterning contrast In most practical cases, optical microcavities satisfy this condition.

III. DIELECTRIC LATTICES DEFINED BY REFRACTIVE INDEX VARIATIONS
The approach developed in Sec.II is focused on lattices of optical resonators realized by mirror reflectivity patterning (Fig.1). For completeness of the Hamiltonian (25), it has to be extended to similar lattices of coupled microcavities or parallel waveguides defined by periodic variations of the dielectric constant (Fig.2). Examples of such coupled microcavities with positive [ Fig.2(a)] or negative [ Fig.2(b)] refractive index contrast are, respectively, etched VCSEL arrays 53 or matrices of antiguided VC-SEL resonators fabricated using a two-step organometallic chemical vapor deposition (OMCVD) growth. 54 As shown below, the approach based on Hamiltonian formalism (25) applies equally well to arrays of parallel waveguides 55 [ Fig.2(a)] and low-index-core antiguides [ Fig.2(b)]. These photonic crystal materials are employed in the cladding of photonic band gap fibers 56 and photonic liquid crystal fibers. 57 In this section, the Hamiltonian (25) is thus extended to almost entire subclass of 2D paraxial photonic crystal structures characterized by the light propagation mostly normal to the periodic crystal plane. As it can be expected, the Hamiltonian for such dielectric lattices [see Eq.(52)] is analytically similar to the scalar paraxial wave equation but differs from it by the use of two-component spinor wave functions for photons.
Note that a general form of the Hamiltonian for dielectric lattices can be obtained directly from Eq. (25), assuming that R=1 and introducing the refractive index variations n(x, y) −1 − n −1 and effective mass m 0 =hk z /c n −1 in the first term of Eq. (25). These yield the effective periodic potential of U =m 0 c 2 1 n ( 1 n(x,y) − 1 n ). However, such calculus does not reveal the maximum refractive index contrast δn satisfying the conditions of paraxial Hamiltonian approximation. As shown below, the allowed refractive index variations assume the condition (59), which for low-order photonic bands, can be represented in the form where F F is the lattice cell fill factor, Λ is the lattice pitch, λ is the wavelength (in vacuum) and ξ = n −1 λ/2Λ is the order parameter (ξ≪1) corresponding to a propagation angle of state at the boundary of the first (2D) Brillouin zone. Unlike the second-order variations δR in reflectivity-patterned structures, the refractive index variations δn can thus be treated by the model as first-order perturbations [compare Eqs. (41) and (42)].
As discussed in Sec. IV, this allows the paraxial Hamiltonian to be applied to some of the structures exhibiting high contrast of refractive index variations (δn∼1), as in the case of holey photonic crystal fibers 10 or arrays of micropillars. 58 To analyze separately the impact of refractive index variations on photonic band structure, we assume that g=0 in constitutive equations (6). In the case of coupled microcavities, this suggests that the uniform loss and gain distributions are neglected and the reflectivity of the output coupling mirror is R=1. As in the case of mirror-patterned structures, the vertical cavity composition is accounted for by using an effective (complex) refractive index profile 46 n(x, y). The Bragg scattering effects set in, conditioned by the periodic pattern ε(x, y)=n(x, y) 2 /µ. For the structures considered here, the magnetic permeability µ is assumed constant (µ = 1 in the optical spectrum range). Thus, in accordance with the relationships ε(x, y) = n(x, y)/Z(x, y) and µ = n(x, y)Z(x, y), both the refractive index n(x, y) and the impedance Z(x, y) exhibit periodic variations.
The cavity-unfolded representation (Sec. II) can be applied to the matrices of dielectric microcavities as well. In this particular case, it effectively translates the cavity into a structure that exhibits a translational symmetry along cavity z-axis, allowing thus the correspondence between an array of microcavities and equivalent structure of parallel dielectric waveguides. Electromagnetic waves propagating in such waveguiding structures (in the z-axis direction) can be expressed using the paraxial wave approximation (12) with separated fast (longitudinal) and slow (lateral) wave oscillations.
Since R=1, an electromagnetic wave propagating in a cavity-unfolded structure does not show sudden amplitude variations at g-layers [see Sec. II B]. Therefore, as in the arrays of parallel waveguides as in the cavityunfolded lattices of dielectric microcavities, the fast lon-gitudinal component of the wave is a uniform plane wave showing no periodic modulation in the z-axis direction [see Eq. (14)]: The electric and magnetic field components then read After substitution of (44), Maxwell's equations for the curl of E and H yields the equations for the slowly varying wave components in the lateral xy-plane direction: In order to convert these equations into the Hamiltonian eigenproblem, an appropriate gauge transformation has to be applied for the fields e qk and h qk . A general form of the gauge transformation, which takes into account the variations of Z(x, y) and n(x, y) = µ/Z(x, y) and, at the same time, encompasses the particular case of a uniform dielectric medium (21), reads: where, as in (21), ψ qk is a two-component spinor wave function, the tensor operatorÊ is defined in (19) and the parameters a and b will be adjusted to fit the spinor transversality condition (ẑψ qk ) = 0. The corresponding inverse transformation reads whereÊ −1 is the tensor operator introduced in Eq. (24).
Note that 1/n can be substituted for Z in Eqs.(46)- (47) since for considered here structures, the magnetic permeability µ is constant. In these expressions, the spatial derivatives of Z (or 1/n) should match the paraxial approximation conditions, which thus limit the contrast of dielectric materials composing the lattice. Introducing the average impedance Z and variations δZ(x,y)=Z(x,y)− Z , one should assume that the relative variations δZ(x,y)/ Z are of the order of ξ. In that case, the gauge transformation (46)-(47) allows spatial derivatives of δZ(x, y) to be taken into account by the model: Otherwise (if δZ(x, y)/ Z ∼ξ 2 ), in conditions of the paraxial approximation, the effect of periodic variations δZ(x, y) cancels out in the gauge transformation (46)- (47), yielding thus expressions for a uniform dielectric material [Eqs. (18) and (23)]. The conditions in terms of refractive index variations δn(x, y) follow from Eq. (48) by substitution n → Z since Z(x, y)=µ/n(x, y). Following the analogy between photonic crystals and semiconductors, it is interesting to observe that the ratio between two expressions in (48) reads This corresponds to the condition 50,59,60 of smooth potential variations Λ∇δU δU ≪1 in the effective mass method widely used in solid-state physics. In the case of paraxial light propagation considered here, the condition (49) verifies the use of effective mass m 0 =hk z /c n −1 in analysis of photonic bands in the lateral propagation direction.
By substituting the gauge transformation (46) and applying the operators −ihe 3αβ Z −bÊ −1 βγ Z b and ihZ aÊ −1 βα Z −a in the first and second equations (45), both equations are converted into a similar form with respect to spinor functionsẑ×ψ qk and ψ qk : The operators applied toẑ×ψ qk and ψ qk (in the lefthand side) show different longitudinal components (third term in each equation). At the same time, the gauge transformation (46) implies transversality of spinor functions (ẑψ qk =0), such that these longitudinal terms cancel out by adjusting the gauge transformation parameters a and b : where µ is the magnetic permeability (a= 3 2 and b=− 1 2 in the optical spectrum range).
The vector cross productẑ×ψ qk in the first equation (50) is the result of spin operator 20ŝ z =iẑ× applied to the spinor ψ qk . In the paraxial gauge transformation considered here (withẑψ=0), the spinor functions are invariant under the operatorŝ 2 z . Therefore, substituting parameters (51) and taking a cross product ofẑ and first equation, one obtains in (50) two identical eigenproblems that read where m 0 =hk z /c n −1 is the effective mass andp ⊥ = −ih∇ ⊥ is the momentum operator in the lateral (xyplane) direction. The third term in the Hamiltonian (52) is the effective crystal potential induced by variations of the (complex) refractive index n(x, y). This is the stationary Schrödinger equation for photons in dissipative dielectric lattices of coupled microcavities or parallel waveguides. The orthogonality relationship can be established between the concomitant partners of biorthonormal set of its solutions [see Eq.
In the particular case of parallel dielectric waveguides, the conventional waveguiding theories have established the mode orthogonality relationship in terms of the electric and magnetic field components. The biorthonormal set of solutions in mirror-patterned structures is shown to satisfy such relationship [Sec. II, Eq. (37)]. However, the paraxial gauge transformation in periodic dielectric structures [Eq. (46)] differs from the one in mirrorpatterned microcavities with uniform dielectric material in the cavity [Eq. (18)]. Therefore, it is crucial to verify that the gauge transformation (46) and the biorthonormal orthogonality relationship (53) are in agreement with the well-established results of conventional waveguiding theories.
An electromagnetic wave associated with the concomitant partnerψ qk [in (53)] propagates in the same structure as the wave (46) but with the complex conjugated refractive index and impedance: Respectively, the inverse transformation reads Substituting −ẑ×[ẑ×ψ q ′ k ′ ] [Eq. (55)] and ψ qk [Eq. (47)] in the orthogonality relationship (53), we obtain The second and third terms in the right-hand side of this equation are null. The first term coincides with expression in Eq. (34), yielding the orthogonality relationship [see Eq. (36)] The fieldsẼ q ′ k ′ andH q ′ k ′ of the wave associated with the concomitant partnerψ q ′ k ′ are defined by relation-ship (44) with the complex conjugated dielectric function. Taking this fact into account, one can verify that where integration runs over the entire structure of parallel dielectric waveguides (or cavity-unfolded array of microcavities.) The expression (58) corresponds to conventional orthogonality relationship between the modes of parallel dielectric waveguides. This result verifies the theoretical treatment presented in this paper.
For low-order photonic bands, the structure parameters satisfying the paraxial approximation conditions can be defined more precisely. Thus, for a particular photonic band, the conditions (48) expressed in terms of intraband matrix elements read where the relationship Z(x, y)=µ/n(x, y) is taken into account. Being applied to the low-order bands, these expressions yield the relationships (42) [see also the discussion in Sec. II C].

A. Generalized Hamiltonian
Combining the results of Sec. II and III [Eqs. (25) and (52)], we obtain the Hamiltonian for entire subclass of 2D structures characterized by paraxial light propagation in the direction normal to periodic lattice plane The photonic state wave function ψ qk (spinor) is related with the corresponding electromagnetic wave via paraxial gauge transformation where operatorÊ was introduced in Eq. (19) and parameters a and b were obtained in Eq. (51). The effective crystal potential is defined by the third and fourth terms in the left-hand side of Eq. (60). In the case of 2D arrays of coupled microcavities, it takes into account the effects of mirror reflectivity patterning, as in metal-patterned VCSEL arrays 11 (Fig. 1), and dielectric material variations, as in periodically etched VCSEL structures 53,54 (Fig. 2). Such simple expression for the effective potential was obtained by unfolding the cavities along the optical axis [z-axis in Fig. 1(d)] and representing the standing optical modes in the form of propagating Bloch waves (61) in equivalent 3D structure.
In a periodic array of microcavities, the longitudinal component of wave vector (k z ) is fixed by the selfrepetition condition at the cavity roundtrip. At the same time, Eq. (60) cannot reproduce the cavity resonance condition. Since Eq. (60) takes into account a phase shift at the cavity mirrors, k z has to be evaluated from the roundtrip condition in a cavity with perfectly reflecting mirrors. In particular, k z =2π/L in the case of one-wavelength microcavities (L=λ/n). In addition, analyzing a symmetry of the group of k, the z-axis nonreciprocity of equivalent cavity-unfolded structure has to be taken into account (see Sec. II). Obviously, for an array of parallel dielectric waveguides or antiguides (e.g., photonic band gap fibers, 56 photonic liquid crystal fibers 57 ), these restrictions of the model do not apply. The effective potential is uniquely defined by the refractive index variations [fourth term in the Hamiltonian of Eq. (60)].
For most important (in practical applications) loworder bands, and within the range of parameters limited by conditions (41) and (42), the Hamiltonian is suitable for structures with high contrast of refractive index variations (e.g., etched arrays of pillar microcavities, 58 holey photonic crystal fibers 10 ). The class of 2D photonic crystal materials encompassing valid solutions of the non-Hermitian Hamiltonian eigenproblem (60) is further illustrated below with several structure examples.
Arrays of semiconductor microcavities (n=3.5) defined by mirror reflectivity patterning (Fig.1) and operating at the optical wavelength λ∼1 µm typically employ lattices of the pitch Λ∼5 µm. These parameters assume that ξ 2 ∼10 −3 [Eq. (41)] and for any lattice fill factor F F , the Hamiltonian can be applied to the structures with mirror reflectivity contrast δR up to 10 −2 . Note that in a typical VCSEL array with mirror reflectivity patterning, the reflectivity contrast is in the range from ∼10 −4 to ∼10 −3 .
The first working silica photonic crystal fiber 10 has a cladding material consisting of 300 nm air holes arranged in a hexagonal lattice of 2.5 µm pitch. The refractive index of fused silica varies from 1.55 to 1.44 in the wavelength range of 0.2 − 1.5 µm. The lattice fill factor and the contrast are thus F F ∼0.02 and |δn|∼0.5.
In the case of arrays of parallel dielectric waveguides (or antiguides) with low-contrast refractive index variations (|δn|≪1), the Hamiltonian (60) is valid for any fill factor of the lattice, in the optical wavelength range from λ∼4|δn|Λ [Eq. (42)] to λ∼0.2nΛ. (The long wavelength range is limited by the paraxial approximation condition ξ≪1 with the critical value of ξ∼0.1 .) The silica photonic band gap fibers reported in Ref. 56 (n∼1.5, δn=0.015 and Λ=6 µm) satisfy the paraxial approximation conditions in the visible (λ>360 nm) and near infrared regions of the optical spectrum, up to the upper transparency edge of fused silica (λ∼1.5 µm).
In the case of 2D photonic lattices exhibiting inversion symmetry (e.g., square or triangular lattices), the biorthonormal system of lattice-periodic functions can be readily obtained by applying the P T transformation (time reversal followed by coordinate inversion) to Eq.(60) and noting the degeneracyhω qk↑ =hω qk↓ . This has much in common with the Kramers degeneracy 62,63 of a single-electron Hamiltonian in lattices with inversion symmetry. The invariance of the single-electron Hamiltonian under the P T transformation yields degeneracy of eigenvalues E qk↑ = E qk↓ associated with the states of opposite spin.
In the Hamiltonian (60), the degeneracy on spin variable (hω qk↑ =hω qk↓ ) is caused by the fact that the Hamiltonian contains no spin operators. On the other hand, the Hamiltonian itself is not invariant under the P T transformation. In the lattices with inversion symmetry, its P T transform is uniquely defined by the result of the time reversal operation (complex conjugation followed by substitution −t → t), 43,49 which transformsĤ intoĤ * =Ĥ + .
In (66), the relative phase of co-partners ψ qk↑,↓ and ψ qk↑,↓ is set by the orthogonality relationship (31). This relationship evidences that an arbitrary phase factor e iα can also be introduced in both equations (66), modifying thus the absolute phases of co-partners but preserving their relative phase shift. The fact that the co-partner phases have to be the same (at t=0) can be seen from the following observation as well: There should be no difference between the co-partners in the case of Hermitian Hamiltonian. Therefore, the functions ψ qk↑,↓ (t, r ⊥ ) and ψ qk↑,↓ (t, r ⊥ ) coincide in the limit Im(ω qk )→0.
Eq. (66) defines the biorthonormal system of solutions in the case of lattices with inversion symmetry. It effectively transforms (60) into a scalar Hamiltonian eigenproblem with respect to the amplitude ψ qk (t, r ⊥ ) of nonzero spinor component. The orthogonality relationship for the biorthonormal set of scalar functions ψ qk (t, r ⊥ ) andψ qk (t,r ⊥ )=ψ * qk (−t,−r ⊥ ) follows from Eq. (31): Finally, note that in the case of significant difference between reflection coefficients of a mirror for s-and p-polarized waves in microcavities ( | ln R s /R p | > ∼ ξ 2 in (1)) or large polarization anisotropy in periodic dielectric lattices (|n x −n y | n −1 > ∼ ξ 2 ), the spin-orbit coupling effects set in, rendering invalid the scalar approximation (66). Another example of spin degeneracy removal, for which the scalar approximation (66) is unsuitable, is the Coriolis-Zeeman splitting of photonic energy bands in nonpermanent gravitational field. 20

C. Biorthonormal plane wave expansion in square lattices
For square lattice structures (Sec. IV D), it is more convenient to express the biorthonormal set of spinor functions (66) in the Cartesian coordinates basis, in the form of functions ψ= ψ (x) ψ (y) . For a first rank tensor ψ representing a state of spin s=1, the relationships 34 ψ (x) =i(ψ 1,1 −ψ 1,−1 )/ √ 2 and ψ (y) =(ψ 1,1 +ψ 1,−1 )/ √ 2 provide a transformation between the Cartesian coordinates and |s, m s functions bases. (For paraxial photonic states considered here, the spinor component ψ (z) =−iψ 1,0 is null and therefore not indicated explicitly.) The degenerate spinor functions (66) can be represented as linear combinations of positive-and negativespin states yielding the biorthonormal set of solutions associated with the x-and y-polarized states of electromagnetic field Here, the spinor components are expressed in the basis of Cartesian coordinates. Like the functions (66), degenerate spinor functions ψ qk,x and ψ qk,ŷ (hω qk,x =hω qk,ŷ =hω qk↑,↓ ) convert Eq. (60) into an eigenproblem with respect to scalar amplitudes ψ qk (t, r ⊥ ), which assume the orthogonality relationship (67).
For nonzero spinor components, the orthogonally polarized states (69) show equal distributions |ψ qk (t, r ⊥ )| 2 associated with the energy flux along the z axis. In experiment, such states are observed as orthogonally polarized modes showing indistinguishable intensity patterns. 22,23 The stationary Schrödinger equation (60) is solved here using a biorthonormal plane wave expansion of the periodic crystal potential and Bloch functions (69). In the stationary case, the time evolution of wave functions can be omitted, yielding the biorthonormal set of scalar amplitudes ψ qk (r ⊥ ) andψ qk (r ⊥ )=ψ * qk (−r ⊥ ). The relationship between co-partners of the set implies complex con-jugated coefficients of expansion: where G is reciprocal lattice vector in the x-y plane. The difference in expansion coefficients in non-Hermitian and Hermitian Hamiltonian cases can be appreciated by examining the Parseval theorem for coefficients C qk (G): It can be seen that G C 2 qk (G)=1, as opposed to the usual expression G |C qk (G)| 2 =1 (Parseval's theorem). Note that for the lattices with inversion symmetry discussed here, the two sums are consistent in the limit Im(ω qk )→0 (Hermitian Hamiltonian case). In this limit, ψ qk →ψ qk and the phases of wave functions can be adjusted to obtain real expansion coefficients C qk (G).
Thus, in photonic structures exhibiting inversion symmetry of the lattice, the biorthonormal plane wave expansion series differ from the usual OPW series by normalization condition for the amplitudes of spatial harmonics. The Schrödinger equation (60) can be then readily converted into a matrix equation for C qk (G) coefficients. Furthermore, the inversion symmetry of the lattice assumes a simple relationship between the the matrix elements of operatorsĤ andĤ + . Thus, the periodic crystal potential and its adjoint operator are represented by series of lattice harmonics with complex conjugated coefficients where V −G =V G due to the symmetry of crystal potential U (−r ⊥ ) = U (r ⊥ ). Substituting (70) and (72) into the Schrödinger equation (60) and its P T -transform, multiplying by exp(−i(k ⊥ + G ′ )r ⊥ ) and integrating over a lattice cell in the x-y plane, one obtains the matrix equations These are the two complex conjugated matrix eigenproblem equations. Within the accuracy of complex eigen values and normalization condition for expansion coefficients, the eigenproblem (60) is thus converted into the usual form encountering in conventional OPW expansion method.
D. Band structure of square-lattice arrays of optical microcavities or parallel waveguides To study the light propagation behaviour in paraxial photonic structures incorporating loss and gain distributions, the non-Hermitian Hamiltonian (60) is used here to analyze the structures with simple lattice symmetry and cell topology. Its numerical solutions are reported for square-symmetry lattices depicted in Figs.1 and 2. The results apply both to the structures defined by mirror reflectivity patterning (e.g., metal-patterned VCSEL arrays 11 ) and to the ones exhibiting periodic variations of refractive index (e.g., etched VCSEL arrays 53 or matrices of antiguided VCSEL microcavities 54 ).
All matrices of microcavities considered here have similar lattice cell topology indicated in Figs.1 and 2. In the case of cavities with mirror reflectivity patterning, the position of optical microresonators is defined by highreflectivity square pixels separated by low-reflecting cavity mirror domains forming a grid pattern. The fill factor of such square lattice is the area ratio F F =a 2 /Λ 2 , with a being the square pixel width and Λ being the lattice pitch. In similar dielectric lattices defined by refractive index variations, the position of microcavities is set by square dielectric waveguide (antiguide) cores embedded in the background of the cladding material. Obviously, the same expression for the fill factor can be used to characterize such lattices of dielectric waveguides, with a being the waveguide (antiguide) core width.
The Hamiltonian (60) and the gauge transformation (61) utilize equivalent, cavity-unfolded (3D) representation of microcavities. As discussed in Sec. IV A, the longitudinal component of propagation vector k is defined by condition k z =2π/L (in case of one-wavelength microcavities). Furthermore, by virtue of the dissipative effects rendering the opposite z-axis directions nonequivalent, the group of k in Eq. (61) contains only symmetry operations preserving the z-axis direction.
In the particular case of square lattice, the cavityunfolded 3D structure (2L-periodic in the z-axis direction) has a tetragonal symmetry (Λ =2L) with symmorphic space group Γ q associated with D 4h point subgroup of rotations. 49,64,65 The reciprocal lattice is of the tetragonal symmetry as well with the first Brillouin zone (BZ) being of the rectangular prism shape [ Fig.1(e)]. However, the crystal z-axis nonreciprocity implies that rotations and reflections of the D 4h group altering the zaxis direction are not allowed. Therefore, for square arrays of microcavities, the group of k has a reduced symmetry characterized by C 4v point group rotations in the ∆ and T points of the BZ [at k=(0, 0, k z ) and (± π Λ , ± π Λ , k z ), respectively] and C 2v symmetry in the Z point [at k = (± π Λ , 0, k z ) or (0, ± π Λ , k z )]. Figs.3-5 show the results of band structure computations for square lattices of microcavities with the pitch Λ = 5 µm and lattice cell fill factor F F = 0.5. The photon energy dispersion and lifetime energy broadening curves along the high-symmetry lines ∆-Z-T -∆ in the BZ are obtained here using 25 lattice harmonics in Eq.(73). In numerical simulations, the cavity structure is assumed to be optimized for 960 nm wavelength (operation wavelength range of GaAs/AlGaAs VCSEL structures with InGaAs quantum wells in the optical gain region).
Figs. 3 and 4 detail the impact of mirror patterning contrast δ ln R = δ|R| + iδϕ and take contributions of the periodic phase (δϕ) and amplitude (δ|R|) variations of the mirror reflectivity into account. Fig. 3 shows the energy dispersion curves for low-oder photonic bands in the case of a lattice defined by phase modulation of mirror reflectivity with parameters δ|R|=0 and δϕ=10 −2 (black curves). The reflectivity at pixel positions is R = 1 while the periodic pattern R(r ⊥ ) is introduced by a phase shift −δϕ at reflections form the cavity mirror domains forming a grid. Comparison with the band structure of empty lattice (Fig.3, gray curves) indicates a blueshift of the optical modes due to contribution of phase-advancing domains of the grid. The main outof-phase mode T 5 (doubly degenerate by polarization) 23 has maxima of probability amplitude |ψ| located at the pixel positions and zeros located at the grid. Therefore, it shows the smallest energy shift. The contribution of the grid is particularly pronounced for T ′ 5 states exhibiting |ψ| distribution with maxima located at the cross points of the grid and zeros at the pixels.
Numerically calculated wave functions of these states (not shown in the figures) are in good quantitative agreement with the results 20 obtained by means of group the- oretical analysis. ? The intensity patterns of the optical modes associated with these states are considered in Refs. 21,22 both theoretically and experimentally.
Even for such low-contrast reflectivity patterns as the ones considered here, a complete 2D band gap can be opened in the lateral direction, between the T 5 and Z 2 states. Thus, in Fig.3, the gap in the spectrum of optical modes is of 0.1 meV width. Since the effect of the grid on T 5 states is small, the width of the energy gap is mostly defined by the blueshift of Z 2 states of an empty lattice.
Wave functions of the Z 2 and Z 3 states in the first BZ have zeros only along one crystalline direction (either x or y-axis directions) and differ by the position of probability amplitude |ψ| maxima. Along this lattice direction, the Z 3 states are localized to the pixels while Z 2 states are localized to the grid, like, respectively, T 5 and T ′ 5 states. However, wave functions of the states originating from the Z points in the first BZ show large probability amplitude oscillations (and localization) along one lattice direction, as opposed to the states from the T point with wave functions oscillating in both lateral crystalline directions. For completeness of the discussion note that the wave functions of the lowest energy states in the ∆ point (doubly degenerate by polarization ∆ 1 states) show no zeros of probability amplitude |ψ|. Therefore, due to the effect of the grid pattern, the Z 3 (Z 2 ) states exhibit an intermediate blueshift as compared with the energy shifts of the ∆ 1 and T 5 states (∆ 1 and T ′ 5 states, respectively) of an empty lattice. Fig. 3 shows other states in the BZ (Z ′ 3 states) exhibiting small energy blueshift, which is comparable to that one of the T 5 states with wave functions localized to the pixel positions. The Z ′ 3 states originate from the next nearest equivalent Z-points of reciprocal lattice at (± π Λ , ±2 π Λ , k z ) and (±2 π Λ , ± π Λ , k z ). Therefore, in addition to the unidirectional oscillation features of the Zstates in the first BZ, their wave functions show large oscillations of probability amplitudes in the second crystalline direction as well. The probability densities of the Z ′ 3 states thus have a better overlap with the pixels of reflectivity pattern, which explains the smaller blueshift energy of Z ′ 3 states as compared to their counterparts in the first BZ. In Fig. 3, the energy broadening of the bands due to the optical cavity loss has not been taken into account. Fig.4 illustrates the effects of optical loss distributions in the structures with the same lattice parameters as in Fig.3 (Λ = 5 µm and F F =0.5). All structures in Fig. 4 have the same amplitude reflectivity patterning |R(r ⊥ )| and differ only by the phase contrast of the pattern. The amplitude reflectivity |R(r ⊥ )| is 0.999 and 0.989 for the pixels and grid domains of the pattern, respectively, yielding the amplitude contrast of δ|R|=10 −2 . Note that the pixel reflectivity corresponds to the cavity loss of 0.2%, which is in the range of optical losses in a typical VCSEL structure. As in the case of the structures in Fig. 3, in Fig. 4, reflections at pixel domains of reflectivity pattern introduce no additional phase shift into the cavity roundtrip phase accrual of optical modes. Variations in the photon energy and lifetime broadening dispersion curves in Fig. 4 are thus introduced by different phase-advancing shifts at reflections from the grid domains in these structures.
At no phase modulation of reflectivity pattern [δϕ = 0], the energy band structure is close to the one of an empty lattice but also shows a set of new peculiar features [ Figs.4 (a), black curves]. Thus the degeneracy in energy of photonic states in high symmetry points of the BZ is partially removed (the states are doubly degenerate by polarization) and a set of partial flat bands appears in the Z-and T -points of the BZ. These features of energy dispersion curves are uniquely defined by the dissipative effects in the structure.
The energy broadening curves (or loss-dispersion curves) of optical modes in this structure are shown in Fig. 4 (b). (Only the losses related with the disperssive features of the patterned cavity mirror are taken into account.) The doubly degenerate photonic states T 5 have the lowest cavity loss (and lifetime energy broadening). They are of particular interest since they define the main lasing modes in coupled laser arrays. 22,23 The Z ′ 3 states define the next lowest-loss modes in the high-symmetry points of reciprocal lattice. Finally, the states T ′ 5 are associated with the highest-loss optical modes. Using the same considerations as in the case of energy dispersion curves in Fig. 3, one can readily explain these features of the loss dispersion curves in terms of the overlap between photonic state wave functions and high reflectivity pixels.
Introduction of phase variations at the array grid does not affect the optical losses of modes associated with T 5 , Z ′ 3 and T ′ 5 states [see Figs.4 (b), (c) and (d)]. Fig.4 shows the energy-and loss-dispersion curves calculated for the phase variations contrast δϕ of 10 −1 (gray curves) and 10 −2 (light gray curves). Note that in the last case δϕ=δ|R|. For small phase variations δϕ<δ|R| (gray curves), the structure of energy dispersion bands is close to that one in the case of pure amplitude modulation of reflectivity pattern (see Fig.4 (a), gray and black curves). However, for δϕ∼δ|R| (light gray curves), it approaches the energy band structure in the case of pure phase modulation of mirror reflectivity (Fig.3, black curves) and exhibits a band gap between T 5 and Z 2 states.
A common feature of the loss-dispersion curves in Fig. 4 is the opened photonic band gap in the loss domain 9 (or in the domain of photon lifetime in the cavity). Thus, it is impossible to excite an optical mode showing longer cavity lifetime than in the T 5 state. By properly varying the lattice cell fill factor, 9,15 this effect has been used to implement photonic crystal heterostructure wells capable to confine photonic envelope wave functions to the regions of lower band gap material. The confined states show usual features with cosine envelope functions in the well and exponential tails in the cladding material. The dissipative photonic crystal materials have been used so far in such photonic crystal heterostructures show no band gaps in the photon energy domain.
For a structure with complex parameter of reflectivity patterning contrast δϕ=δ|R|=10 −2 (the same parameters as in Fig. 4, light gray curves), Fig. 5 shows variations of the two photonic band gaps (in the photon energy domain and in the optical loss spectrum) as a function of the lattice cell fill factor F F . The lowest loss state T 5 (gray curve, right axis) defines the band edge in the cavity loss domain. The gap below this edge is opened at any fill factor of the lattice. In the range of lattice cell fill factor 0.16−0.66, the energy of Z 2 state is higher than the band edge T 5 , such that a second band gap originating at the band edge T 5 exists in the photon energy domain, in parallel with the gap in the photon lifetime domain.
The notion of double photonic crystal band gap illustrated in Fig. 5 opens new possibilities for tailoring photonic envelope wave function and controlling quantization features of confined photonic states in photonic crystal heterostructures. Thus, the eigenvalues of the Hamiltonian (60) at a photonic crystal heterostructure barrier assume the dispersion relationship where indexes A and B distinguish photonic crystal materials at the heterostructure barrier and the parabolic band approximation is used for both photonic crystal materials. For complex band edge parametershω A,B and effective masses m A,B , this condition assumes that both propagation constants k ⊥,A and k ⊥,B are complex, independent of particular photonic crystal heterostructure configuration. Thus, for an N -dimensional photonic crystal heterostructure well, Eq.(74) envisages a possibility of 2Ndimensional confinement of photonic envelope wave functions by introducing quantization of both real and imaginary parts of propagation vector k ⊥ = k ′ ⊥ + ik ′′ ⊥ of confined photonic states. For such states, k ′′ ⊥ = 0 even in the region of lower band gap material (at the well core), allowing the confined states to be excited at the energies within forbidden energy gaps of the heterostructure materials (well core and barrier materials).
In the numerical examples of band structure computations presented here, photonic structures utilizing reflectivity patterning for definition of periodic crystal lattice are considered. The effective crystal potential of these structures is governed by the third term in the Hamiltonian of Eq. (60) and the magnitudes of its matrix elements are bounded to the second-oder perturbations (of the relative order ∼ ξ 2 ). Therefore, in considered here case of paraxial light propagation with photon energy of 1.3 eV and ξ 2 ∼10 3 , the Hamiltonian (60) yields accurate estimates of photonic bands splitting up to 1 meV . At the same time, for dielectric lattices defined by periodic variations of refractive index, the effective crystal potential [fourth term in the Hamiltonian] is bounded to the first-order perturbations (∼ ξ). In the case considered here (hω = 1.3 eV , ξ∼0.03), the Hamiltonian (60) allows energy dispersion curves with band structure splitting up to 40 meV to be analyzed.

V. CONCLUSION
In this paper, a simple non-Hermitian Hamiltonian formalism is developed for a subclass of two dimensional photonic crystal structures characterized by paraxial light propagation in the direction normal to the lattice plane, like in the case of coupled microcavity arrays or microstructured photonic crystal fibers. It allows the optical loss or gain distributions to be taken into account in the band structure analysis and envisages the effect of double photonic band gap opened both in the photon energy and lifetime domains. Predicted novel features of optical mode behaviour at double photonic crystal heterostructure barriers with band edge discontinuities in the energy and lifetime domains offer new possibilities for photonic crystal applications in optoelectronic devices and integrated photonic circuits.