Decoupling the energy and momentum of photons in the quasistatic limit
M J A Smith P Y Chen
physics.optics


APS/123-QED

 


 School of Mathematics, The University of Manchester,  Manchester M13 9PL, United Kingdom 
 michael.j.smith@manchester.ac.uk
 

Unit of Electro-optic Engineering, Faculty of Engineering Sciences, Ben-Gurion University, Beer Sheva, Israel

 School of Physics and Astronomy, Raymond and Beverly Sackler Faculty of Exact Sciences, Tel Aviv University, Tel Aviv, Israel

 

 

 


We  theoretically show that  the  frequency and momentum of a photon are not necessarily proportional to one another at low frequencies in     photonic crystals  comprising materials with     positive- and negative-valued   material properties.    We rigorously determine    closed-form  conditions   for  the light cone  to emanate from points other than the origin of k space, ultimately decoupling the first band from the origin and demonstrating light propagation at zero energy with nonzero crystal momentum.   We also numerically show that  first bands   can   originate from an arbitrary Bloch coordinate as well as from multiple   coordinates simultaneously.  


 
Violating the Energy-Momentum Proportionality of Photonic Crystals in the Low-Frequency Limit
    Parry Y. Chen
    December 30, 2023
=============================================================================================

When a photon propagates through a dielectric medium at low frequencies, it satisfies the energy-momentum (E-k) relation E = c ħ ||, where c is the phase velocity in the medium <cit.>. This relation ensures that at zero energy, the photon possesses zero momentum. Fundamental relations of this type are prevalent throughout nature and are not isolated to photons, for example, electrons propagate through a crystal lattice as E = ħ^2 /(2 m_eff) ||^2 at low energies, where m_eff is the effective mass <cit.>.   This proportionality is fundamental for the study of particles and fields in relativistic mechanics, particle physics, and quantum mechanics.   


In this Letter, we break the conventional  low-frequency E-k proportionality for photons,   obtaining     relations of the form E = C  | - Ξ|, where Ξ denotes a high-symmetry point of the reciprocal lattice and  is   the crystal momentum (see Figs. <ref>-c).
This is achieved in  two-dimensional  photonic crystals comprising materials with positive-definite and negative-definite <cit.> optical properties.  We present explicit    conditions on the constituent properties and  explicit forms for the proportionality constants C,   for all high-symmetry coordinates of a   square lattice. 
Furthermore, we numerically demonstrate the existence of other novel low-frequency behaviors, including   photonic crystals with E = C_1   || + C_2  | - |, where C_j are proportionality constants and  X is a high-symmetry point (see Fig. <ref>). Such  unconventional behavior contrasts the   standard outcomes for light in photonic crystals, where  either E-k proportionality is supported,  or there exists a complete   band gap <cit.>,  at low frequencies.
 
 
In the  nonstandard  photonic settings we describe,   massless photons are predicted to  propagate as massive polaritons which  travel superfluidically through the medium   <cit.>. Consequently, our findings have the potential to  motivate the  development of new photonic devices, and to deepen our understanding of light in structured  media. The behaviors   we describe     complement  existing observations in optical systems   incorporating negative-definite materials,  such as  folded band surfaces with infinite group velocities <cit.>, cloaking and superresolution  <cit.>, and   new  types of band gaps <cit.>.   Analogies to our  low-frequency E-k relations  may be found  in the electronic properties of transition-metal perovskites,  where the first   band is centered about   high symmetry points other than the origin <cit.>.   
 
 
 


 
We begin by considering the modes of  the time-harmonic form  of the source-free Maxwell equations in a   nondispersive and lossless system with  Bloch vectors   = (k_x, k_y,0).  This    wave vector restriction reduces Maxwell's equations to    the   Helmholtz equation
 
    ∇_⊥·( ε_^-1∇_⊥ H_) +  ^2 c_0^-2μ_ H_ = 0,

 for   fields polarized as   = (0,0,H_). Here, ∇_⊥≡(∂_,∂_), ε_   is the relative permittivity, μ_ the relative permeability,    the angular   frequency, and c_0 is the speed of light in vacuum.   We consider an array of infinitely extending isotropic cylinders, periodically positioned in the (x,y) plane at the coordinates of a square lattice, and which are embedded in an infinitely extending isotropic background material. In the background and cylinder domains,   material constants   are allowed to be negative-valued.  At the  cylinder edges  we impose      continuity conditions, and     between unit cells we impose   Bloch–Floquet conditions.  This admits the system <cit.>
 
    _l  ℬ_l + ∑_m=-∞^∞ (-1)^l+m S_m-l^ℬ_m = 0,

  where S_m^ = S_m^(,_;n_) denotes lattice sums (see   Supplemental Material <cit.>mcphedran1996low,poulton2000eigenvalue,mcphedran2000lattice,linton1998greens,twersky1961elementary,chin1994greens,bensoussan1978asymptotic,jikov2012homogenization,bergman1979dielectric,mcphedran1980electrostatic,busch1998photonic,linton2010lattice,abramowitz1964handbook,nicorovici1996analytical),  = /c_0,     _ = (k_x, k_y), ℬ_m are     amplitudes of the cylindrical-harmonic basis functions, and N_m   are   inverse cylindrical-Mie coefficients. Where applicable,   subscripts  and  denote  the  background and cylinder properties, respectively.  The dispersion equation  for the crystal  is given by the vanishing determinant   of (<ref>), which we     truncate to dipolar order.
 
   
   
  We now outline the procedure for determining when a  low-frequency band surface     emerges from the Γ point. First, we    evaluate   expansions  for  N_m in ω_ (these     are lengthy, see  Supplemental Material <cit.>).  Next, we determine   closed-form expressions    for   the S_l^   in  (<ref>) at low frequencies and about the Γ point; these are obtained following <cit.> (also extensive, see Supplemental Material <cit.>). Assuming that   = α k_, where α is real and positive-valued, we subsequently obtain series coefficients for S_l^ in  alone.  Substituting  the  expansions for N_m and S_l^ into (<ref>), the zero determinant condition is satisfied to the lowest order  
 
 for  α such that 
 
    = {1/ε_( 1+f  τ/ 1-f τ)
       1 /μ_+f    (μ_-μ_) }^1/2 k_,

 where τ = ( ε_ -   ε_ )/( ε_ +   ε_ ), f = π a^' 2/a^2 is the filling fraction, a^' denotes the radius of the cylinders, and a is the lattice period. Thus,  the constituent permittivity and permeability values can be negative, but provided α>0 then a band surface will emerge from Γ.   Here α = 1/n_eff where n_eff   is the effective refractive index <cit.>.  
    
   
   
 
 
    
   
   
   


   Next, we determine the conditions   and asymptotic behavior   of a   band that emerges from the M point at low frequencies. As before, we derive  asymptotic forms for    the  S_l^  sums near M (see   Supplemental Material <cit.>) and assume  = α^' k_^', where α^' is real and positive-valued,  _^' = _ - and = (π/a,π/a). This assumption yields series coefficients for  S_l^  in    (given in Supplemental Material <cit.>).   Substituting  these expansions  for   S_l^ and N_m  into     (<ref>),   the  zero determinant condition is satisfied to the lowest orders   when  

    ε_ = -ε_,

 and for α^'   such that
  
    = {1/  8 π^ 2| 64 π^4 ^4 θ_^' + Γ(14 )^8 |^1/2[ ε_μ_ + 2 ε_μ_log( 16   π^2/f   Γ(14)^4) ]^-1/2. 
    .  ×[ ε_μ_ + 2 ε_μ_log( 16   π^2/f   Γ(14)^4) ]^-1/2} k_^',

where Γ(z) is the Gamma function.    That is, a band surface is supported from   M   at   low frequencies   provided    (<ref>)    and  α^'>0 
  are satisfied.  The condition ε_ = -ε_    corresponds to an anomalous resonance   in      quasistatic problems <cit.>   (discussed below). 


Likewise,  for the X point at low frequencies, having derived       asymptotic  forms for the   S_l^ sums near X (see Supplemental Material <cit.>)  we assume that = α^'' k_^, where α^'' is real and positive-valued,  _^ = _ -, and = (π/a,0).     Substituting    the resulting   expansions   for   S_l^, and N_m,    into (<ref>),    the   zero determinant condition is satisfied to leading order      when  

    ε_ = ( ζ - 16   π^2 /ζ+ 16  π^2) ε_,

where ζ =   Γ(14)^4  f. At   the next order, provided

    μ_ = ( ζ ( ζ - 64     π^2 )  
     + 512   π^4 log( ζ/(32   π^2) )   /(ζ - 16 π^2)^2)μ_,

  then we obtain  the low-frequency dispersion relation


    = {( 16 Γ(14)^4 π^2+ 64 π^4 ^2θ_^ -Γ(14)^8 ^-2θ_^) /16 Γ(14)^4 π^2  ε_μ_}^1/2 k_^.

