Nonlocal homogenization of PT-symmetric multilayered structures

Unique and highly tunable optical properties of PT-symmetric systems and metamaterials enable a plenty of entirely new linear and nonlinear optical phenomena with numerous applications, e.g., for designing sub-diffraction lenses, non-reciprocal devices etc. Therefore, the artificial media with the PT symmetry attract ever-increasing attention and are now a subject for intensive investigations. One of the commonly used methods providing information about optical response of artificial nanostructural media is a so called effective medium theory. Here, we examine the possibility of utilizing the effective medium theory for a comprehensive analysis of PT-symmetric multilayered systems composed of alternating loss and gain slabs. We show that applicability of local effective material parameters (or Maxwell Garnett approximation) is very limited and cannot be exploited for a prediction of exceptional points marking a PT symmetry breaking. On the contrary, nonlocal bianisotropic effective medium parameters can be reliably used, if the thickness of a unit cell is much smaller than the radiation wavelength. In the case of obliquely incident plane waves, we reveal the limitation on the loss/gain coefficient, which should not be too large compared with the real part of permittivity. We believe, our findings can improve fundamental understanding of physics behind PT-symmetric systems and make a new step in development of auxiliary tools for analyzing their peculiar optical response.


I. INTRODUCTION
PT symmetry entered physics as a realization of the non-Hermitian quantum mechanics keeping eigenvalues real [1,2]. It soon turned out that the properties of PT -symmetric systems can be relatively easily proved in practice in the optical [3,4] and microwave [5] domains (see also recent review articles [6,7]). Generally, a PT -symmetric optical system is a periodic structure with balanced loss and gain (the same value of loss and gain coefficients): the field of the eigenmodes is equally distributed between the loss and gain components. When the loss/gain coefficient increases, the balance may spontaneously disappear at the point of symmetry breaking (exceptional point), where the eigenstates become degenerate. In the non-PT -symmetric regime, the field is asymmetrically distributed between the loss and gain parts of the system producing amplifying and decaying eigenmodes. The system may return to the PT -symmetric state at the following exceptional point (see, e.g., [8]).
Optical PT symmetry is the basic concept for various prospective applications including lasing [9,10] and coherent perfect absorption [11,12], enhanced sensing [13][14][15], effects of asymmetric light propagation such as unidirectional invisibility [16], nonreciprocity [17], localization [18], and so on. Optical PT -symmetric systems have usually either waveguide or multilayer configuration, although there are more exotic variants, such as metasurfaces [19,20] or graphene-based structures [21]. PT waveguides can be described in the coupled-mode approximation and can be both active [22] and passive [23]. Multilayer structure as an important system to understand basics of the PT symmetry without using any ap- * dvnovitsky@gmail.com † andreyvnovitsky@gmail.com proximations will be analyzed in this paper.
A multilayer can be considered as a photonic crystal or a simplest metamaterial (artificial periodic subwavelength structure) [24][25][26] depending on relation between the size of the unit cell and the radiation wavelength. Key problem of the metamaterial theory is a homogenization that allows us to treat the metamaterial as a quasicontinuous medium with a set of effective parameters, such as dielectric permittivity and magnetic permeability. Metamaterials of the multilayer geometry can be homogenized using the standard techniques, such as first-principle homogenization [27,28], nonlocal effective-medium theory [29][30][31], Whitney interpolation [32], etc. However, to the best of our knowledge, there has been no thorough investigation of the homogenization of PT -symmetric systems before. In the most relevant article [33], the PT -symmetric system is composed of the hyperbolic metamaterial and gain medium, the former being described using the Maxwell Garnett approach.
In this paper, we employ the operator effective medium approximation (OEMA) [29,31] to investigate its area of validity in description of PT -symmetric multilayered systems. OEMA juxtaposes a homogeneous nonlocal bianisotropic effective medium to the multilayer, thus allowing us to accurately find the transmission/reflection spectra [31] and surface-wave propagation [34]. Neglecting the nonlocal effects, the OEMA provides the Maxwell Garnett approximation and the well-known mixing formulae for the components of effective permittivity tensor [35]. In Section II, we write out the effective medium tensors derived in Ref. [31] for the PT -symmetric multilayer in the zeroth, first, and second orders of the OEMA. In Section III, we discuss the transmission/reflection characteristics and description of PT symmetry breaking using the local and nonlocal material parameters and find out the limits of their applicability. Section IV sums up the main results of the article.

