## Abstract

The article demonstrates uncommon manifestation of spatial dispersion in low refractive index contrast 3D periodic dielectric composites with periods of about one tenth of the wavelength. First principles simulations by the well established plane wave method reveal that spatial dispersion leads to appearance of additional optical axes and can compensate anisotropy in certain directions.

© 2017 Optical Society of America

## 1. Introduction

Spatial dispersion effects in elastic light interaction with 3D periodic structures result in various phenomena unusual to isotropic media, and strongly influence propagation of surface waves [1]. Currently manifestations of the spatial dispersion are well-known both for natural crystals [1,2] and artificial periodic structures being either photonic crystals or metamaterials [3].

Nonlocality in natural crystals is weak at optical wavelengths, and results in directional variations of mode propagation constant at scales of about 10^{−5} −10^{−6} [4]. However, cubic crystals appear to possess seven optical axes, when such deviations are considered. This effect called spatial-dispersion-induced birefringence, was already known by Lorentz [5], and nowadays is taken into account by manufacturers of projection systems for deep ultraviolet lithography set-ups [6,7].

On the other hand, nonlocality in artificial photonic structures – photonic crystals and metamaterials – is much more prominent [3,8] and draws a lot of attention due to promising applications. Strong spatial dispersion in photonic crystals results in superprism effect, possibility for perfect imaging, and generally provides advanced potential for spatial filtering and beam manipulation [9,10]. Interpretation of phenomena related to spatial dispersion appears to be handy with representation by isofrequency contours and surfaces [11–15]. For certain classes of metamaterials isofrequency contours were demonstrated to be dramatically distorted in comparison with those of natural crystals highly depending on structure parameters [16,17]. Spatial-dispersion-induced birefringence being small for natural materials appeared to have prominent values for metamaterials with resonant metallic inclusions [18, 19]. Authors of [20] demonstrated a theoretical possibility to stop the light in metamaterial waveguide coming from an interplay between the spatial dispersion and anisotropy.

An intermediate position between weak spatial dispersion in natural crystals and prominent effects in photonic crystals and metamaterials can be filled with low refractive index periodic dielectric structures of period-to-wavelength ratio of the order of 0.1, so that no diffraction occurs. Paper [21] provided an evidence that such composites can exhibit quite unusual optical properties by demonstrating 2D periodically patterned waveguide to manifest trirefringence. Current technology level of the direct laser writing technique allows fabricating even 3D periodic high-quality dielectric composites with periods as small as tenths of a micron for operating at mid-infrared and telecom wavelengths [22]. To estimate an order of magnitude of spatial dispersion effects in such structures consider permittivity decomposition under assumption that structure is non-gyrotropic:

**q**is wavevector. The second term scales as (Λ/

*λ*)

^{2}with Λ and

*λ*being characteristic period of a composite and the wavelength respectively [1]. For structures of interest set Λ/

*λ*~10

^{−2}−10

^{−1}, which yields |

*ε*

^{−1}− (

*ε*

^{−1})

^{(0)}| ~ 10

^{−4}− 10

^{−2}. Taking into account that optical anisotropy in such structures in the limit of infinitely small period amounts to 10

^{−3}−10

^{−2}for low index dielectrics and to 10

^{−1}for high index dielectrics [23,24], it becomes clear, that spatial dispersion terms can compete with anisotropy of the first term ${\epsilon}_{\alpha \beta}^{(0)}$, which may lead to interesting effects for low refractive index contrast composites operating away from resonances.

## 2. Plane wave expansion method

In order to investigate spatial dispersion in the described parameter range we utilized a variation of the plane wave expansion method for first principle simulation of photonic structures [25–27]. By using such first principle approach, from the one hand, we avoid restricting ourselves by decomposition of Eq. (1) and consequent loss of accuracy due limited applicability of Eq. (1) for finite period structures. From the other hand, this method provides means for estimation of a range of applicability of truncated series in Eq. (1). The plane wave expansion method can be formulated on the basis of both differential and integral equations. In the current research we use the electric field formulation, which is based on the volume integral equation solution of the Maxwell’s equations [28]

*ε*is some constant, and

_{b}*µ*

_{0}is the permeability of vacuum),

*α*,

*β*= 1, 2, 3 enumerate Cartesian axes,

*δ*is the Kronecker symbol, and the multiplier before source