However,   the slope  in (<ref>) is only real-valued  along   Γ X and XM. Numerical investigations confirm     elliptical contours at low frequencies; interpolating between these   paths with the ansatz     ^2 = α_^ k_^ 2 + α_^ k_^ 2   we   obtain  
 
    ^2 = (  16 Γ(14)^4 π^2 + 64 π^4 - Γ(14)^8 /16 Γ(14)^4  π^2 ε_μ_) k_^ 2
    
    + (  16 Γ(14)^4 π^2 - 64 π^4 + Γ(14)^8 /16 Γ(14)^4  π^2 ε_μ_)  k_^ 2,

  for the first band surface as → 0 and  as _→ X. 

Hence, a band surface is supported from X at low frequencies provided  (<ref>),   (<ref>),  and ε_μ_ >0   are satisfied.  In  (<ref>), the proportionality factor  is      negative-valued for f≲ 0.914, and the proportionality factor in  (<ref>) is negative-valued for all f, demonstrating that highly restrictive sign-changing conditions must be satisfied in both ε_ and μ_ so that  the first band   emerges from X.  The number of conditions for each high-symmetry Bloch coordinate is entirely due to the different asymptotic behaviors of   S_l^.   For arbitrary Bloch coordinate origin, we anticipate that the number of conditions  will change significantly.
  
     