II. OPERATOR EFFECTIVE MEDIUM APPROXIMATION FOR A PT -SYMMETRIC MULTILAYER SYSTEM
We consider a PT -symmetric multilayered structure illustrated in Fig. 1. It consists of alternating loss and gain layers of the same thickness d/2, the total number of layers being 2N. To ensure the optical PT symmetry, the permittivity in the multilayer should be distributed as ε(z) = ε * (z), that is realized using the loss ε L = ε ′ + iε ′′ and gain ε G = ε ′ − iε ′′ permittivities (loss/gain coefficient ε ′′ > 0). Such a simple structure can be fully described using the transfer-matrix method (TMM) including its transmission and reflection (scattering) properties. We will exploit the TMM for homogenized slabs as well.
To homogenize the multilayered system, we use a recently developed operator effective medium approximation [29,31]. The idea behind the OEMA is based on writing the fundamental solution for the layered periodic structure as a series over the size parameter k 0 d, where k 0 is the vacuum wavenumber and d is the thickness of the unit cell. Then one is able to introduce the m-th order of approximation as a solution containing the terms up to (k 0 d) m . The larger the order, the better the approximation. We do not derive equations for the effective material parameters here, but just borrow them from the previous research [31]. Let us start with the zeroth-order (m = 0) approximation known as the Maxwell Garnett approximation.

A. Maxwell Garnett approximation
Maxwell Garnett approach is believed to be valid, if k 0 d ≪ 1. In this case, the loss/gain multilayer can be represented as a homogeneous uniaxial medium characterized by the permit-tivity tensorε Since the size parameter enters the above expression as (k 0 d) 0 , the Maxwell Garnett approach is the zeroth-order approximation. It should be stressed that this oversimplified technique may not be satisfactory, even if k 0 d ≪ 1 (see the works on the breakdown of the effective medium theory, e.g., [36,37]).

B. First-order OEMA
Material parameters in the first-order approximation imply nonlocality and bianisotropy. Bianisotropic (magnetoelectric) terms emerge in the first order, (k 0 d) 1 , together with the permittivity tensor (1). Magnetoelectric coupling tensor was derived in [31] and for the loss/gain multilayered system equalŝ Here the nonlocality appears as evidenced by the dependence of the material parameters on the tangential wavenumber k t , i.e. the projection of the wavevector on the interface between slabs. The constitutive equations for such a bianisotropic medium read D =ε (0) E+αH and B = H +α T E, where E and H (D and B) are the electric and magnetic field strengths (inductions) and superscript T stands for the transposition. It is worth noticing that the magnetoelectric coupling (2) depends on the wave propagation direction, that is the order of layers (ε L ↔ ε G is equivalent to ε ′′ ↔ −ε ′′ ).

C. Second-order OEMA
The second-order corrections influence the effective permittivity and permeability tensors keeping the effective magnetoelectric coupling tensors as it has been defined according to Eq. (2). According to Ref. [31], the permittivity and per-meability tensors for the loss/gain multilayer take the form The permittivity correction is electric quadrupolar and nonlocal (k t -dependent). The permeability is caused by the artificial magnetic moment due to the displacement currents. The effective medium in the second-order OEMA is characterized by the constitutive equations D =ε (2) E +αH and B = µ (2) H +α T E.

III. WAVE PROPAGATION IN HOMOGENIZED
PT -SYMMETRIC SLAB