_{αβ}*J*(

_{β}**r**′) is the free space dyadic Green function. ${E}_{\alpha}^{inc}(\mathbf{r})$ denotes components of some pre-defined incident field, and

*E*(

_{α}**r**) is unknown field. Taking the source as to be generated by spatial permittivity changes

**J**= −

*iω*(

*ε*−

*ε*)

_{b}**E**one obtains a self-consistent equation on the unknown electric field. In the discreet 3D Fourier space Eq. (2) becomes

**m**= (

*m*

_{1},

*m*

_{2},

*m*

_{3}) enumerates Fourier harmonics, possible

**k**component values are ${k}_{\alpha \mathbf{m}}={k}_{\alpha}^{inc}+{m}_{\alpha}{K}_{\alpha},{K}_{\alpha}=2\pi /{\mathrm{\Lambda}}_{\alpha}$, and Λ

*are periods of the composite structure along coordinate axes. Notation for wavevectors here is different from Eq. (1) to distinguish wavevectors q of propagating modes. Propagation constants and polarizations of eigen modes in a given 3D periodic composite can be retrieved from solution of either an eigenvalue problem coming from Eq. (3) with zero incident field ${E}_{\alpha \mathbf{m}}^{inc}$, or a resonant analysis of Eq. (3). We implemented the latter possibility through the use of an efficient complex pole search algorithm [23].*

_{α}Special precautions were taken to control simulation accuracy at each step of the algorithm, as we are interested in effects which can occur for relatively small effective refractive index contrasts starting from values of about 10^{−4}. To satisfy this requirement we first set a sufficiently low convergence criterion for the Generalized Minimal Residual method for solution of the linear equation system (3), and, second, apply the second order Richardson extrapolation to improve calculated propagation constants and harmonic field amplitudes for subsequent numbers of Fourier harmonics used in simulations. This allowed obtaining mode propagation constants with 5-7 significant digit precision.

## 3. Spatial dispersion effects

In what follows we provide three illustrative examples of spatial dispersion manifestation. First example concerns the above mentioned existence of seven optical axes in cubic crystals, and aims at quantitatively estimating an impact of the second and higher order terms in Eq. (1). Consider a scaffold structure shown in Fig. 1(a). For all examples let the refractive index of constituting rods be *n _{c}* = 1.6, rod size-to-period ratio be 1/3, and surrounding medium be the air. Figure 1(b) shows isofrequency surfaces for two propagating modes and directions of optical axes in crystal with cubic lattice and Λ/

*λ*= 0.1. In order to clearly represent small differences we choose the scale along each axis

*X*to be [

_{α}*n*(

**ŝ**) −

*n*

_{0}]ŝ

*, with*

_{α}*n*

_{0}being the same constant for all three Cartesian axes, and

*n*(

**ŝ**) being calculated propagation constant in direction

**ŝ**=

**q**/|

**q**|.

