Inverse dispersion method for calculation of complex photonic band diagram and $\cal{PT}$-symmetry
Abstract
We suggest an inverse dispersion method for calculating photonic band diagram for materials with arbitrary frequency-dependent dielectric functions. The method is able to calculate the complex wave vector for a given frequency by solving the eigenvalue problem with a non-Hermitian operator. The analogy with -symmetric Hamiltonians reveals that the operator corresponds to the momentum as a physical quantity and the singularities at the band edges are related to the branch points and responses for the features on the band edges. The method is realized using plane wave expansion technique for two-dimensional periodical structure in the case of TE- and TM-polarization. We illustrate the applicability of the method by calculation of the photonic band diagrams of an infinite two-dimension square lattice composed of dielectric cylinders using the measured frequency dependent dielectric functions of different materials (amorphous hydrogenated carbon, silicon, and chalcogenide glass). We show that the method allows to distinguish unambiguously between Bragg and Mie gaps in the spectra.
I Introduction
In the past decades, a significant progress has been achieved in the fabrication of light-controlling nanostructures including metamaterials (1); (2) and photonic crystals (3); (4). Hence the ability to accurately calculate photonic properties is required for designing novel devices, while the photonic band diagram, i.e., frequency vs. wave vector , being the most important inherent property of periodic media (3). The common approach to the problem consists in reducing Maxwell’s equations to the standard eigenproblem for the frequency as an eigenvalue (5); (6); (7); (8); (9). The first successful approach (5) assumed that the dielectric function does not vary with frequency. The authors used the plane wave expansion technique to solve the wave equation with a Hermitian operator for the magnetic field. The subsequent papers suggested numerical improvements such as application of the fast Fourier transform method, simultaneous calculations of a block of eigenvectors, etc (10); (8). In the case of a cylindrical or spherical symmetry of the scatterer the Green’s function method was proposed (11); (6); (12). With these methods the photonic band diagram of the most important structures as well as the conditions for the appearance of a complete photonic band gap were determined (7); (13); (14).
At the same time, most materials strongly interacting with light have a frequency dependence of dielectric function and accurate methods should take this circumstance into account. A number of approaches suggesting that is a Lorentz-like single pole function have been proposed (15); (16); (17); (18); (19); (20); (21); (22); (23). The following techniques have been used to calculate the photonic band diagram: excluding the frequency from the operator to yield a generalized eigenproblem (16); considering the polarization field and solving it in the time domain (19); including the mechanical problem for charged particles along with the electric and magnetic fields into the operator (23). However, the optical range of electronic interband transition is not described correctly by the Lorentz-type functions and other advanced methods should be employed. In addition, the traditional methods calculate only propagating modes with a real wave vector. They cannot provide any information about the light propagation within the bandgap frequencies, and about the rate at which the evanescent mode disappears. Also a study of modes within the forbidden band is important for understanding the transition between two regimes of light propagation in a dielectric periodical structure, i.e. the regime where the propagation of light is mainly influenced by the interference between multiple scattered Bragg waves (photonic crystals), and the regime where the propagation of light is mainly determined by the properties of each single element of the materials (Mie modes of metamaterials). The transition is accompanied by complicated interactions of multiple bands and the correct interpretation is a key to creating a phase diagram photonic crystal/all-dielectric metamaterial (24). Therefore we need a tool to explore evanescent modes within a bandgap and distinguish whether it is the Bragg gap related to the non-localized Bragg resonance or the Mie gap originated from a localized Mie-like resonance.
An alternative to the conventional techniques is to consider methods for calculating complex wave vector while the frequency is treated as parameter. We can find it from transcendent KKR equation (Green’s function method) (11); (6); (12) using so called on-shell approach (the term is borrowed from scattering theory that means solutions on the energy shell (25)). Another method is based on Wannier functions (26). Using these methods we eliminate two limitation of the direct approaches. The first one is related to the strong frequency dependence of . For traditional numerical techniques to be applicable, should be excluded from the operator but this can only be done for restricted classes of function. The second limitation is that the direct methods are intended for calculating the propagating modes. By contrast, in implementation of the inverse method the wave vector should be excluded from the operator. Fortunately, in most cases the spatial dispersion is negligible and the dielectric function is independent of the wave vector considered to be the eigenvalue. Also we can study prohibited evanescent modes by setting within the bandgap and finding whether or not the mode originates from the localized resonance and the structure can be homogenized (27); (28); (29). However both Green’s function method(11); (6); (12) and the Wannier method (26) are quite complicated techniques relative to straightforward plane wave expansion techniques. On the other hand plane wave based methods (so-called rigorous coupled wave analysis) are widely used for calculating transport properties of periodic slab by fixing tangential component of the wave vector and the frequency (30); (31); (32); (33); (34).
In this paper we propose an alternative to the conventional methods by introducing an inverse dispersion method for calculating the photonic band diagram. We reduce Maxwell’s equations to a problem with the eigenvalue , while is considered to be a real parameter. We name this method the inverse problem instead of the direct problem. These method has straightforward implementations in plane wave basis as well as in real space. By way of illustration, we consider three periodic structures with the same geometry but different materials, namely amorphous hydrogenated carbon a-C:H (low dielectric contrast with air in VIS), silicon ( in VIS), and chalcogenide glass GeSbTe (high dielectric contrast in IR). We note that the 2D square lattice of dielectric cylinders acts as a photonic crystal at a low dielectric contrast (35), whereas in a special structural condition the same structure with a high enough contrast behaves as a metamaterial with (36). In the present study we use dielectric function measured in Refs. (37); (38); (39); (40) to calculate the photonic band diagram.
Ii Inverse problem
In this section we formulate the electrodynamic equations in the form suitable for the inverse dispersion approach. In general, the electromagnetic field is described by the Maxwell’s equations which for the frequency harmonics read
(1a) | |||
(1b) |
Here and are respectively the electric and magnetic vectors, and are the dielectric permittivity and the magnetic permeability. In periodic structures the solutions can be written as the Bloch waves and . By substituting these into Eqs. (1) we obtain the generalized eigenproblem for the frequency
(2) |
Instead of applying the operator at r.h.s. once again to get the direct problem we extract the wave vector from the operator at l.h.s. and move the matrix with the frequency to l.h.s. Now, we obtained generalized eigenproblem for the wave vector
(3) |
where is the direction of the wave vector and .
It is known that the direct eigenproblem corresponds to the Schrödinger equation (3). For the loss-free photonic structure the operator is Hermitian and all eigenfrequencies are real. However, in the past several years non-Hermitian -symmetric Hamiltonians have caught a lot of attention (41). An interesting property of -symmetric operators is the transition between the completely real spectrum and the complex spectrum, which are called exact/spontaneously-broken -symmetry. The common situation of the transition is the coalescence of two eigenvalues at a branch point, which is referred to as an exceptional point in the operator theory (42). If we consider a finite dimensional space and a generalized eigenproblem with a perturbed operator in the form the eigenvalues satisfy the characteristic equation having the form of an algebraic equation in or . Let us return to Eq. (3). For a material with a real dielectric constant we suppose the wave vector to be a given real parameter and write the direct problem with real . But now we consider to be a given parameter. In the periodic structure for which falls within the allowed frequency bands we have a real , whereas for the forbidden bands the wave vector becomes complex. Hence we can regard the bandgap boundaries as the branch points and there is an analogy between the band theory and the problem of a spontaneously broken -symmetry.
Let us describe a complex band diagram that accounts for both the propagating and evanescent modes. We start from the empty-lattice approximation (i.e., consider homogeneous media) to give a branch notation. For simplicity we consider a two dimensional reciprocal space. For the constant and there is a conical relation . However the periodicity results in that this light cone has many repetitions with the origin at any reciprocal lattice vector [see. Fig. 1(a)]. The Bragg resonances that lead to the scattering from one cone to another arise in vicinity of the intersection of the cones, whereas outside the intersection we can identify the cones by indices. As an example we consider a scanning plane in the direction of the Brillouin zone of the square lattice. Only the cones with the origin lying along this direction intersect the scanning plane at the low frequency range [Fig. 1(b)]. This is a common situation for the direct method. If, however, we use the inverse method we obtain the evanescent modes because the cone apex is out of the scanning plane and the conical relation gives complex solutions. However in the empty-lattice approximation all the modes are orthogonal and there is no coupling between any modes, and hence the evanescent modes can be neglected, but for the case of a nonzero dielectric contrast we cannot do this anymore. Note that in this paper we plot band diagrams in the scale of the dimensionless frequency , where is the wavelength in a vacuum and is the lattice constant. Below we consider the case of the two-dimensional structure however the general case is quite similar.
ii.1 TE-polarization
In the TE-polarization we deal with the -component of the magnetic Bloch’s wave , as well as with the tangential and the normal components of the electric Bloch’s wave. Equation (3) can be rewritten as
(4a) | |||||
(4b) | |||||
(4c) |
Note that the functions , and are periodic. To proceed further we have to reduce these equations to those of the algebraic type. There are two relatively simple abilities. The first one is to represent operators in plane wave basis, and the second is to use a finite difference scheme in real space. The latter ability is simpler in realization, however it is not work perfectly if the calculating grid is too rare being pure numerical technique In contrast plane wave basis makes it possible to evaluate mode properties within two-wave approximations if the modulation of dielectric function is moderate. On the other hand if we consider periodic structure with metal components, plane wave basis is not an optimal choose (31); (32), since zero-order approximation (approximation of free-photons) cannot describe the system correctly, because of free-photon averages dielectric function: positive air with negative dielectric index of metallic inclusions. As a result averaged dielectric index is negative and propagation of light is prohibited that is completely disagree with the correct answer. Hence, to calculate structures with metallic inclusions we should take into account a very large number of plane waves and this technique is worse than finite difference scheme. However, even if we deal with a pure dielectric structure that possesses strong dielectric contrast we have to consider convergence problems as well (43), and the usage of only few number of plane waves is inappropriate. To be sure in correctnesses of our study we examine convergence of traditional plane wave method and find that 625 (25 by 25) plane waves are sufficient to obtain reliable results.
Because the functions , and are proportional to the Bloch’s amplitudes they can be represented as Fourier series. We expand , and in the r.h.s., find Fourier amplitudes of both sides of Eqs. (4), and finally write these equations in the matrix form
(5) |
(6) |
(7) |
Here and are diagonal matrices with elements and , and are Toeplitz matrices with elements , and , which are the Fourier amplitudes of and respectively. On eliminating the component , we obtain the conventional eigenproblem
(8) |
ii.2 TM-polarization
Following similar procedure we write linear eigenproblem for the electromagnetic field in TM-polarization
(9) |
Here and define tangential component of magnetic field and -component of electric field. Matrices , , , and was defined above in TE-polarization subsection. In most cases material does no posses magnetic properties and is the identity matrix. Hence Eq. (9) is somewhat simpler than Eq. (8) since in case of TM polarization the matrix does not contain the inverse of .
Iii Band diagram for photonic crystals and dielectric metamaterials
To illustrate the inverse dispersion method we consider here three types of structures with the same geometry. Periodic all-dielectric systems are known to be in two scattering regime either as a photonic crystal with the Bragg fundamental gap or as a metamaterial with the Mie fundamental gap (24). By the fundamental gap we mean the gap between the first and second photonic bands (3). Here we consider 2D square lattices of parallel dielectric circular rods infinitely long in the direction in air with lattice constant and radius of . We study three different dielectric functions for the rods corresponding to a-C:H, Si, and GeSbTe.
iii.1 Constant dielectric function
As a reference point we calculated the photonic band diagram by the traditional method in the implementation of Ref. (9) by using constant values for the dielectric function: , , and (the thin magenta dash-and-dot curves in Fig. 2). There are bandgaps within the frequency spectra and we have no information about the evanescent modes. However these modes define the scattering, transmission and other light transport properties. Next, we calculate the photonic band diagram by the inverse dispersion method which yields both propagating and evanescent modes (Fig. 2). The bands for the propagating modes (with real wave vector) obtained by using and demonstrate a complete agreement, but now we obtain the dispersion for the evanescent modes. First we consider TE-polarization.
TE-polarization
The spatial modulation of the dielectric function perturbs eigenmodes calculated in the empty-lattice approximation. If the perturbation is small enough (photonic crystal regime) we can identify the branches as and . Their interaction produces a pair of branch points at the band gap edges marked by the green circles in Fig. 2(a,d). We describe the frequency dependence of the branches and . In the low frequency range the wave vector is real. In the interval defined by a pair of the branch points the real part of the wave vector is constant while the imaginary part varies in amplitude, with of about 0.1. As a result the branches form an narrow elliptical feature (NEF) that is prolate along the frequency axis. The feature corresponds to the bandgap. With increasing dielectric contrast, the higher frequency modes cannot be adequately identified with the modes due to the strong interband interaction. In the intermediate case of [Fig. 2(b)] we can emphasize three NEF: one centered at [, ] and a pair of similar NEF at [, ], which do not intersect the real plane . In the metamaterial regime(24) [Fig. 2(c)] we recognize the NEF at [, ]. In addition, we observe two wide features of complicated shape (WCF). The first WCF is bounded by with a width , and the second lies at with a width . We notice that all the NEF and WCF in Fig. 2 originate at the surface of the Brillouin zone (X point). Analysis of the band diagrams reveals that the wave vector of the branch may vary in its real or imaginary part only and these two types of behavior are switched by the branch point. In addition, the NEF keeps modes connected while the WCF splits the modes to be disconnected at the X point (i.e, we cannot join points at the lower and upper bands).
TM-polarization
Now we describe band diagrams for TM-polarization. When the dielectric contrast is small [Fig. 2(d)] the band diagram is similar to the case of TE-polarization [Fig. 2(a)]. We identify NEF with parameters and being wider than in Fig. 2(a)]. Now let us consider band diagram for , which differs from TE-polarization. Figure 2(e) demonstrates just one NEF shifts down and increases its width and a WCF ( and ). With increasing of dielectric index up to 30 the NEF moves down to the interval () and does not disappear (as in TE-polarization) remaining fundamental gap. Also we identify another NEF ( and ) and a pair of slightly overlapped WCF (, ) and (, ). Being similar to the case of TE-polarization the NEF keeps modes connected while the WCF splits the modes to be disconnected at the X point (i.e, we cannot join points at the lower and upper bands).
iii.2 Experimental dielectric function
We can calculate photonic band diagrams for the three systems under consideration by using the experimentally measured dielectric response of a material (37); (38); (39); (40) plotted in Fig. 3, since the inverse dispersion method is applicable for calculation of an arbitrary dielectric function. As far as Maxwell’s equations are not frequency-scalable for non-constant dielectric functions we have to specify both the material and the geometrical sizes of the 2D structures: (i) for a-C:H the lattice constant is set to nm and the dielectric function is taken from Ref. (39); (ii) for Si the lattice constant is set to nm, the dielectric function is taken from Refs. (37); (38); and (iii) nm for a metamaterial fabricated form GeSbTe (fcc phase of a GeSbTe chalcogenide glass)(40). Calculated photonic band diagrams are presented in Fig. 4 and comparison with Fig. 2 reveals a number of essential differences in the mode behavior. The imaginary part of the dielectric function lifts the degeneracy at the branch points marked by the green circles in Fig. 2. It is convenient to introduce the density of states per interval of a real part of wave vector for the branches. In the case of real dielectric function both types of features NEF and WCF lie at the surface of the Brillouin zone [] with the density of states becoming a delta-function at the X point. With increasing imaginary part of the dielectric function, the density of states broadens away from the surface of the Brillouin zone.
Also, there is a difference in the mode behavior at NEF and WCF when it is fundamental gap. We notice that WCF becomes fundamental gap only in TE-polarization. Figure 5 shows the band gaps related to different features NEF and WCF on a detailed scale. The group velocity does not change its direction at NEF, while at WCF the group velocity becomes zero and changes its direction to the opposite. By contrast, in the case of a real dielectric function there is an interval of a vanished group velocity hence the dielectric function needs to possess at least a small amount of an imaginary part to distinguish the type of the gap.
Iv Discussion
The inverse dispersion method is applicable for accurate calculating of the photonic band diagram for periodic structures composed of materials with an arbitrary dielectric response whereas the conventional techniques cannot describe the photonic properties correctly. Also the method provides a deep insight into the physics of propagation of traveling of electromagnetic waves in periodic media including photonic crystals and metamaterials. The photonic band diagram calculated by the direct method hides the branch points because this method is inapplicable to finding of the evanescent modes. At the same time the linear operator theory describes singular properties of eigenmodes at the branch point (42).
In the imaginary plane of the band diagram, the NEF is unambiguously associated with the Bragg resonance, because it increases at very small dielectric contrast and involves a pair of modes [Fig. 5(b)]. The wide complicated feature is associated with the Mie resonance because of being related to the Mie gap in all-dielectric structures as shown in Ref. (36). The complex photonic band diagram in Fig. 2(c) demonstrates that the Mie resonance feature extends from to (see in the imaginary plane) being much wider than the Mie-band gap (thin dash-and-dot curves in the real plane). Analysis of the frequency ranges of the resonant Mie modes (44); (45) shows that the lower Mie feature in Fig. 2(c) is associated with the TE Mie mode while the higher feature is relater to the TE mode. In a similar way WCF in Fig. 2(c) are originated from the TM and TM Mie modes. We also demonstrate that the Mie-related WCF make the band diagram disconnected at the X point [Fig. 2(c)].
Turning to the Eqs. (8) and (9), we notice that they are non-Hermitian operators with the wave vector length as an eigenvalue. Therefore we can consider these operators to be a wave vector operator for the periodic media. Inspired by the results obtained from non-Hermitian -symmetric Hamiltonians (41) we assume that the wave vector operator possesses a certain symmetry. The symmetry is related to the time-reversal symmetry, while the parity should be applied to the fields represented in the -space (we notice that the real space representation is no better than the -space representation). The important consequence of the description of operator (8) and (9) as operators with -symmetry is that we should not consider these operators as a mathematical trick only and reject they as being unphysical because they violate the axiom of Hermiticity. Conversely we do consider these operators as physical quantities that are related to momentum of excitation in periodical media. At the branch point the symmetry changes from exact to spontaneously-broken that makes this excitation prohibited.
In addition, using the results of the inverse dispersion method we can separate the Bragg and Mie gaps by introducing a small amount of loss. The modes forming Bragg gaps do not change their direction, while modes forming the Mie gaps show a repulsive behavior. Using this criteria we explicitly identify the lower bandgap as a Bragg gap (Fig. 5).
V Conclusion
We have suggested the inverse dispersion method for calculating the complex photonic band diagram for periodic media with an arbitrary dielectric function. The method makes it possible to calculate both the propagating and evanescent modes. We have calculated the band diagram of periodic structures composed of the following materials: low-refractive index amorphous hydrogenated carbon (constituting a photonic crystal), high-refractive index silicon, and chalcogenide glass (constituting an all-dielectric metamaterial) by using the experimentally measured properties of the material. Our analysis demonstrated the essential difference between the fundamental gaps that originates from Bragg or Mie evanescent modes. This result is of fundamental importance for studies of the light scattering regimes of photonic crystals and metamaterials(24). Moreover we note that the inverse dispersion method is not restricted to only photonic diagrams, being rather applicable to any waves in periodic media: electrons, phonons, magnons etc.
Acknowledgements.
We acknowledge fruitful discussions with A.A. Kaplyanskii, Yu.S. Kivshar, S.F. Mingaleev, and K.B. Samusev. This work was supported by the Russian Foundation for Basic Research (15-02-07529).References
- N. Enghata and R. Ziolkowski, eds., Electromagnetic Metamaterials: Physics and Engineering Exploration (Wiley-IEEE Press, New York, 2006).
- D. Smith and R. Liu, eds., Metamaterials: theory, design, and applications (Springer, New York, 2010).
- J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light (Princeton Univ. Press, 2008), 2nd ed.
- M. F. Limonov and R. M. De La Rue, eds., Optical properties of photonic structures: interplay of order and disorder (CRC Press, Taylor & Francis Group, 2012).
- K. M. Ho, C. T. Chan, and C. M. Soukoulis, Phys. Rev. Lett. 65, 3152 (1990).
- X. Wang, X.-G. Zhang, Q. Yu, and B. N. Harmon, Phys. Rev. B 47, 4161 (1993).
- K. Busch and S. John, Phys. Rev. E 58, 3896 (1998).
- S. G. Johnson and J. D. Joannopoulos, Opt. Express 8, 173 (2001).
- A. V. Moroz, M. F. Limonov, M. V. Rybin, and K. B. Samusev, Phys. Solid State 53, 1105 (2011).
- R. D. Meade, A. M. Rappe, K. D. Brommer, J. D. Joannopoulos, and O. L. Alerhand, Phys. Rev. B 48, 8434 (1993).
- K. M. Leung and Y. Qiu, Phys. Rev. B 48, 7767 (1993).
- A. Moroz, Phys. Rev. B 51, 2068 (1995).
- A. Moroz, Phys. Rev. Lett. 83, 5274 (1999).
- M. V. Rybin, I. I. Shishkin, K. B. Samusev, P. A. Belov, Y. S. Kivshar, R. V. Kiyan, B. N. Chichkov, and M. F. Limonov, Crystals 5, 61 (2015).
- J. B. Pendry, J. Phys.: Cond. Matt. 8, 1085 (1996).
- V. Kuzmiak, A. A. Maradudin, and A. R. McGurn, Phys. Rev. B 55, 4298 (1997).
- V. Kuzmiak and A. A. Maradudin, Phys. Rev. B 58, 7230 (1998).
- P. Halevi, , and F. Ramos-Mendieta, Phys. Rev. Lett. 85, 1875 (2000).
- K. Sakoda, N. Kawai, T. Ito, A. Chutinan, S. Noda, T. Mitsuyu, and K. Hirao, Phys. Rev. B 64, 045116 (2001).
- E. Moreno, D. Erni, and C. Hafner, Phys. Rev. B 65, 155120 (2002).
- K. C. Huang, P. Bienstman, J. D. Joannopoulos, K. A. Nelson, and S. Fan, Phys. Rev. Lett. 90, 196402 (2003).
- O. Toader and S. John, Phys. Rev. E 70, 046605 (2004).
- A. Raman and S. Fan, Phys. Rev. Lett. 104, 087401 (2010).
- M. V. Rybin, D. S. Filonov, K. B. Samusev, P. A. Belov, Y. S. Kivshar, and M. F. Limonov, Nature Commun. 6, 10102 (2015).
- R. G. Newton, Scattering theory of waves and particles, International series in pure and applied physics (McGraw-Hill, 1966).
- D. Hermann, M. Diem, S. F. Mingaleev, A. García-Martín, P. Wölfle, and K. Busch, Phys. Rev. B 77, 035112 (2008).
- S. Datta, C. T. Chan, K. M. Ho, and C. M. Soukoulis, Phys. Rev. B 48, 14936 (1993).
- A. Kirchner, K. Busch, and C. M. Soukoulis, Phys. Rev. B 57, 277 (1998).
- L. Tartar, The general theory of homogenization (Springer, 2009).
- M. Moharam and T. K. Gaylord, JOSA 71, 811 (1981).
- L. Li and C. W. Haggans, JOSA A 10, 1184 (1993).
- P. Lalanne and G. M. Morris, JOSA A 13, 779 (1996).
- S. Shi, C. Chen, and D. W. Prather, Appl. Phys. Lett. 86, 043104 (2005).
- I. I. Shishkin, M. V. Rybin, K. B. Samusev, V. G. Golubev, and M. F. Limonov, Phys. Rev. B 89, 035124 (2014).
- F. Quiñónez, J. W. Menezes, L. Cescato, V. F. Rodriguez-Esquerre, H. Hernandez-Figueroa, and R. D. Mansano, Opt. Express 14, 4873 (2006).
- S. O’Brien and J. B. Pendry, J. Phys.: Cond. Matt. 14, 4035 (2002).
- H. H. Li, J. Phys. Chem. Ref. Data 9, 561 (1980).
- D. E. Aspnes and A. A. Studna, Phys. Rev. B 27, 985 (1983).
- E. D. Palik, Handbook of optical constants of solids, vol. 3 (Academic press, 1998).
- P. Němec, A. Moreac, V. Nazabal, M. Pavlišta, J. Přikryl, and M. Frumar, J. Appl. Phys. 106, 103509 (2009).
- C. M. Bender, Rep. Prog. Phys. 70, 947 (2007).
- T. Kato, Perturbation theory for linear operators (Springer, 1966).
- H. S. Sözüer, J. W. Haus, and R. Inguva, Phys. Rev. B 45, 13962 (1992).
- M. V. Rybin, K. B. Samusev, I. S. Sinev, G. Semouchkin, E. Semouchkina, Y. S. Kivshar, and M. F. Limonov, Opt. Express 21, 30107 ( 2013).
- M. V. Rybin, D. S. Filonov, P. A. Belov, Y. S. Kivshar, and M. F. Limonov, Sci. Rep. 5, 8774 (2015).