A. Normal incidence
Here we compare the scattering properties of homogenized versus multilayered PT -symmetric media. Strictly speaking, the homogenized slab cannot be treated as the PTsymmetric system, but a wave propagation in both systems can be quite similar. Homogenized material is considered as anisotropic (Maxwell Garnett approximation) or nonlocal bianisotropic (first and second order of OEMA) media. In Fig. 2, we demonstrate the transmission T of a normally incident plane wave. In the zeroth order OEMA, the transmission of the normally incident wave is defined by the in-plane permittivity ε (0) || . That is why it does not depend on the imaginary part of permittivity ε ′′ [see the horizontal straight line in Fig. 2(a)]. It should be noticed that the Maxwell Garnett approach reproduces the correct value of T only for small ε ′′ , when the second-order corrections are negligible. The region of validity of the first-order OEMA is wider, but it is anyway very limited. In the second order, the transmission spectra well follow the curve corresponding to the transmission through the multilayer. However, this approximation is valid only for the thin unit cells like d = 100nm ≪ λ . For instance, when we double the thickness of the unit cell, the second order OEMA only qualitatively describes the transmission reproducing peaks, but not their positions and heights [see Fig.  2 Since the thicknesses of the layers are limited by the validity of the effective medium theory, the large values of the loss/gain coefficient ε ′′ should be used for reproducing the symmetry breaking in PT -symmetric systems, when loss and gain are no longer balanced. We exploit the standard means for studying the PT symmetry breakdown, i.e., the calculation of eigenvalues and eigenvectors of system's scattering matrix. The scattering matrix for the multilayered sys- tem and the homogenized medium under consideration can be derived from the corresponding transfer matrix [38]. The eigenvalues of the scattering matrix are unitary (|s 1 | = |s 2 | = 1) in the PT -symmetric state and inverse (|s 1 | = 1/|s 2 | > 1) in the PT -broken state. The eigenvectors of the scattering matrix coalesce at the exceptional point (point of the non-Hermitian singularity), where a phase transition from the PT -symmetric to the non-PT-symmetric state occurs.
Let us examine, whether the local or nonlocal homogenization is able to catch the exceptional points of the PT symmetry breaking. Scattering matrix, transmission, and reflection are expected to be simultaneously well predicted, since the scattering matrix is defined in terms of the amplitude transmission and reflection coefficients. In Fig. 3(a), we show dependence of the eigenvalues of the scattering matrix on ε ′′ for d = 100 nm. Similar to the transmission curves reported in Fig. 2(a), the second-order OEMA correctly describes behavior of the eigenvalues and reproduces the breaking of the PT symmetry observed in the TMM calculations of the multilayer. The second-order OEMA well suits also for describing the coalescence of the eigenvectors (not shown here). Detailed analysis shows that neither the first-order nor zeroth-order OEMA give a hint of the phase transition. The range of applicability for the Maxwell Garnett approximation is limited by the small loss/gain coefficients ε ′′ , where the PT symmetry breaking does not occur. The symmetry breaking could be reached for smaller ε ′′ , if the unit cells were thicker. However, the larger d would also ruin the Maxwell Garnett approximation. It is worth mentioning that the homogeneous slab characterized by the permittivity tensor Eq. (1) is inherently inappropriate for description of the exceptional points, because its transmission and reflection properties do not depend on the direction of wave incidence. For instance, in the case of the normal incidence, the properties of the homogeneous slab are equivalent to those of an isotropic dielectric slab of permittivity ε ′ .
The breakdown of the PT symmetry cannot be also obtained, if we skip the terms proportional to (k 0 d) 1 and leave the second-order terms (k 0 d) 2 . This means that only the mix of the first-order and second-order terms properly reproduces  the eigenvalues. As previously, the thicker layers ruin the correspondence between the exact and homogenized description of the structure [ Fig. 3(b)]. The difference between the accurate and approximate eigenvalues is mainly the shift along the ε ′′ -axis. The eigenvalues in the second-order OEMA s (2) are roughly related to the eigenvalues of the PT -symmetric multilayer s as s(ε ′′ ) ≈ s (2) (ε ′′ + ak 0 d), where the shift is proportional to the size parameter k 0 d and a is a constant. Hence, we can estimate the higher-order corrections caused by exci- tation of the higher-order multipoles as [ds (2) (ε ′′ )/dε ′′ ]ak 0 d.

B. Oblique incidence
Now we face to the dependence of the multilayer response on the incidence angle of the plane waves. In Fig. 4, we observe a strong dependence of the correctness of the homogenization on the value of ε ′′ . For ε ′′ = 1, both transmission and reflection dependencies on the incidence angle are reliably described by the second-order OEMA, but this is not already the case at ε ′′ = 5. This discrepancy grows with the increase of ε ′′ : the second-order OEMA may even predict a number of artifact peaks absent in the curves calculated for the multilayer.
Validity of the nonlocal effective medium approximation of the second order is determined by the loss/gain coefficient ε ′′ , but not the angle of incidence. We relate the behavior shown in Fig. 4(b) to the appearance of the PT -symmetrybroken states for the larger incidence angles. At the oblique incidence, the light passes a longer path compared to the normal incidence, therefore, the effective thickness of the structure is enlarged and exceptional points of the PT symmetry breaking may appear at the smaller ε ′′ . In Fig. 5 we plot the angular dependencies of the scattering matrix eigenvalues for ε ′′ = 10 and ε ′′ = 20. The results for the multilayered and homogenized structures are qualitatively different in these cases. The OEMA cannot catch the PT -symmetry-broken state for ε ′′ = 10, because it does not properly enlarge the effective thickness at the oblique incidence. Indeed, the effect of the tangential wavenumber k t is negligible for the great ε ′′ [see Eqs. (2) and (3)], i.e., the results are quite close to those for k t = 0. As well as for the normal incidence, the PTsymmetry-broken states appear at the larger ε ′′ = 20, but they do not correspond to the accurate calculations.
According to Eqs. (2) and (3), the wavenumber k t makes significant effect for greater real part of the permittivity ε ′ . Figure 6 shows the angular dependencies of transmission and reflection at ε ′ = 10 and small loss/gain coefficient ε ′′ = 1. Comparing with the good matching of the accurate and approximate results for ε ′ = 2 and ε ′′ = 1 (see Fig. 4), one can note that there is some discrepancy between the OEMA and accurate calculations, but only in the region of large incidence angles. We conclude that the nonlocal homogenization may be applicable, if ε ′′ is not much greater than ε ′ .

IV. CONCLUSION
We have studied the local and nonlocal homogenization of the PT-symmetric multilayered system. It has been carried out using the previously developed operator effective medium approximation providing the successive approximations for effective parameters including the usual mixing (Maxwell Garnett) formulae and the nonlocal bianisotropic material tensors. We have found that the Maxwell Garnett approach is entirely inappropriate for description of the PT symmetry. The nonlocal model takes into account the distribution of the loss and gain materials and can be applied in the limit of electrodynamically small thicknesses of the unit cells. Such a behavior is equally related to the transmission and scattering matrix' eigenvalues spectra. For the obliquely incident plane waves, the validity of the nonlocal homogenization is limited by the loss/gain coefficients ε ′′ comparable or smaller than the real part of the permittivity ε ′ . Thus, though the nonlocal homogenization is applicable, it is strongly restricted by the thickness of the unit cell and the value of the loss/gain coefficient.