Dependence of absolute values of maximum effective refractive index difference for two propagating modes in the structure, which occurs in 〈110〉 directions, from the period-to-wavelength ratio is shown in Fig. 2(a) in doubly logarithmic scale. The curves reveal that series truncation by second order terms in Eq. (1) is valid roughly up to periods about Λ/*λ* ≈ 0.2 independently of composite dielectric constant within low contrast range *n _{c}* −1 ~0.5. Then, for periods lower than 0.2 we can extract numerical values of coefficients before second order terms in permittivity decomposition. For cubic crystals there are only three independent components of tensor

*ζ*, namely,

_{αβγδ}*ζ*

_{1111},

*ζ*

_{1212}, and

*ζ*

_{1122}[29]. This permits reducing Eq. (1) to [7,29]:

*p*

_{1},

*p*

_{2},

*p*

_{3}depending on Λ/

*λ*. Our simulations do not allow to account for longitudinal term

*n*

^{2}

*p*

_{3}

*ŝ*since its multiplication by the electric field yields values of order of 10

_{α}ŝ_{β}^{−5}at maximum (because

*p*

_{3}~ (Λ/

*λ*)

^{2}~ 10

^{−2}−10

^{−3}, and |π/2 − ∠(

**ŝ**,

**E**)| < 10

^{−3}). So we put

*p*

_{3}= 0. The other two coefficients

*p*

_{1,2}are obtained from the Helmholtz equation

*n*

^{2}(

**ŝ**)

**ŝ**×

**ŝ**×

**E**+

*ε*

**E**= 0 by substitution of calculated mode propagation constants

*n*(

**ŝ**). This yields the expected quadratic dependence of

*p*

_{1,2}, as Fig. 2(b) shows, with

*p*

_{1}= −0.014(Λ/

*λ*)

^{2}, and

*p*

_{2}= 0.153(Λ/

*λ*)

^{2}. Value

*ε*

^{(0)}= 1.31866 was precalculated as a limit of effective refractive index along coordinate axis at Λ/

*λ*= 0 under assumption of its quadratic behaviour. Besides, one can calculate decomposition of inverse permittivity

*ε*

^{−1}similarly. This example also demonstrates that maximum propagation index difference in cubic crystals practically does not exceed 10

^{−2}when the spatial dispersion is described by the second order term for structures of interest.

For the second example modify the considered scaffold structure so as one of periods to be different from the two others: Λ_{1} = Λ_{2}, Λ_{3} = 1.1Λ_{1} and consider tetragonal lattice. In this case we get effectively uniaxial medium in the limit Λ_{1} →0 with optical axis coinciding with the third Cartesian axis. Figure 3(a) demonstrates the well-known refractive index sphere and ellipsoid in the above described coordinates. Increase of periods results in appearance of four additional optical axes in planes (110) analogously to the first case of cubic lattice, and another four optical axes in planes (100) and (010), which come from splitting of cubic lattice axes in directions 〈100〉 and 〈010〉 due to symmetry reduction. These additional axes rise from 〈001〉 direction at infinitely small periods. Then, as Λ_{1}/*λ* increases, axes of the first set lying in planes (110) tend to approach 〈111〉 directions, while axes of the second set approach 〈100〉 and 〈010〉. This means that spatial dispersion terms tend to compensate anisotropy in 〈100〉 and 〈010〉 directions.

In case of tetragonal lattice tensor *ζ _{αβγδ}* has seven independent components [29], and one of them defines the longitudinal input. Therefore, six constants

*p*

_{1}, …,

*p*

_{6}define the permittivity tensor without the longitudinal term:

_{1}/

*λ*)

^{2}, and shows, that for periods up to Λ/

*λ*~ 0.15 they write

*p*

_{1}= −0.023(Λ/

*λ*)

^{2},

*p*

_{2}= −0.017(Λ/

*λ*)

^{2},

*p*

_{3}= −0.024(Λ/

*λ*)

^{2},

*p*

_{4}= 0.163(Λ/

*λ*)

^{2},

*p*

_{5}= 0.17(Λ/

*λ*)

^{2}, and

*p*

_{6}= 0.16(Λ/

*λ*)

^{2}. Diagonal components are ${\epsilon}_{11}^{(0)}=1.3179$, and ${\epsilon}_{33}^{(0)}=1.32019$.

In the third example we put all three periods to be different Λ_{2} = 0.9Λ_{1}, Λ_{3} = 1.1Λ_{1} to consider an orthorhombic lattice, so the crystal behaves as effectively biaxial medium in the limit Λ_{1} → 0, which is shown in Fig. 5(a). With increase of periods there remain two optical axes up to Λ_{1}/*λ* ≈ 0.15. At larger values of Λ_{1}/*λ* the lattice isofrequency surface transits to a state with 10 optical axes. In comparison with the first example the rotational symmetry is broken along all coordinate axes, and consequently all three optical axes of the cubic crystal are split in two at sufficiently large periods (Fig. 5(b)). Longitudinal term being extracted, permittivity tensor tensor writes in this case via 11 parameters:

_{1}/

*λ*ratio is shown in Fig. 4(b). Corresponding numerical coefficients for small period values Λ

_{1}/

*λ*< 0.15 write

*p*

_{1}= −0.003(Λ/

*λ*)

^{2},

*p*

_{2}= 0.152(Λ/

*λ*)

^{2},

*p*

_{3}= 0.155(Λ/

*λ*)

^{2},

*p*

_{4}= 0.037(Λ/

*λ*)

^{2},

*p*

_{5}= 0.157(Λ/

*λ*)

^{2},

*p*

_{6}=−0.003(Λ/

*λ*)

^{2},

*p*

_{7}= 0.158(Λ/

*λ*)

^{2},

*p*

_{8}= 0.019(Λ/

*λ*)

^{2},

*p*

_{9}= 0.151(Λ/

*λ*)

^{2},

*p*

_{10}= 0.15(Λ/

*λ*)

^{2}, and

*p*

_{11}= −0.003(Λ/

*λ*)

^{2}with diagonal components ${\epsilon}_{11}^{(0)}=1.31887$, ${\epsilon}_{22}^{(0)}=1.316136$, and ${\epsilon}_{33}^{(0)}=1.321055$.

## 4. Conclusion

To conclude, our work reveals existence of nontrivial spatial dispersion effects in low refractive index contrast 3D periodic dielectric composites operating in the effective medium regime. For sufficiently large period values artificial crystals possessing uniaxial and biaxial properties in the limit of infinitely small periods are demonstrated to have 9 and 10 optical axes respectively. In addition, quantitative values of effective parameters of composites retrieved from results of first-principle simulations yield a tendency for the spatial dispersion to compensate anisotropy in certain directions. These effects can be regarded as strong in a sense that spatial dispersion terms have the same order of magnitude as anisotropy in the considered structures.

We verified that small absorption (Im(*n*) ≲ 10^{−2}) in a dielectric material constituting a composite does not qualitatively affect the results, and only leads to appearance of proportionally small imaginary parts of effective refractive indices. The described effects may be used in applications which require advanced spectral and spatial filtering of the electromagnetic radiation. One may expect, that a plane wave interacting with the investigated structures should be tolerant to small defects inevitably present in experimental samples: due to small period-to-wavelength ratio the wave phase appreciably changes at scales of dozens of periods, thus, making only average period filling and length to play a role.

Presented retrieval of *ζ _{αβγδ}* values is a kind of a numerical homogenization procedure for periodic dielectric composites with determined range of validity. This procedure is a direct use of effective propagation constants and corresponding polarizations, obtained through the first-principles calculations, in the Maxwell’s equations under assumption of validity of decomposition in Eq. (1). Probably, analogous results can be obtained from various procedures elaborated within the homogenization theory, which is still being intensively developed (e.g., see review [30] and references therein). However, we did not intend to evaluate all possible approaches. Main purpose of the paper is to demonstrate the effects appearing in periodic structures and to quantitatively characterize them.

## Funding

Russian Ministry of Education and Science (16-19-2014/K); Russian Foundation for Basic Research (16-29-11747-ofi-m and 16-07-00837).

## References and links

**1. **V. M. Agranovich and V. L. Ginzburg, *Spatial Dispersion in Crystal Optics and the Theory of Excitons* (Wiley-Interscience, 1966).

**2. **J. H. Burnett, Z. H. Levine, and E. L. Shirley, “Intrinsic birefringence in calcium fluoride and barium fluoride,” Phys. Rev. B **64**, 241102 (2001). [CrossRef]

**3. **R. C. Rumpf, “Engineering the dispersion and anisotropy of periodic electromagnetic structures,” Sol. State. Phys. **66**, 213–260 (2015).

**4. **V. L. Ginzburg, “Electromagnetic waves in isotropic and crystalline media characterized by dielectric permittivity with spatial dispersion,” Sov. Phys. JETP **34**(7), 1096–1103 (1958).

**5. **H. A. Lorentz, *Collected Papers II* (Martinus Nijohff, The Hague, 1936), pp. 1–119.

**6. **A. Serebriakov, E. Maksimov, F. Bociort, and J. Braat, “The effect of intrinsic birefringence in deep UV-lithography,” Proc. SPIE **5249**, 624 (2004). [CrossRef]

**7. **J. H. Burnett, Z. H. Levine, E. L. Shirley, and J. H. Bruning, “Symmetry of spatial-dispersion-induced birefringence and its implications for CaF2 ultraviolet optics,” J. Micro/Nanolith. MEMS MOEMS **1**(3), 213–224 (2002). [CrossRef]

**8. **K. Vynck, D. Felbacq, E. Centeno, A. I. Cčbuz, D. Cassagne, and B. Guizal, “All-dielectric rod-type matamaterials at optical frequencies,” Phys. Rev. Lett. **102**, 133901 (2009). [CrossRef]

**9. **M. Notomi, “Manipulating light with strongly modulated photonic crystals,” Rep. Progr. Phys. **73**, 096501 (2010). [CrossRef]

**10. **L. Maigyte, V. Purlis, J. Trull, M. Peckus, C. Cojocaru, D. Gailevičius, M. Malinauskas, and K. Staliunas, “Flat lensing in the visible frequency range by woodpile photonic crystals,” Opt. Lett. **38**, 2376–2378 (2013). [CrossRef] [PubMed]

**11. **H. Kosaka, A. Tomita, T. Kawashima, T. Sato, and S. Kawakami, “Splitting of triply degenerate refractive indices by photonic crystals,” Phys. Rev. B **62**(3), 1477–1480 (2000). [CrossRef]

**12. **M. Notomi, “Negative refraction in photonic crystals,” Opt. Quantum. Electron. **34**, 133–143 (2002). [CrossRef]

**13. **T. Baba and M. Nakamura, “Photonic crystal light deflection devices using the superprism effect,” IEEE J. Quantum. Electron. **38**(7), 909–914 (2002). [CrossRef]

**14. **Y. Loiko, C. Serrat, R. Herrero, and K. Staliunas, “Quantitative analysis of subdiffractive light propagation in photonic crystals,” Opt. Commun. **269**, 128–136 (2007). [CrossRef]

**15. **L. Maigyte and K. Staliunas, “Spatial filtering with photonic crystals,” Appl. Phys. Rev. **2**, 011102 (2015). [CrossRef]

**16. **M. A. Gorlach and P. A. Belov, “Nonlocality in uniaxially polarizable media,” Phys. Rev. B **92**, 085107 (2015). [CrossRef]

**17. **C. Fietz, Y. Urzhumov, and G. Shvets, “Complex k band diagrams of 3D metamaterial/photonic crystals,” Opt. Expr. **19**(20), 19027–19041 (2011). [CrossRef]

**18. **A. V. Chebykin, M. A. Gorlach, and P. A. Belov, “Spatial-dispersion-induced birefringence in metamaterials with cubic symmetry,” Phys. Rev. B **92**, 045127 (2015). [CrossRef]

**19. **M. A. Gorlach, S. B. Glybovsky, A. A. Hurshkainen, and P.A. Belov, “Giant spatial-dispersion-induced birefringence in metamaterials,” Phys. Rev. B **93**, 201115 (2016). [CrossRef]

**20. **K. L. Koshelev and A. A. Bogdanov, “Interplay between anisotropy and spatial dispersion in metamaterial waveguide,” Phys. Rev. B **94**, 115439 (2016). [CrossRef]

**21. **M. C. Netti, A. Harris, J. J. Baumberg, D. M. Whittaker, M. B. Charlton, M. E. Zoorob, and G. J. Parker, “Optical trirefringence in photonic crystal waveguides,” Phys. Rev. Lett. **86**(8), 1526–1529 (2001). [CrossRef] [PubMed]

**22. **S. Maruo and J. T. Fourkas, “Recent progress in multiphoton microfabrication,” Laser Photon. Rev. **2**(1–2), 100–111 (2008). [CrossRef]

**23. **A. A. Shcherbakov and A. V. Tishchenko, “3D periodic dielectric composite homogenization based on the Generalized Source Method,” J. Opt. **17**, 065101 (2015). [CrossRef]

**24. **M. Che, Z-Y. Li, and R.-J. Liu, “Tunable optical anisotropy in three-dimensional photonic crystals,” Phys. Rev. A **76**, 023809 (2007). [CrossRef]

**25. **K. M. Leung and Y. F. Liu, “Full vector wave calculation of photonic band structures in face-centered-cubic dielectric media,” Phys. Rev. Lett. **65**(21), 2646–2649 (1990). [CrossRef] [PubMed]

**26. **K. M. Ho, C. T. Chan, and C. M. Soukoulis, “Existence of a photonic gap in periodic dielectric structures,” Phys. Rev. Lett. **65**(25), 3152–3155 (1990). [CrossRef] [PubMed]

**27. **S. G. Johnson and J. D. Joannopoulus, “Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis,” Opt. Expr. **8**(3), 173–190 (2001). [CrossRef]

**28. **P. M. Morse and H. Feshbach, *Methods of theoretical physics* (McGraw-Hill, 1953).

**29. **J. F. Nye, *Physical Properties of Crystals* (Clarendon, 1985).

**30. **C. R. Simovski, “On electromagnetic characterization and homogenization of nanostructured metamaterials,” J. Opt. **13**, 013001 (2011). [CrossRef]