Concurrency of anisotropy and spatial dispersion in low refractive index dielectric composites

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.


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][12][13][14][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: where 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 ε (0) αβ , which may lead to interesting effects for low refractive index contrast composites operating away from resonances.

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][26][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] where k 2 b = ω 2 ε b µ 0 (ε b is some constant, and µ 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 inc α (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 Here vector index m = (m 1 , m 2 , m 3 ) enumerates Fourier harmonics, possible k component values are k αm = k inc α + m α K α , K α = 2π/Λ α , 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 inc αm , 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.

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-towavelength 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]: with parameters 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   Fig. 1(a) has tetragonal lattice with period relations Λ 2 = Λ 1 , Λ 3 = 1.1Λ 1 . Axis scales are the same as in Fig. 1(b).
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.

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 firstprinciples 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.