We now compare the  asymptotic forms  above  against results from a fully numerical treatment of (<ref>) within the dipole truncation.    We begin by validating   (<ref>) for a regular photonic crystal; in Fig. <ref> we present the band diagram of a  representative   crystal  with    ε_ = 1 and μ_ = 2 embedded in   air (ε_ = μ_ = 1). As expected, the first band   emanates from  the Γ point, and  (<ref>) shows excellent agreement.  In Fig. <ref> we consider  ε_ = -1 and μ_ = 2, where  the first band emanates from the M point as described by (<ref>) at low frequencies, also with excellent agreement.   In Fig. <ref> we consider  ε_≈ -0.53 and μ_≈ -10.22  satisfying (<ref>) and  (<ref>). Here, the  slope differs along   Γ X and X M, demonstrating   twofold symmetry  as → 0. The   asymptotic estimate  (<ref>)    shows  excellent agreement near X at low frequencies. The   diagram also possesses folded bands <cit.> at high frequencies.
  
  
 In Fig. <ref> we present     iso-frequency contours for the first band surfaces of    the   photonic crystals considered in Fig. <ref>, over a quarter   of the first Brillouin zone. In Fig. <ref> we see  a cone  (∞-symmetric) as → 0,   whereas in Fig. <ref> we observe fourfold symmetric contours, as expected from the exp(4θ_^') dependence in (<ref>).   In Fig. <ref> we see that the first    band  has twofold symmetric contours,   as expected from (<ref>). These low-frequency symmetries contrast with the electronic band diagrams of graphene (and photonic analogues to graphene)  where the high-energy band emanates from    the K point   and about the Fermi energy E_F as an ideal cone    <cit.>. 
  
   
Having   numerically validated the  E-k relations (<ref>), (<ref>), and (<ref>), we now briefly demonstrate that magnetic constituents are not necessary to observe exotic low-frequency behaviors. Low-frequency descriptions of nonmagnetic photonic crystals   emanating from Γ  and M are obtained by the replacements μ_, μ_↦ 1 in  (<ref>) and (<ref>)    above. This is despite the fact that the N_m coefficients for  nonmagnetic crystals exhibit  different leading order behavior for small  (see Supplemental Material <cit.>).  However,   the new leading-order behavior of N_m   yields a nonmagnetic analogue to (<ref>)  of the form

    ζ + 8π^2 - 16 π^2 log( ζ/(32 π^2) ) = 0,

 which is not satisfied for any f. As such, low-frequency emanation from   X   as   = α^'' k_^ or  ^2 = α_^ k_^ 2 + α_^ k_^ 2 is not supported for nonmagnetic crystals.


 In Fig. <ref> we present first band(s)    for   a selection of crystals  comprising  nonmagnetic cylinders (μ_ =1) in   air, and describe their evolution as ε_ is varied from 0 < ε_≤ -2.1.  For values   0 > ε_≳ -0.5274   (<ref>), a  single first band   emanates from the Γ point,   analogously to the Γ-emerging band in Fig. <ref>. At  ε_≈ -0.5274,   a band emerges from the X point, giving rise to two first band surfaces at  low frequencies, as in Fig. <ref>. This X-emergent surface eventually  forms  a double degeneracy at the Γ point at ε_≈ -0.56, and thereafter, as shown in Fig. <ref>,   becomes the first band surface,  pushing the existing Γ centered band to higher frequencies. As we proceed in ε_, the    origin of the first band travels   along the Γ M and Γ X symmetry planes to the M point at ε_ = -1 (the anomalous resonance condition), as shown in Fig. <ref>. Thereafter, a new first band emerges at the M point whose origins move along   high symmetry planes towards the Γ point. At the same time, this emergent band pushes the existing M centered band to higher frequencies. Both of these behaviors are demonstrated in Fig. <ref>.   The new first band eventually takes the form given in Fig. <ref> where it emanates from both the Γ and X points simultaneously at ε_ = (ζ + 16   π^2)/(ζ - 16   π^2) ≈ -1.896.  As we proceed in ε_, the first eigenfrequency at X becomes nonzero and the first band    emanates from the Γ point alone, as demonstrated by Fig. <ref>. Undoubtedly further behaviors are observed beyond ε_<-2.1,  however, this falls  outside the scope of the present work.  Note that in Figs. <ref>  and <ref>   the dipolar approximations (<ref>)   and (<ref>) are superposed, respectively,  and  show    excellent  agreement at low frequencies.  The   examples above, in consideration with (<ref>), emphasize that  the absence of a band from a single Bloch coordinate does not preclude the   emergence of   bands from multiple Bloch coordinates simultaneously; this observation has important implications for determining the existence of   band gaps at low frequencies. We emphasize that the band structure smoothly transitions between the examples shown above. 


     In summary, we have     determined new low-frequency E-k relations for photons in   2D  photonic crystals. These   relations, and the conditions for their existence, are given explicitly for  first bands with  origins at the Γ, X, and M points of a square lattice. In general, we have found that    sign changes in the properties of the constituents are required for the first band to originate from coordinates away from Γ at low frequencies. We have also demonstrated that photonic crystals can possess low-frequency E-k relations   with origins   at one or more arbitrary Bloch coordinates.     Given that all conventional  photonic crystals   possess either a low-frequency band gap or a band surface emanating from   Γ, this work has significant implications for the homogenization   of periodic media,  theoretical descriptions of light propagation in complex media, as well as  future  photonic  crystal designs.    In the latter case, the closed-form conditions we obtain represent a powerful design tool for determining filling fractions 
and background materials for a given   cylinder material, and vice versa.  However, an important consequence of the M-point condition coinciding with   the anomalous resonance condition is that   the slope of the     M-emerging band (<ref>) is not necessarily accurate beyond a dipolar truncation;  further investigations are    required to accurately determine   behaviors      near resonance (see Supplemental Material <cit.>).     Preliminary results for   hexagonal lattices at the anomalous resonance reveal  a band surface emerging from  the K point, implying  that     the first band   originates from the furthest edge of the irreducible Brillouin zone when on resonance.  Away from the anomalous resonance condition, we believe that experimental validation is feasible for the crystals we describe, as all   emergence conditions   (i.e., (<ref>) and (<ref>))      are valid for complex-valued ε_ and μ_. When ε_ and μ_ possess     moderate loss,  and ε_ and μ_ satisfy the necessary emergence conditions, we find    the   band diagrams to be unchanged  (see Supplemental Material <cit.>).     A natural next step for this work is    1D and 3D photonic structures; 2D photonic crystals were only chosen  for analytical convenience. Finally, we emphasize that our approach  extends  readily to phononic and other   systems. 
 
M.J.A.S. acknowledges discussions with R.C. McPhedran, S. Guenneau, R.V. Craster, support from the ERC (279673), and support from the EPSRC (EP/L018039/1). 
 
    

§ SUPPLEMENTAL MATERIAL FOR “VIOLATING THE ENERGY-MOMENTUM PROPORTIONALITY OF PHOTONIC CRYSTALS IN THE LOW-FREQUENCY LIMIT”


 

§   EXTENDED PROBLEM FORMULATION
 
 We begin by considering the time-harmonic form  of the source-free Maxwell equations, in a   non-dispersive and lossless system. The domain comprises a two-dimensional array of infinitely extending isotropic cylinders that are periodically positioned at the coordinates of a square lattice in the (x,y) plane and embedded in an isotropic background material.    That is, we  examine the band structure of  the Maxwell  equations
 
    ∇× = -,      ∇× = ,
    ∇· = 0,       ∇· = 0,

 having assumed that fields possess the dependence  exp(- t). We then assume the linear constitutive relations

    = ε_0ε_,      = μ_0 μ_,

where the relative permittivity and permeability tensors take the scalar values
 
    ε_ = {[  ε_ ∈Ω_;  ε_ ∈Ω_ ]. ,
       μ_ = {[  μ_ ∈Ω_;  μ_ ∈Ω_ ].,

 where Ω_ denotes the region inside the cylinder (of radius r = a^') positioned in the centre of the fundamental unit cell Ω_WSC, and Ω_ represents the remainder of the unit cell (i.e., Ω_ = Ω_WSC\Ω_). Substituting the constitutive relations (<ref>)  into the Maxwell system (<ref>) admits the vector wave equation

    ∇×ε_^-1∇× -   ^2 c_0^-2μ_ = 0.

We then  impose that all electromagnetic fields have the dependence   exp( k_z z) and that  wave propagation in the direction of the fibres is forbidden, i.e.    k_z=0.  Subsequently,   the  Bloch vector takes the form = (k_x, k_y,0) where we introduce _ = (k_x, k_y) as the in-plane Bloch wave vector.   This restriction on the wave vector  permits a decomposition of the    Maxwell wave equation (<ref>) above  into two decoupled   scalar problems for the   two   orthogonal polarisation states; the first with the electric field  polarised in the direction z (along the cylinders), and the second with the magnetic field  polarised in the direction z.  

  In this work, we examine the latter polarisation problem and subsequently    the Maxwell wave equation (<ref>) reduces to the study of   the two-dimensional Helmholtz equation
 
    ∇_⊥·( ε_^-1∇_⊥ H_) +  ^2 c_0^-2μ_ H_ = 0,

 where ∇_⊥≡(∂_,∂_) denotes the in-plane gradient,  denotes the angular   frequency, and c_0 is the speed of light in vacuum. Similarly,    the   electromagnetic continuity conditions at the edges of the cylinders   take the form
 
 
    ( H_^ - H_^)|_r=a^' =0,   ( ε_^-1∂_r H_^ - ε_^-1∂_r H_^) |_r =a^' = 0,

 corresponding to continuity of H_z and E_θ at the   boundary, respectively.  Additionally, we   impose the Floquet–Bloch conditions  
 
    H_( + _p)  = H_() ^_·_p ,

  
 where _p= (ma,na) for m,n ∈ℤ is the real-lattice vector  for a square lattice of period a, where p is multi-index notation for the pair (m,n). The system for the   E_ polarised field is given by     exchanging ε_ and μ_ in the above, and subsequently all band diagrams shown are thus also attainable under the E_ polarisation. The   general solution to the homogeneous two-dimensional Helmholtz equation    (<ref>) is well-known and takes the form 
  
    H_ = {[ H_^  ,   ∈Ω_,; H_^  ,   ∈Ω_, ].

 where  
 
 
 
    H_^   = ∑_m=-∞^∞[ 𝒜_m J_m(   n_ r) +  ℬ_m Y_m( n_ r) ]^ m θ,  
    
     H_^   = ∑_m=-∞^∞𝒞_m J_m( n_ r)  ^ m θ,


 with (r, θ) denoting polar coordinates in the (x,y) plane, = /c_0, n_ = √(ε_)√(μ_), and n_ = √(ε_)√(μ_).  Note that in this work, we choose the positive branch of the square root function as the principal root and that √(ε)√(μ)≠√(εμ) in general; such a result only holds when both parameters are positive-valued. From our definition of n above, we observe that both ε and μ must be negative-valued in order to obtain a negative-index material (n<0).  The functions J_m(z) and Y_m(z) denote Bessel functions of the first and second kind. Substituting these general solutions  (<ref>) into the boundary conditions (<ref>) we obtain


    𝒞_m = 𝒜_m J_m( n_  a^') +  ℬ_m Y_m( n_ a^') /J_m( n_ a^') ,

and

    𝒜_m = -_m ℬ_m,
 
where
  
    _m =     Z_^-1 Y_m(n_ a^') ∂_r J _m(r) |_r=n_ a^'   - Z_^-1 J_m(n_ a^') ∂_r Y_m(r) |_r=n_ a^'/   Z_^-1 J_m(n_ a^') ∂_r J_m (r) |_r=n_ a^'    -Z_^-1 J_m(n_ a^') ∂_r J_m (r) |_r=n_ a^',
 

having introduced  the impedances  Z_ = √(μ_) / √(ε_) and Z_ = √(μ_) / √(ε_).
  

In the presence of an incident plane wave, we determine the scattered fields  from all   other cylinders  in the array  (i.e., the multiple scattering contribution) using    Green's second identity,     admitting the  Rayleigh identity   <cit.>
  
    𝒜_l =  ∑_m=-∞^∞ (-1)^l+m S_m-l^ℬ_m,

where S_m^ denote lattice sums  (defined in Section <ref>) and we have omitted the incident field as we seek the Bloch modes of the structure.   Substituting (<ref>) into (<ref>) above we obtain the   system
 
    _l  ℬ_l + ∑_m=-∞^∞ (-1)^l+m S_m-l^ℬ_m = 0,

 whose vanishing determinant   gives the dispersion equation for the homogeneous periodic problem. To obtain a numerically tractable form of (<ref>),   it is necessary to   truncate parameters, i.e. to evaluate (<ref>) for l,m = -L, …, L.   For conventional photonic crystals, the convergence of the band surfaces   improves with larger truncation values L. However, at low frequencies    a monopole (L=0) or dipole truncation (L=1) is     regarded as  a good approximation for the fully converged  dispersion equation.
  

 

§ ASYMPTOTIC ANALYSIS OF DISPERSION EQUATION
 
 In the preceding section we determined the form of the dispersion equation (<ref>) for a two-dimensional photonic crystal of infinitely extending cylinders, which to monopolar order (L=0) takes the form
  
    _0 + S_0^  =0,

 and to dipolar order (L=1) takes the form
 
    ( _0 + S_0^) ( [ _1 + S_0^]^2 - | S_2^|^2 ) + ( S_1^)^2  [ S_2^]^∗  
    - 2 | S_1^|^2 ( _1 + S_0^) 
    +  ( [ S_1^]^∗)^2    S_2^ =0,

where we have used the identities  _-m = _m  and S_-l^ = [ S_l^]^∗,  where ∗ denotes   complex conjugation  <cit.>.  In this section, we  obtain  conditions necessary for the existence of    band surfaces, through a close examination   of (<ref>) and (<ref>)   in the low-frequency   regime.  
 	


 §.§ Asymptotic behaviour of  dispersion relations near the Γ point

We begin by  evaluating asymptotic expressions for the scattering coefficients N_m defined in  (<ref>) in the limit as = /c_0 → 0.  These terms exhibit the leading order behaviour 

    lim_→ 0_m = 
    {[   O(^-2)     m =0; O( ^-2m)     m >0;          ]. ,
 
and subsequently we  write  

 
    lim_→ 0{_0 }≈∑_m=-2^∞ D_m ^m,     lim_→ 0{_1 }≈∑_m=-2^∞ E_m ^m,


where   the first few terms of these series  are given by


    D_-2   =    4/π   a^' 2ε_1/μ_ - μ_    , 
    
    D_0   =     (4 μ_  -3   μ_)   ε_μ_ -μ_^2  ε_/2 πε_ (μ_-μ_)^2 + 2    /π[  γ_ +log ()  ] 
    + 1/πlog(a^' 2ε_μ_/4) ,     
    
    E_-2   =  4/π  a^' 2ε_μ_ε_ + ε_/ε_ - ε_  , 
    
    E_0   =   μ_ (ε_ + ε_)(5 ε_ - 3 ε_)- 4  μ_ε_^2/2πμ_ (ε_ - ε_)^2   
       + 2/π[  γ_  + log() ]  + 1/πlog(  a^' 2ε_μ_/4) ,


and γ_≈ 0.577215 … denotes the Euler–Mascheroni constant. 
 
 Next we consider the   leading-order behaviours of the S_l^ terms present in (<ref>) and (<ref>). Semi-analytical representations for these sums were previously     determined by  <cit.>, using  numerical estimates for certain terms. Using   the approach in <cit.> (details in Section <ref>) we obtain the fully closed forms 



    lim_k_→ 0lim_→ 0 S_0^Y    ≈ - 4/a^21/k_^2 - (n_)^2  - 2/πlog(n_ )    
    -   2 γ_/π  - 1/πlog(  a^2 π/Γ(14)^4 )  ,  
    lim_k_→ 0lim_→ 0 S_1^Y    ≈  - 4 /a^2k_/n_1/k_^2 -     (n_)^2^θ_ + /π k_/n_^θ_   ,  
    lim_k_→ 0lim_→ 0 S_2^Y    ≈4/a^2k_^2/(n_)^21/k_^2 - (n_)^2^2 θ_  
    +  k_^2/(n_)^2[   ( Γ(14)^8 /384 π^5) ^-2 θ_  - 1/2π^2 θ_],


where (k_,θ_) is the polar representation of the Bloch vector _, and Γ(z) is the Gamma function.  We now assume      that a band surface   emanates from the Γ point at low frequencies in the form = α k_, where α is real and positive-valued,   and substitute this    into  (<ref>) to obtain  
 
    lim_k_→ 0lim_→ 0 S_l^Y(,_;n_) ∼ O(^-2),
 
for l = 0,1,2, and where we treat log() terms as O(^0). Subsequently, we expand the   lattice sums in the form
  
    lim_k_→ 0lim_→ 0{ S_0^Y }≈∑_m=-2^∞ A_m ^m ,     lim_k_→ 0lim_→ 0{ S_1^Y }≈∑_m=-2^∞ B_m ^m,     lim_k_→ 0lim_→ 0{ S_2^Y }≈∑_m=-2^∞ C_m ^m,

 where summation is made   over non-zero orders, all of which are even due to symmetry. Substituting  (<ref>) and (<ref>)  into the monopolar dispersion equation (<ref>) and collecting terms  in advancing orders of  reveals a hierarchy of equations describing the leading-order behaviour of the dispersion relation at low frequencies and the conditions necessary for the existence of the   band. 


 The first non-vanishing order  for the monopolar dispersion equation (<ref>) is given by

    O(^-2):     A_-2 + D_-2 = 0,

where (<ref>) and (<ref>) admit  the dispersion relation  

    = {1/ε_[  μ_ + f   (μ_ - μ_)  ]}^1/2 k_,

 Subsequently,   provided  ε_μ_ + f ε_ (μ_ - μ_)  >0, a low-frequency band surface   emanates from the Γ point within a monopolar approximation.   
  

  Proceeding to the dipolar truncation (<ref>), the first non-vanishing order is 
 
    O(^-6):     (A_-2)^3  + (A_-2)^2 (D_-2 + 2   E_-2 ) + 
     A_-2[ (E_-2)^2 - |C_-2|^2 + 2   D_-2 E_-2 -2 |B_-2|^2      ] 
    
    + (C_-2)^∗ (B_-2)^2 + [ (B_-2)^∗]^2 C_-2  
      - 2 |B_-2|^2 E_-2  + D_-2[ (E_-2)^2 -|C_-2|^2 ]= 0 ,

which gives the low-frequency dispersion relation 
 

    = { ( ε_ + ε_)+f    ( ε_ -   ε_ ) /[ε_μ_+f    ε_ (μ_-μ_)]  [  (ε_+ε_)-f   (ε_-ε_)]  }^1/2 k_,

 where f = π a^' 2/a^2 is the filling fraction, provided

    ( ε_ + ε_)+f    ( ε_ -   ε_ ) /[ε_μ_+f    ε_ (μ_-μ_)]  [  (ε_+ε_)-f   (ε_-ε_)]  >0.

 Note that (<ref>) is recovered under the substitution   ε_ = ε_ in (<ref>)  above,  and that an effective refractive index is readily obtained via n_ = 1/α. In this setting,  the expressions for the effective permittivity and permeability are  decoupled, motivating the representations
 
    ε_ = ε_( (ε_+ε_)-f   (ε_-ε_)/ ( ε_ + ε_)+f    ( ε_ -   ε_ ) )   = ε_ + 2ε_ (ε_ - ε_)f/  (ε_ + ε_) -    (ε_ - ε_)f,     μ_ =    μ_+f      (μ_-μ_).

 Having determined the low-frequency dispersion relation   in the vicinity of the Γ point, we now consider the asymptotics of the dispersion relation both at low frequencies and in the vicinity of other high symmetry points.




 §.§ Asymptotic behaviour of  dispersion relations near the M point
 
The asymptotic forms of the first few  S_l^ sums   in the limit of vanishing frequency and  in the vicinity of the M point     are given by
 


    lim_k_→ Mlim_→ 0{ S_0^Y }   ≈    -2 γ_/π - 2/πlog(  n_) 
    - 1/πlog( 4 a^2 π/Γ(14)^4),  
    lim_k_→ Mlim_→ 0{ S_1^Y }   ≈/π k_^'/ n_^θ_^'   , 
    lim_k_→ Mlim_→ 0{ S_2^Y }   ≈ 
      -  k_^' 2/ ( n_)^2 [   1/π^2( Γ(14)^8/128 π^3)  ^-2 θ_^'
      +  1/2π^2 θ_^'],


 where (k_^' ,θ_^' ) is the polar representation of _^' = _ -, with = (π/a,π/a).    A derivation of these expression  is presented in  Section <ref>.   As in  Section <ref> we now construct an ansatz assuming that a low frequency band surface emanates from the M point in the form  = α^' k_^', where α^' is real and positive-valued. From the lattice sum asymptotics in (<ref>) this dependence yields

    lim_k_→ Mlim_→ 0 S_l^Y(,_;n_) ∼   O(^0) ,

 and so   we introduce the expansions
 
    lim_k_→ Mlim_→ 0{ S_0^Y }≈∑_m=0^∞ A_m^'^m ,   lim_k_→ Mlim_→ 0{ S_1^Y }≈∑_m=0^∞ B_m^'^m,    lim_k_→ Mlim_→ 0{ S_2^Y }≈∑_m=0^∞ C_m^'^m,

where summation is made over even orders, as the odd orders remain identically zero.  Next we substitute the expansion for _0 in (<ref>) and S_0^ in  (<ref>)   into the monopolar dispersion equation (<ref>) and collect terms  in  to obtain a hierarchy of conditions. The first two of these equations are  



    O(^-2):       D_-2 = 0, 
    
    O(^0):       D_0 + A_0^'  = 0,


where from (<ref>), it follows that the O(^-2) equation is never satisfied for finite permittivity and permeability values, and subsequently a band surface is not supported from the M point at low frequencies within a monopolar approximation. This finding is consistent with the fact that a monopole field is   symmetric and subsequently  cannot satisfy quasi-periodic boundary conditions at the unit cell edges.


In contrast, from  the dipolar dispersion equation (<ref>) we   obtain  



    O(^-6):      (   E_-2)^2  D_-2 = 0, 
    
    O(^-4):        E_-2[  (A_0^'  + D_0) E_-2   + 2  (A_0^'  + E_0 )  D_-2] = 0,   
      
    O(^-2):        [t]  D_-2[ (A_0^'  + E_0) ^2  -|C_0^' |^2    ]  
    +  (E_-2)^2 (A_2^'  + D_2)  
    
      + 2 E_-2[   (A_0^')^2     -  |B_0^' |^2  + A_2^'  D_-2 +D_0   E_0 + A_0^'  ( D_0  +       E_0) + D_-2  E_2 ] =0,


 where from inspection,   the O(^-6) and O(^-4) equations are     satisfied when E_-2 = 0, or  ε_ = -ε_ following (<ref>).
 Furthermore, the O(^-2) equation is simplified by  the  E_-2 = 0 requirement, reducing it to the form 

    O(^-2):     D_-2(  A_0^'  + E_0   - |C_0^' | )  (  A_0^' + E_0   + |C_0^' | ) = 0,

where   A_0^'  + E_0   + |C_0^'| = 0 ultimately gives the dispersion relation 
 
 
    = {1/  8 π^ 2| 64 π^4 ^4 θ_^' + Γ(14 )^8 |^1/2[ ε_μ_ + 2 ε_μ_log( 16   π^2/f   Γ(14)^4) ]^-1/2} k_^'.

 
  Thus within a dipolar approximation, a band surface is supported from the M point at   low frequencies (with a slope described by (<ref>)) provided  both
  
  
    ε_ = -ε_,

   and   

    μ_ 
            >[ log( f^2  Γ(14)^8/2^8   π^4) ]  μ_ ,
 
       
  are satisfied, where the prefactor in (<ref>) is negative-valued for 0<f≲0.914 (and for reference, f = π/4 ≈ 0.785 represents   dense packing   for a square array of cylinders).  The expression  (<ref>) follows from the condition that  α must be real and positive-valued.  
      
 
 

 §.§ Asymptotic behaviour of  dispersion relations near the X point
 
The asymptotic  forms of the first few    S_l^at both low frequencies and  in   the vicinity of the X point   are given by



    lim_k_→ Xlim_→ 0 S_0^Y     ≈   -   2 γ_/π  - 2/πlog(n_ ) 
    - 1/πlog(  8 a^2 π/Γ(14)^4 ) ,  
    lim_k_→ Xlim_→ 0 S_1^Y     ≈/π k_^/n_(  ^θ_^  +  Γ  (14 )^4 /8 π^2^-θ_^) , 
    lim_k_→ Xlim_→ 0 S_2^Y     ≈1/ (n_)^2( Γ(14)^4/4 a^2 π^2) 
      + Γ(14)^4 /8 π^3(1- k_^ 2/  (n_)^2 )  
      +   k_^ 2/ (n_)^2[ 1/π^2( Γ(14)^8 /128  π^3)^-2 θ_^ 
      - 1  /2 π^2 θ_^]  ,


 where (k_^,θ_^) is the polar representation of _^ = _ -. The derivation for the   lattice sum  expressions above is given in  Section <ref>.  In an identical manner to Sections <ref> and <ref>,  we  now assume that a linear band surface emanates from the X point at low frequencies. That is, we  substitute the form = α^'' k_^, where α^'' is real and positive-valued,   into (<ref>) and observe that 
 
    lim_k_→ Xlim_→ 0 S_l^Y(,_;n_)    ∼ O(^-floor(l^2/2)),

for l=0,1,2.   Such limit behaviour     admits  the   expansions

  
    lim_k_→ Xlim_→ 0{ S_0^Y }≈∑_m=0^∞ A_m^''^m ,   lim_k_→ Xlim_→ 0{ S_1^Y }≈∑_m=0^∞ B_m^''^m,    lim_k_→ Xlim_→ 0{ S_2^Y }≈∑_m=-2^∞ C_m^''^m,


where summation is made over even orders, as the odd orders remain identically zero. For the monopolar dispersion equation,  we   obtain an analogous hierarchy of equations to (<ref>) which does not permit a low frequency band at X. However, for the dipolar dispersion relation we obtain  


    O(^-6):       D_-2( E_-2 - |C_-2^''| )  ( E_-2+ |C_-2^''| ) = 0, 
    
    O(^-4):          |C_-2^''|^2 (A_0^'' + D_0) + 2 |C_-2^''| D_-2 |C_0^'' |   - E_-2 (D_0 E_-2 + A_0^'' (2 D_-2 + E_-2)  
    + 2 D_-2 E_0) = 0,


where the O(^-6) system is satisfied provided  |C_-2^''|=  -E_-2 or equivalently

    ε_ = (  f   Γ(14)^4 - 16   π^2 / f   Γ(14)^4+ 16  π^2) ε_.

The above condition   simplifies  the O(^-4) expression (<ref>) to the form

    2 D_-2 E_-2(   A_0^'' + E_0 +  |C_0^''| ) =0.

   Substituting (<ref>), the coefficients for (<ref>), and (<ref>) into the reduced condition A_0^'' + E_0 +  |C_0^''|  =0      yields  a lengthy expression of the form

    β_1  + β_2   | β_3 + β_4  α^ 2/α^ 2| = 0,

where β_j are constants.   To ensure that   β_1 = 0 it is necessary that

    μ_ = (    f   Γ(14)^4 ( f    Γ(14)^4  - 64     π^2 )  
     + 512   π^4 log( f   Γ(14)^4/(32   π^2) )   /(f   Γ(14)^4  - 16 π^2)^2)μ_,

which then simplifies  (<ref>) to reveal  the low-frequency dispersion relation


    = (  ( 16 Γ(14)^4 π^2+ 64 π^4 ^2θ_^ -Γ(14)^8 ^-2θ_^) /16 Γ(14)^4 π^2 ε_μ_)^1/2 k_^.

However, we remark that the α^''  in (<ref>) is only real-valued  along the paths Γ X and XM (i.e., for θ_^ = π and θ_^ =π/2, respectively). Numerical investigations confirm the existence of a band surface across the entire Brillouin zone,  where the contours are elliptical at low frequencies. After amending our dispersion relation ansatz to   ^2 = α_ (k_-π/a)^2 + α_ k_^2 and using (<ref>) above we finally obtain  
 

    ^2 = (  16 Γ(14)^4 π^2 + 64 π^4 - Γ(14)^8 /16 Γ(14)^4  π^2 ε_μ_) k_^ 2 + (  16 Γ(14)^4 π^2 - 64 π^4 + Γ(14)^8 /16 Γ(14)^4  π^2 ε_μ_)  k_^ 2,

 for the first band surface as → 0 and  for _ in the vicinity of X. Subsequently, the requirement that α_,α_>0 is equivalent to ε_μ_ >0.

 



§ LOW-FREQUENCY DESCRIPTIONS FOR NON-MAGNETIC MEDIA
 


In this section, we provide a condensed outline of the conditions for a band to emanate from the Γ, X, and M points of the crystal  when it comprises non-magnetic media, within a dipolar approximation.   As mentioned in Section <ref> the leading-order   behaviour of the scattering coefficients _m differ  for photonic crystals made from non-magnetic and magnetic materials. In the former instance, the scattering coefficients take the form
 
    lim_→ 0_m = 
    {[  O(^-4)    m =0; O(^-2m)    m >0;         ]. ,
 
in contrast to (<ref>). Expanding the first two  _n coefficients    as
 
    lim_→ 0{_0 }≈∑_m=-4^∞D_m ^m,    lim_→ 0{_1 }≈∑_m=-2^∞E_m ^m,

we obtain the first few terms 



    D_-4   =    32/π   a^' 4ε_1/ε_ - ε_    , 
    D_-2   = 16 /3 π a^' 2ε_ε_ - 2 ε_/ε_ - ε_,     
    E_-2   =  4/π  a^' 2ε_ε_ + ε_/ε_ - ε_  , 
    E_0   =  1/2π[ ε_ + 3 ε_/ε_ - ε_ + 4  γ_      +2   log( a^' 2ε_/4) + 4  log() ].


Substituting the _m expansions (<ref>)   and  the  Γ centred   lattice sum expansions (<ref>) into the dipolar dispersion relation (<ref>) gives  an O(^-8) condition, and not an O(^-6) condition as before, revealing

    = {1/ε_(ε_ + ε_) + f (ε_ - ε_)/(ε_ + ε_) - f (ε_ - ε_)}^1/2 k_,

which is identical to  (<ref>) above with the replacement μ_, μ_↦ 1.


 Similarly, for the M point, substituting (<ref>) and (<ref>) into the dipolar dispersion equation we obtain a system of conditions analogous to (<ref>). For non-magnetic crystals, the O(^-8) and O(^-6) conditions are   satisfied when 
 

    ε_ = -ε_,

 with the O(^-4) condition giving

    = {1/  8 π^ 2| 64 π^4 ^4 θ_^' + Γ(14 )^8 |^1/2[ ε_   + 2 ε_log( 16   π^2/f   Γ(14)^4) ]^-1/2} k_^',

 
which is identical to (<ref>) but with μ_,μ_↦ 1. The analogue to α>0 in (<ref>) now takes the simpler form  ε_  >0.  For the X point,  the leading-order O(^-8) condition (c.f.  O(^-6) before) is satisfied for  
 

    ε_ = (  f     Γ(14)^4 - 16   π^2 / f       Γ(14)^4 + 16  π^2) ε_,

(which is identical to (<ref>)). However, the O(^-6) condition  (c.f.  O(^-4) before) takes the form  (<ref>) where β_1 = 0 requires that 

    f  Γ(14)^4   + 8π^2 + 16π^2log( 32 π^2/f   Γ(14)^4) = 0,

 
which does not hold for any f. Hence, non-magnetic photonic crystals do not support first bands emerging from the X point as  = α^'' k_^ or ^2 = α_ (k_-π/a)^2 + α_ k_^2. To clarify,  (<ref>) demonstrates that the ansatz = α^'' k_^ in a non-magnetic crystal cannot satisfy β_1 = 0 in one of the orders, and subsequently, a dispersion relation cannot be extracted. Hence, that particular ansatz does not describe a low-frequency behaviour of non-magnetic crystals. From our numerical study, we have observed     = C_1   |_| + C_2  |_ - | in non-magnetic crystals,  which   is an entirely different ansatz, with different low-frequency behaviour, that will undoubtedly be supported. However, we believe that obtaining    descriptions for the latter case   falls outside the scope of the present work. Note that for non-magnetic crystals, if we have negative permittivities in the constituent materials, then the refractive index and impedance are purely imaginary.
  



§ NUMERICAL RESULTS
 
In this section, we provide further details on the numerical calculation of the  band structures     of a two-dimensional  photonic crystal. For this task, we     search for the vanishing determinant of the system <cit.>
      
    ( δ_lm +  sgn(_l) / (|_l| |_m|)^1/2 (-1)^l+m S_m-l^) 𝒟_m = 0,

            where δ_lm denotes the Kronecker delta function and  𝒟_m = (  |_m|)^1/2ℬ_m,  in place of   the original system (<ref>).  The above representation scales the value of the determinant in absolute value terms, as the scattering coefficients and lattice sum terms   become increasingly large for l,m →∞. The lattice sums S_l^ are evaluated using  the expressions    given in  <cit.>, which  explicitly relate  one-    and two-dimensional array sums for any Bravais lattice configuration. For the one-dimensional array sums,  we use   <cit.>, which are   accelerated forms of the   expressions first derived in <cit.>, and then evaluate the necessary correction factors in <cit.> to obtain two-dimensional sums for a square lattice. The approach of evaluating grating sums and then correction factors  <cit.> is preferred to evaluating     accelerated lattice sum expressions   <cit.> directly, as we have found   that the former   approach   is much more numerically stable    for → 0.

 

  

 
 For reference, in Fig. <ref> we demonstrate the shape of the modes in the vicinity of their emanation points. In Fig. <ref> we present the mode for ε_ = 1 and μ_ = 2 in air with _≈(0,0) which shows that the mode is approximately constant as → 0. This behaviour is typical  for a crystal with such a low material contrast <cit.>. In Fig. <ref>, we have ε_ = -1 and μ_ = 2 in air with _≈(π/a,π/a) which shows that the mode takes the form of a line dipole largely concentrated to the boundary of the cylinder. The dipole response is  also oriented along θ = 3π/4 to ensure that    the   Bloch  conditions are met  along all edges of the unit cell (i.e., anti-periodicity). Finally, in Fig. <ref> we give the mode for ε_≈ -0.53 and μ_≈ -10.22 in air with _≈(π/a,0) which also shows a    dipole response concentrated to the boundary of the cylinder. This dipole response      is oriented along θ = π/2 which ensures that  the Bloch conditions are satisfied.
 

In Fig. <ref> we present  band diagrams for a crystal satisfying the conditions for Γ, M, and X point emanation (given by (<ref>),   (<ref>) and (<ref>), and (<ref>) and (<ref>), respectively), but with a higher truncation  value L=5. In all instances, the dashed red line represents the dipolar estimate of the first band surface. In Fig. <ref>  we observe that there   are negligible differences between the few first band surfaces of a conventional photonic crystal, upon comparing  Fig. <ref> (L=5) with Fig. 1a of the letter (L=1).  In Fig. <ref>,   we impose the M-point conditions for L=1 and evaluate the dispersion equation for L=5, and see   that  the   conditions correctly predict the material parameters for which M emanation occurs, but that the asymptotics are unable to correctly predict the slope.  The higher bands are also different, but this is not unexpected as  higher order multipole terms are generally required to reproduce both the dispersion relation and the modal fields at higher frequencies.  We emphasize that this   sensitivity is observed close to the  M point conditions alone (and is not observed      away from ε_ =-ε_), due to its correspondence with the anomalous resonance condition  <cit.>.  We discuss this in  further detail below. In Fig. <ref>  we impose  the X-point conditions corresponding to L=1 for   L=5, where we observe that the first two bands are well-approximated, but that the origin of the first band moves slightly away from X. Note that the  first band is   re-centred about X after slightly  perturbing the values of ε_ and  μ_. These variations in the band diagrams   emphasise  that  the conditions and descriptions we obtain for Γ  and X point emanation are good approximations for the full systems, but that our description of the first band emerging from M are specific to the dipolar    dispersion equation and are not necessarily  accurate for higher truncations.      Analytically  determining high-truncation conditions   poses a significant  challenge, as  the   dispersion equation becomes highly intractable   for large L.  


 
    

Returning to the discussion of anomalous resonances, these may be understood by analogy to (<ref>); the zeros of N_l and N_m correspond to poles in the system at the  angular orders l and m, which at low frequencies, accumulate about ε_ =-ε_  with increasing truncation L (hence why anomalous resonances are also known as accumulation points). As such, in the vicinity of  ε_ =-ε_, the band diagrams for L=1 do  not necessarily reflect the band diagrams for the full system, and we have found that results using other numerical tools, such as plane-wave expansion methods <cit.> and finite-element methods, exhibit extremely strong numerical instability with variations in the number of plane waves and the maximum  element size, respectively.   Preliminary results suggest that as we increase the truncation parameter L →∞, the slope of the first band  tends to α^'→ 0. We also find that   the first band is only present at truncation values L=1,2,5,6,9,10,… (OEIS A042963) otherwise a low-frequency gap is observed; it is possible that the truncation order is tied to the existence of the first band   (i.e., the truncation must be chosen to ensure that the system possesses the correct   symmetries).  This was observed earlier, where M point emanation was not supported within a monopolar truncation.   We remark that the issues observed with the M point are not overcome  through the use of other solution procedures, such as  finite-element method (FEM) solvers (i.e.,   Comsol).

    In fact, conventional FEM solvers are unable to validate the bands in Fig. 1 (Letter) centered about M or X at low frequencies. For the M point, this is because the   band diagrams returned by FEM solvers exhibit   strong instability near the anomalous resonance condition. For the X point,  most commercially available solvers   cannot evaluate bands when   both the permittivity and permeability are negative-valued. However,  it is possible to examine all band diagrams in Fig. 3 (Letter) using FEM solvers,  excluding Figs. 3c and 3d, which cannot be validated due to their closeness to the anomalous resonance condition. We find that   the  conditions we derive herein give an excellent approximation to  those obtained using a full FEM solver approach  (i.e., the condition for Fig. 3e (Letter) is ε_≈ -1.9 using FEM solvers and not ε_≈ -1.896). Overall, we find  that the dipole approximation gives a good   description of both the conditions and the first band(s) when compared to those obtained with FEM solvers. For reference, the anomalous resonance condition ε_ = -ε_ is identical for   square and hexagonal lattices of cylinders   <cit.>, and   that   no analogy was found for the hexagonal lattice under the X point conditions for a square lattice.

Regarding the issue of loss, if we incorporate loss in the background material (i.e., consider complex ε_ and μ_), and then use the  vanishing-determinant conditions to calculate  the optical properties of the cylinder   (thus obtaining complex  ε_ and μ_), then the corresponding band diagram (under the complex ω and real k representation) is robust to moderately large imaginary values. Satisfying only the real parts of the background or cylinder materials constants   may cause a  gap to emerge about a band origin   at low frequencies.  Determining materials which satisfy the vanishing determinant conditions for complex-valued constants  may pose something of a challenge, however it should be possible to engineer values for the permittivity and permeability with  metamaterials as the cylinder media, for example. 


  
 

  As a final comment, we remark that   large condition numbers are intrinsic   for our cylindrical-Mie theory method, i.e. the discretised version of (<ref>) or (<ref>), even for conventional photonic crystals. Such  numerical ill-conditioning poses a challenge  for accurately determining   band structures   (as a numerically singular matrix  may return   spurious band surfaces). For L=1, we find that   condition numbers are not prohibitively  large, but due to the large values taken by _l and S_l^ as L→∞, that ill-conditioning may pose a significant challenge for larger truncations (and is not corrected by the use of (<ref>) in place of  (<ref>)).

  


 

§  ASYMPTOTIC FORMS OF  LATTICE SUMS FOR A SQUARE ARRAY 

 The   cylindrical-Mie  solution for a two-dimensional lattice of cylinders embedded in a uniform background material with permittivity ε_ and permeability μ_ features      sums  of the form <cit.>

    S_l^(,_;n_) ='∑_(m,n)∈ℤ^2  Y_l ( n_ R_p) ^ l ϕ_p^_·_p,

where Y_l(z) is a Bessel function  of the second kind, = /c_0 is the scaled angular frequency, n_ = √(ε_)√(μ_), _ is the in-plane Bloch vector, (R_p, ϕ_p) represents polar coordinates for the real lattice generator _p = (ma,na), for m,n∈ℤ, and prime notation denotes summation over the entire lattice excluding  (m,n)=(0,0), where the summand is singular. Here, the subscript p represents a double index ranging over all m and n.  The lattice sum above is conditionally convergent, however, numerous procedures exist for  accelerating convergence <cit.>. One of these methods is to consider the reciprocal-space representation of (<ref>); this is readily obtained  by comparing forms of  the quasi-periodic Green's function for the Helmholtz operator defined by

    {Δ + (n_)^2 } G( , ^') = -2π∑_pδ( - ^' - _p) ^_·_p.

The comparison procedure is outlined in <cit.> and gives the absolutely convergent representation 
 
    S_l^(,_;n_ ) J_l( n_ξ)  = -Y_0(n_ξ) δ_l0  
    - 4^l/A∑_h J_l(Q_h ξ) ^ l Θ_h/Q_h^2 - (n_)^2,

where A denotes the area of the fundamental unit cell, (K_h,ψ_h) are polar coordinates for the reciprocal lattice generator _h =  2π (m/a,n/a) for m,n ∈ℤ, (Q_h,Θ_h) denote polar coordinates for _h = _h + _, and (ξ,γ) represent  polar coordinates for ξ =  - ^' (where we remark that ψ_h and γ are unused in the above, but are defined for completeness). Here, the subscript h is a double index ranging over all m and n. Note that ξ, whilst defined as the difference between the source and field coordinates for the purposes of evaluating the Green's function G(,^'), now represents an arbitrary vector of finite length in the unit cell, as the lattice sums S_l^ in (<ref>) are independent of spatial coordinates.

 

 §.§ Asymptotic representations of dynamic lattice sums S_l^ near Γ 
  We now present a brief summary of the method for determining the asymptotic forms of S_l^ in the low frequency and vanishing Bloch vector limit, following the approach outlined  in Appendix A of  <cit.>. This  begins  by isolating the h = (0,0) term in   (<ref>), and applying Graf's addition theorem <cit.> to the surviving summand, admitting 
 
    S_l^(,_;n_)   = - Y_0(n_ξ)/J_0(n_ξ)δ_l0 - 4^l/A J_l(k_ξ)/J_l(n_ξ)^ l θ_/ k_^2 - (n_)^2
    
     - 4^l/A 1/J_l(n_ξ)'∑_h∑_m=-∞^∞1/Q_h^2 - (n_)^2{ (-1)^m J_l+m (K_h ξ)   J_m(k_ξ) ^ (l+m) ψ_h -  m θ_},
 
 where (k_,θ_) denotes the polar representation of the Bloch vector _. A Taylor series expansion for small  and k_ gives
 
    lim_k_→ 0 lim_→ 01/Q_h^2 - (n_)^2≈1/K_h^2(  1 -   2 k_/K_hcos(ψ_h - θ_)  + 2 k_^2/K_h^2cos(2{ψ_h - θ_} ) + k_^2 +(n_)^2/K_h^2 + …),

 which after an appropriate truncation of m in (<ref>), and considerable algebraic manipulation, admits the asymptotic representations  
 


    lim_k_→ 0lim_→ 0{ S_0^Y }   ≈ - Y_0(n_ξ)/J_0(n_ξ) - 4/A{( 1/k_^2 - (n_)^2 + S_0,0,2)   J_0(k_ξ)/J_0(n_ξ) + 2 k_  S_1,0,3J_1(k_ξ)/J_0(n_ξ)  + …},
    lim_k_→ 0lim_→ 0{ S_1^Y }   ≈  - 4^θ_/A{( 1/k_^2 - (n_)^2 + S_0,0,2)    J_1(k_ξ)/J_1(n_ξ) - k_  S_1,0,3J_0(k_ξ)/J_1(n_ξ)+ …}, 
    lim_k_→ 0lim_→ 0{ S_2^Y }   ≈4 ^2 θ_/A{( 1/k_^2 - (n_)^2 + S_0,0,2)    J_2(k_ξ)/J_2(n_ξ) - k_  S_1,0,3J_1(k_ξ)/J_2(n_ξ)   + k_^' 2   S_2,0,4J_0(k_ξ)/J_2(n_ξ)}
            +4 ^-2 θ_/A{  S_4,4,2 J_2(k_ξ)/J_2(n_ξ)  + k_  S_3,4,3 J_1(k_ξ)/J_2(n_ξ)    +k_^  2  S_2,4,4 J_0(k_ξ)/J_2(n_ξ)   + …},


 where we define the double Schömilch series <cit.>
 
    S_l,m,n  (ξ;τ,a)= ∑_h^'J_l(K_h ξ)/K_h^n^ m ψ_h.

Using the explicit representations for S_l,m,n in <cit.>, or the recurrence relation procedure in  <cit.>, we obtain    
 

 
    S_0,0,2   = - a^2/2πlog(ξ) + a^2/4πlog(4 π  a^2 /Γ(14)^4) +ξ^2/4,         
      S_2,4,4   =  a^2 ξ^2 Γ(14)^8/3 π^5 2^12    -   ξ^4 Γ(14)^8/9 π^4 2^11  +    ξ^6 Γ(14)^8/15 π^3 a^2 2^13, 
    
      S_1,0,3   = - a^2 ξ/4πlog(ξ) + a^2 ξ/8πlog( 4 π a^2  /Γ(14)^4) + a^2 ξ/8π + ξ^3/16,         
    S_3,4,3   =  ξ^3Γ(14)^8/9 π^4 2^10   -   ξ^5 Γ(14)^8/15 π^3 a^2  2^11,
    
    S_2,0,4   = - a^2 ξ^2/16πlog(ξ) + a^2 ξ^2/8πlog(4 π a^2  /Γ(14)^4)  + 3a^2 ξ^2/64π + ξ^4/96,         
     S_4,4,2   =   ξ^4  Γ(14)^8/15 π^3 a^2 2^10 ,

  
  for a square lattice of period a. The  S_l,m,n expressions above in tandem with the small argument expansions      <cit.>  



    lim_ξ→ 0lim_k_→ 0lim_→ 0Y_0(n_ξ)/J_0(n_ξ)   ≈2/π( γ_ + log( n_ξ/2) )+ (n_)^2 ξ^2/2 π + …, 
    lim_ξ→ 0lim_k_→ 0lim_→ 0J_0(k_ξ)/J_0(n_ξ)   ≈ 1 + (n_)^2 - k_^2/4ξ^2 + …, 
    lim_ξ→ 0lim_k_→ 0lim_→ 0J_1(k_ξ)/J_1(n_ξ)   ≈k_/n_ + k_/n_(n_)^2 - k_^2/8ξ^2 + …, 
    lim_ξ→ 0lim_k_→ 0lim_→ 0J_0(k_ξ)/J_1(n_ξ)   ≈2/n_ξ + (n_)^2 - 2k_^2/4 n_ξ + ( 1/32k_/n_ - 1/16n_/k_ + 1/48(n_)^3/k_^3) k_^3 ξ^3 + …, 
    lim_ξ→ 0lim_k_→ 0lim_→ 0J_1(k_ξ)/J_0(n_ξ)   ≈k_ξ/2 + 2(n_)^2 - k_^2/16 k_ξ^3 + …, 
    lim_ξ→ 0lim_k_→ 0lim_→ 0J_0(k_ξ)/J_2(n_ξ)   ≈8/(n_)^2 ξ^2 + 2( 1/3 - k_^2/(n_)^2) + ( -1/6 + 1/8k_^2/ (n_)^2 + 5/144(n_)^2/k_^2) k_^2 ξ^2 + …, 
    lim_ξ→ 0lim_k_→ 0lim_→ 0J_1(k_ξ)/J_2(n_ξ)   ≈4 k_/(n_)^2ξ + 2(n_)^2 - 3k_^2/6(n_)^2 k_ξ + ( -1/24 + 1/48k_^2/(n_)^2 + 5/288(n_)^2/k_^2) k_^3 ξ^3 + …, 
    lim_ξ→ 0lim_k_→ 0lim_→ 0J_2(k_ξ)/J_2(n_ξ)   ≈k_^2/(n_)^2 +k_^2/(n_)^2(n_)^2 - k_^2/12ξ^2 +…,


  finally admit the representations (<ref>) for S_0^, S_1^, and S_2^ both at low frequencies and in the vicinity of the Γ point, after extensive algebraic manipulation. 
 

 §.§ Asymptotic representations of dynamic lattice sums S_l^ near M 
  To obtain closed-form representations for  S_l^ in the vicinity of  different symmetry points (and also for vanishing ) we extend  the procedure   above. The approach is identical up to (<ref>), however, to consider behaviour near the M point it is necessary to modify   (<ref>) appropriately.  This is achieved by decomposing  the translated reciprocal lattice generator as _h = _h + _ = (_h + )+ (_ - )  = _h^' + _^', where = (π/a,π/a). Substituting this decomposition into (<ref>), we recover an identical    expression   to before, but with the replacements  _h ↦_h^' and _↦_^'    and an updated   limit argument. Subsequently we obtain analogous expressions to (<ref>) but without the  singular terms 1/(k_^2 - (n_)^2), as   light lines are not present near   M   at low frequencies. That is, for the M point we obtain



    lim_k_→ Mlim_→ 0{ S_0^Y }   ≈ - Y_0(n_ξ)/J_0(n_ξ) - 4/A{   S_0,0,2^ J_0(k_^'ξ)/J_0(n_ξ) + 2 k_^'  S_1,0,3^J_1(k_^'ξ)/J_0(n_ξ)},
    lim_k_→ Mlim_→ 0{ S_1^Y }   ≈  - 4^θ_^'/A{   S_0,0,2^ J_1(k_^'ξ)/J_1(n_ξ) - k_^'  S_1,0,3^J_0(k_^'ξ)/J_1(n_ξ)}, 
    lim_k_→ Mlim_→ 0{ S_2^Y }   ≈4 ^2 θ_^'/A{  S_0,0,2^J_2(k_^'ξ)/J_2(n_ξ) - k_^'  S_1,0,3^J_1(k_^'ξ)/J_2(n_ξ)   + k_^' 2   S_2,0,4^J_0(k_^'ξ)/J_2(n_ξ)}
            +4 ^-2 θ_^'/A{  S_4,4,2^ J_2(k_^'ξ)/J_2(n_ξ)  + k_^'  S_3,4,3^ J_1(k_^'ξ)/J_2(n_ξ)    +k_^' 2  S_2,4,4^ J_0(k_^'ξ)/J_2(n_ξ)},


where we emphasise that  ξ is an arbitrary vector in the first Brillouin zone that remains untranslated, (k_^' ,θ_^' ) is the polar representation of _^' = _ -, and  

    S_l,m,n^ (ξ;τ,a)= ∑_h^'J_l(K_h^'ξ)/K_h^' n^ m ψ_h^',

where (K_h^' ,ψ_h^' ) is the polar representation of _h^' = _h +.  As highlighted in  <cit.>, analytical expressions for  S_l,m,n^, corresponding to      Bloch vectors centred around the M point, are easily obtained by evaluating   a selection of  rectangular and square Γ-centred sums and using  the multi-set identity <cit.>
 

    Ω^M(,a) = Ω(,2a) - Ω(/2,2a) + Ω(,a) - Ω(2,a),

where Ω are  reciprocal lattice sets defined by

    Ω(τ,a) = {2π/A (h_1 b _1^' + h_2 a _2^' )  | h_1,h_2 ∈ℤ},

 
with A = ab representing the area of the unit cell (for a rectangular lattice of periods a and b in x and y respectively), _j^' are   basis vectors for the reciprocal lattice,   and τ = (b/a).  The   multi-set identity (<ref>) gives the following expressions



    S_0,0,2^   = - a^2/2 πlog( ξ)  + a^2/4 πlog(  16 π a^2 /Γ  (14 )^4) ,          S_2,4,4^  =   - ξ^6 Γ(14 )^8 /5 π^3 a^2  2^14 + ξ^4 Γ(14 )^8 /9 π^4 2^10   -a^2 ξ^2 Γ(14 )^8  /π^5 2^12 ,  
    
      S_1,0,3^   =  a^2 ξ/8 π -  a^2 ξ/4 πlog( ξ) +  a^2 ξ/8 πlog(  16 π a^2 /Γ  (14 )^4) ,        S_3,4,3^  =    ξ^5 Γ(14 )^8 /5  π^3 a^2 2^12  -ξ^3 Γ(14)^8 /9 π^4 2^9 ,  
     
    S_2,0,4^   =   3 a^2 ξ^2 /64 π  - a^2 ξ^2 /16 πlog(ξ) + a^2 ξ^2 /32 πlog(  16 π a^2 /Γ  (14 )^4)  ,          S_4,4,2^   =  -ξ^4 Γ(14 )^8 /5  π^3 a^2 2^11.


Substituting the S_l,m,n^ expressions   (<ref>) and the Bessel function expansions (<ref>) into the expansions for S_l^   (<ref>) above  ultimately gives the final expressions (<ref>) used in Section <ref>.
   
  


 §.§ Asymptotic representations of  S_l^ at the X point
 As with the M point approach, in order to consider behaviour near the X point it is necessary to modify   (<ref>) appropriately. Decomposing  _h = _h + _ = (_h + )+ (_ - ), where = (π/a,0) into (<ref>)   we obtain
 


    lim_k_→ Xlim_→ 0 
    S_0^≈ -Y_0(n_ξ)/J_0(n_ξ) 
     - 4/AJ_0(k_^ξ)/J_0(n_ξ)(  S_0,0,2^ 
    + (n_)^2 S_0,0,4^  
     +   k_^ 2 S_0,0,4^)
    
     - 4/AJ_0(k_^ξ)/J_0(n_ξ) ^-2θ_^( k_^ 2  S_0,2,4^) 
     - 4/AJ_0(k_^ξ)/J_0(n_ξ)^2θ_^(  k_^ 2  S_0,2,4^),

 
    lim_k_→ Xlim_→ 0 S_1^≈
    -4^θ_^/AJ_1(k_^ξ)/J_1(n_ξ)( S_0,0,2^  
    +(n_)^2 S_0,0,4^  
    + k_^ 2  S_0,0,4^) 
    +4^θ_^/AJ_0(k_^ξ)/J_1(n_ξ)( 
      k_^   S_1,0,3^) 
    
    +4^-θ_^/AJ_0(k_^ξ) /J_1(n_ξ) ( 
     k_^   S_1,2,3^)  
    -4^-θ_^/AJ_1(k_^ξ)/J_1(n_ξ)( k_^ 2   S_0,2,4^) ,

 
    lim_k_→ Xlim_→ 0 S_2^≈4  /AJ_0(k_^ξ)/J_2(n_ξ)( S_2,2,2^ + (n_)^ 2 S_2,2,4^  +  k_^ 2 S_2,2,4^) 
     -4  /AJ_1(k_^ξ)/J_2(n_ξ)( k_^ S_1,2,3^) 
    
    + 4  /AJ_2(k_^ξ)/J_2(n_ξ)( k_^ 2 S_0,2,4^)  
    + 4 ^2θ_^/AJ_2(k_^ξ)/J_2(n_ξ)(  S_0,0,2^ + (n_)^2 S_0,0,4^ +    k_^ 2 S_0,0,4^) 
    
         +  4 ^2θ_^/AJ_1(k_^ξ) /J_2(n_ξ) (  -  k_^ S_1,0,3^) 
         +   4 ^2θ_^/AJ_0(k_^ξ) /J_2(n_ξ) (k_^ 2 S_2,0,4^) 
    +4  ^-2θ_^/AJ_0(k_^ξ)/J_2(n_ξ)( 
     k_^ 2S_2,4,4^),


where _^ = _ - with polar representation (k_^,θ_^), _h^ = _h + with polar representation (K_h^ ,ψ_h^ ), and

    S_l,m,n^ (ξ;τ,a)= ∑_h^'J_l(K_h^ξ)/K_h^ n^ m ψ_h^,

 Note that  in contrast to the representations for S_l^ in (<ref>) and (<ref>),  the closed-form representations for S_l^ near the X point have a considerably different    structure. This is because  only two-fold symmetry is possessed near the X point at low frequencies.  The  S_l,m,n^ are evaluated using  the multi-set identity   <cit.>

    Ω^X(,a) =  Ω(/2,2a) - Ω(,a)  ,

to obtain



    S_0,0,2^   =  -a^2/2 πlog (ξ )  +a^2 /4 πlog(    32 π a^2   /Γ  (14 )^4  ) ,   S_0,0,4^  =   3 a^4   /2 π^4 ζ(2) β(2) -a^2 ξ^2/8π + a^2 ξ^2/8πlog(ξ) 
    - a^2 ξ^2/16 πlog(32 π a^2   /Γ(14)^4), 
    
    S_0,2,4^   = -a^2 Γ  (14 )^4/128 π^3 ξ^2 + Γ  (14 )^4/1024 π^2ξ^4 - a^4/16 π^4 G_4^(2)(2),     S_1,0,3^  =   a^2 ξ/8π - a^2 ξ/4πlog(ξ) + a^2 ξ/8 πlog( 32 π a^2 /Γ  (14 )^4     ), 
    
    S_1,2,3^   = a^2 Γ  (14 )^4  /64 π^3ξ - Γ  (14 )^4  /256 π^2ξ^3,    S_2,0,4^  =   3 a^2 ξ^2/64π - a^2 ξ^2/16πlog(ξ)  + a^2 ξ^2/32 πlog( 32  π a^2 /Γ(14)^4), 
    
    S_2,2,2^   =  Γ  (14 )^4/128 π^2ξ^2,    
    S_2,2,4^   =  a^2 Γ  (14 )^4/256 π^3ξ^2  - Γ  (14 )^4/3 π^2 2^9ξ^4,   
    S_2,4,4^  =ξ^6 Γ(14)^8 /5 π^3 a^2 2^16  - ξ^4 Γ(14)^8/9 π^4 2^11 + a^2 ξ^2 Γ(14)^8/π^5 2^12,


where ζ(s) denotes the Riemann zeta function and β(s) the Dirichlet beta function. 
 Note that closed form expressions are not generally available for phase-conjugated Eisenstein series G_n^(m) when m<n (for reference,  G_4^(2)(2) ≈ -30.854212880849047 and is given as a rapidly convergent sum in <cit.>). 
Substituting  (<ref>) and (<ref>) into the asymptotic approximations (<ref>) above gives the final expressions (<ref>) presented in Section <ref>.

As a final comment, we remark that obtaining closed-form expressions   for S_l^ about   other points in the Brillouin zone may pose something of a challenge as these do not have multiset representations using origin centered lattices. The X and M points are special cases where multiset representations are possible due to symmetry, despite not being origin centered lattices themselves.

 




  
   




