Aperiodic-Order-Induced Enhancement of Weak Nonlocality in Multilayered Dielectric Metamaterials
AAperiodic-Order-Induced Enhancement of Weak Nonlocality inMultilayered Dielectric Metamaterials
Marino Coppolaro, Giuseppe Castaldi, and Vincenzo Galdi ∗ Fields & Waves Lab, Department of Engineering,University of Sannio, I-82100 Benevento, Italy (Dated: November 6, 2018)
Abstract
Recent studies on fully dielectric multilayered metamaterials have shown that the negligibly small nonlocal effects (spatial dispersion) typically observed in the limit of deeply subwavelength layersmay be significantly enhanced by peculiar boundary effects occurring in certain critical parame-ter regimes. These phenomena, observed so far in periodic and randomly disordered geometries,are manifested as strong differences between the exact optical response of finite-size metamaterialsamples and the prediction from conventional effective-theory-medium models based on mixingformulae. Here, with specific focus on the Thue-Morse geometry, we make a first step towardextending the studies above to the middle-ground of aperiodically ordered multilayers, lying in be-tween perfect periodicity and disorder. We show that, also for these geometries, there exist criticalparameter ranges that favor the buildup of boundary effects leading to strong enhancement of the(otherwise negligibly weak) nonlocality. However, the underlying mechanisms are fundamentallydifferent from those observed in the periodic case, and exhibit typical footprints (e.g., fractal gaps,quasi-localized states) that are distinctive of aperiodic order. The outcomes of our study indicatethat aperiodic order plays a key role in the buildup of the aforementioned boundary effects, andmay also find potential applications to optical sensors, absorbers and lasers. ∗ [email protected] a r X i v : . [ phy s i c s . op ti c s ] N ov . INTRODUCTION One key feature that distinguishes optical “metamaterials” [1–3] from other artificialmaterials such as photonic crystals [4] is the possibility to describe their macroscopic responsein terms of effective parameters (e.g., permittivity and permeability), along the lines of whatis conventionally done with natural materials. From the mathematical viewpoint, rigorousimplementations of this modeling process, typically referred to as “homogenization”, rely onfirst-principle concepts such as field averaging [5]. From the experimental viewpoint, sucheffective parameters can be retrieved via suitable measurements of the scattering matrix[6, 7].The basic, intuitive rationale underlying homogenization is that, as long as the electricalsizes of the material inclusions are very small on the wavelength scale, and their interactionsare weak , the fast field fluctuations inside the metamaterial are averaged out, and an electro-magnetic wave effectively “sees” a continuum whose constitutive properties are dictated bymixing formulae [8] which essentially depend on the inclusions’ material properties, shapes,orientations and proportions, but not on their sizes and spatial order. To give an examplethat is especially relevant for the present study, in a multilayered metamaterial composedby stacking two types of deeply subwavelength material layers (with distinct constitutiveproperties and thicknesses, labeled, e.g., with “ a ” and “ b ”), the effective parameters shoulddepend on the filling fractions (i.e., proportions of the a - and b -type constituents in the mix-ture) but not on the specific order and/or arrangement of the layers, so that configurationsassociated with sequences such as abababab , babababa and abbaabba should all be effectivelyequivalent, and should all differ from, e.g., aaabaaba [8].The inherent limitations and range of applicability of the simple “effective-medium the-ory” (EMT) above are well known, and more complex extensions have been developed tocapture the spatial-dispersion (nonlocal) effects [9, 10] which may become non-negligible,e.g., in the presence of electrically thick and/or metallic inclusions (see, e.g., Refs. [11–13]and [14–17] for general and multilayer-specific approaches, respectively). For instance, inmultilayered metamaterials, the presence of metallic layers (albeit deeply subwavelength)may induce strong nonlocal effects, due to the coherent interactions of surface-plasmon-polaritons [18] propagating at the metal-dielectric interfaces, which can manifest as theappearance of additional extraordinary waves [19] not predicted by the EMT.2uch less expectable and counterintuitive is the “breakdown” of the EMT in periodicmultilayered metamaterials with fully dielectric , deeply subwavelength layers, which was re-cently predicted on theoretical grounds by Herzig Sheinfux et al. [20], and experimentallyobserved by Zhukovsky et al. [21]. Basically, it was shown that, under specific illuminationsettings, the optical response (transmittance or reflectance) of finite-thickness samples mayexhibit substantial differences from the EMT prediction, accompanied by an ultrasensitivityto the spatial arrangement, size and termination of the layers. As also elucidated in follow-upstudies [22–26], these phenomena are not manifested in the bulk (infinite-medium) response,and can be interpreted as boundary effects stemming from the peculiar, interface-dominatedphase-accumulation mechanism in the multilayer, which may strongly enhance the (otherwisenegligibly weak) nonlocality. These effects can be captured by suitable nonlocal extensions[23, 26]. Related theoretical [27] and experimental [28] studies in similar parameter regimes,but characterized by random spatial disorder, have evidenced the possibility to attain An-derson localization, likewise in stark contrast with the EMT prediction, and once again withultrasensitivity to changes of features on a deeply subwavelength scale. These results havesparked considerable interest, both in terms of implications for the homogenization theory,and potential applications to extreme optical sensing and switching.Against the background above, this study explores the possibility to observe similar effectsin aperiodically ordered geometries, i.e., the vast middle ground separating perfect periodicityand random disorder. Originally inspired by the concept of “quasicrystals” in solid-statephysics [29, 30], aperiodic order has become increasingly relevant in many fields of science andtechnology [31] and, in particular, in optics and photonics [32] (see also a related perspectivein a recent roadmap on optical metamaterials [3, Sec. 3]). As a representative geometry, weconsider the Thue-Morse (ThM) sequence [33], which has been extensively studied in thepast in connection with photonic crystals [34–41] and metallo-dielectric multilayers [42], buthas never been explored in the fully dielectric, deeply subwavelength regime of interest here.Accordingly, the paper is structured as follows. In Sec. II, we introduce the problemgeometry and its formulation. In Sec. III, we describe the modeling tools and relatedmathematical formalism based on the trace and antitrace maps [43]. In Sec. IV, we discusssome representative results. Finally, in Sec. V, we draw some conclusions and point tofuture work. 3 I. PROBLEM STATEMENTA. Geometry
The problem geometry is schematized in Fig. 1. We consider a multilayer composed oftwo types of dielectric layers (labeled as “ a ” and “ b ”), with relative permittivity ε a and ε b ,and thickness d a and d b , stacked along the z -direction, and of infinite extent in the x − y plane. The layers are arranged aperiodically according to the ThM sequence, generated byassuming the symbol a as an initiator, and iteratively applying the substitution rules [33] a → ab, b → ba. (1)The first iterations are therefore a , ab , abba , abbabaab , and so on, with the generic n th stageof growth containing N = 2 n layers. As general, well-known traits of this sequence, werecall that at any iteration n ≥ i) the frequency of occurrence of a - and b -type symbols isidentical (and hence exactly the same as for periodic sequences), ii) each half of the sequencecorresponds to the “flipped” version of the other half, and iii) sequences containing morethan two consecutive symbols (e.g., aaa or bbb ) are not possible [33].In what follows, we consider a generic multilayer at stage of growth n , with total thickness D = 2 n − d (with d = d a + d b denoting the thickness of an ab -type bilayer), embedded in ahomogeneous dielectric medium with relative permittivity ε e . For instance, the case depictedin Fig. 1 corresponds to the stage of growth n = 4 (i.e., N = 16 layers). B. Formulation
For illumination, we assume a time-harmonic plane wave, with suppressed exp ( − iωt ) timedependence and transverse-electric (TE) polarization ( y -directed electric field), impingingwith an angle θ i from the z -axis, viz., E ( i ) y = E exp [ i ( k x x + k ze z )] , (2)where E denotes a real-valued amplitude, and k x = k √ ε e sin θ i , k ze = k √ ε e cos θ i (3)are the transverse (conserved) and longitudinal components, respectively, of the wavevector k e (see Fig. 1). In Eq. (3), k = ω/c is the vacuum wavenumber, with c denoting the4orresponding wavespeed. As previously mentioned, we assume to operate in the deeplysubwavelength regime, i.e., d a , d b (cid:28) λ , with λ = 2 π/k denoting the vacuum wavelength.Under these conditions, the optical response of the multilayer is generally well capturedby an EMT model in terms of a homogeneous, uniaxially anisotropic slab characterizedby a relative permittivity tensor whose parallel ( (cid:107) , i.e., x − y ) and orthogonal ( ⊥ , i.e., z )components are given by simple Maxwell-Garnett-type mixing formulae [8]¯ ε (cid:107) = f a ε a + f b ε b , (4a)¯ ε ⊥ = (cid:0) f a ε − a + f b ε − b (cid:1) − , (4b)with f a = d a /d and f b = d b /d = 1 − f a denoting the filling fractions pertaining to a - and b -type constituents, respectively, and the overbar utilized throughout the paper to indicateEMT-based quantities. First, we observe that the mixing formulae in Eqs. (4) are exactlyidentical with those pertaining to a conventional periodic multilayer (repetitions of ab -typebilayers). This should not be surprising, as we have previously recalled that, just like theperiodic ones, ThM sequences exhibit the same distribution of a - and b -type symbols, andthat EMT models are sensitive to proportions, rather than spatial arrangement. It is alsoworth pointing out that, in view of the assumed TE polarization, only the parallel componentin Eq. (4a) is actually relevant to our study.Previous studies on ThM-based optical structures have focused on photonic crystals (i.e.,moderately thick layers) [34–41] and hyperbolic metamaterials (i.e., deeply subwavelengthmetallic and dielectric layers) [42], which exhibit a wealth of interesting effects such asbandgaps, resonant transmission, localization and field enhancement, omnidirectional re-flection, fractal edge-states, multistability, additional extraordinary waves.Conversely, in what follows, we deal with ThM-based metamaterials featuring fully di-electric, deeply subwavelength layers, and study the possible buildup of boundary effectsleading to strong enhancement the (otherwise negligibly weak) nonlocality. To this aim, wesystematically compare the exact optical response of structures at various stages of growth,and under different illumination conditions, with the corresponding EMT-based predictions,in order to identify critical parameter regimes where nonlocality may be strongly enhanced.Moreover, to single out behaviors that are genuinely induced by the underlying aperiodicorder, we also consider the comparison with the well-established periodic-multilayer case[20–26] which, as observed above, shares the same EMT model.5 II. MODELING TOOLS AND FORMALISMA. Transfer-Matrix Model
The optical response of the ThM multilayered metamaterial in Fig. 1 can be rigorouslycalculated by means of the well-established transfer-matrix method [44, Chap. 1]. Basically,the transverse field components at the two interfaces of a generic a - or b -type layer can berelated via E ( L ) y iZ e H ( L ) x = M ν · E ( R ) y iZ e H ( R ) x , (5)where the superscripts ( L ) and ( R ) denote the left and right interfaces, respectively, Z e = ωµ k ze , (6)represents the TE wave impedance in the exterior medium (with µ denoting the vacuummagnetic permeability), and M ν = cos ( k zν d ν ) k ze k zν sin ( k zν d ν ) − k zν k ze sin ( k zν d ν ) cos ( k zν d ν ) , (7)is a 2 ×
2, unimodular, adimensional matrix, where ν = a or b , and k zν = (cid:112) k ε ν − k x = k (cid:112) ε ν − ε e sin θ i (8)denote the longitudinal wavenumbers in the two corresponding media [44, Chap. 1]. Therepresentation above can readily be iterated to deal with multiple cascaded layers, via chainproduct of the single-layer transfer matrices [44, Chap. 1]. Accordingly, we can relate thefields at the input ( z = 0) and output ( z = D ) interfaces of a ThM multilayer at stage ofgrowth n as E y iZ e H x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) z =0 = M ( n ) · E y iZ e H x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) z = D , (9a) M ( n ) = N =2 n (cid:89) j =1 M ν ( j ) = m ( n )11 m ( n )12 m ( n )21 m ( n )22 , (9b)with ν ( j ) = a or b , according to the jth -symbol in the ThM sequence.Based on the model in Eqs. (9), for a given incident field, we can rigorously calculate thereflection and transmission coefficients, as well as the field distribution inside the multilayer.6 . Trace and Antitrace Maps As already pointed out in our previous study dealing with the periodic case [26], somekey observables in the optical response of a generic multilayer can be calculated withoutthe need to actually perform the chain matrix product in Eq. (9b), which, for a largenumber of layers, may become both computationally intensive and prone to numerical-errorpropagation. For instance, by defining the transmission coefficient τ n = E ( t ) y (cid:12)(cid:12)(cid:12) z = D E ( i ) y (cid:12)(cid:12)(cid:12) z =0 , (10)with the superscript ( t ) tagging the transmitted field, we obtain from Eq. (9a) (see AppendixA for details) τ n = 2 m ( n )11 + m ( n )22 + i (cid:104) m ( n )21 − m ( n )12 (cid:105) = 2Tr (cid:104) M ( n ) (cid:105) + i ATr (cid:104) M ( n ) (cid:105) , (11)where Tr [ · ] and ATr [ · ] denote the conventional matrix trace and antitrace operators, respec-tively [45]. Quite remarkably, similar to the periodic multilayer case [26], also for the ThMgeometry of interest here it is possible to compute these quantities iteratively via simplepolynomial maps. More specifically, by letting χ n ≡ Tr (cid:104) M ( n ) (cid:105) , υ n ≡ Atr (cid:104) M ( n ) (cid:105) , ˜ υ n ≡ Atr (cid:104) ˜ M ( n ) (cid:105) , (12)with the tilde denoting a complementary configuration featuring a ThM sequence initiatedwith a b -type (instead of a -type) symbol, it can be shown [40, 42, 43] that the evolutionwith respect to the stage of growth n is ruled by the following intertwined maps χ n +2 = χ n ( χ n +1 −
2) + 2 , (13a) υ n +1 = χ n − [( χ n − υ n − + ˜ υ n − ] , (13b)˜ υ n +1 = χ n − [( χ n −
1) ˜ υ n − + υ n − ] , n ≥ , (13c)where χ = Tr (cid:16) M a (cid:17) , υ = Atr (cid:16) M a (cid:17) and ˜ υ = Atr (cid:16) M b (cid:17) .Though directly related to the evolution of the transmission coefficient [see Eq. (11)],trace and antitrace are not physically meaningful quantities, and hence cannot be used perse in order to ascertain the enhancement of nonlocality. Nevertheless, possible departures ofthe maps in Eqs. (13) from the corresponding EMT (local) predictions effectively quantify7he degree of nonlocality. Within this framework, for the periodic multilayer case [26], weshowed that the buildup of boundary effects leading to the enhancement of nonlocality couldbe effectively interpreted and parameterized in closed-form in terms of error propagation inthe trace and antitrace maps. In the ThM case of interest here, the trace and antitrace mapsin Eqs. (13) cannot be solved analytically in closed form. Nevertheless, the interpretationof the boundary effects in terms of error propagation still holds. It is worth stressing thatthe derivation of the trace and antitrace maps is exact , and therefore the computation ofthe transmission coefficient via Eqs. (11) and (13) is fully equivalent to that arising fromthe chain matrix product in Eqs. (9). In addition, the trace-antitrace-map scheme is alsocomputationally more effective and robust with respect to roundoff errors, as well as moreinsightful. IV. REPRESENTATIVE RESULTSA. Parameters and Observables
To facilitate comparison with previous studies on periodic and random structures, weconsider the same material parameters as in Refs. [20, 24, 26, 27], for the layers ( ε a = 1, ε b = 5, possibly with some small loss-gain perturbations) and exterior medium ( ε e = 4),with identical filling fractions f a = f b = 0 . d a = d b = d/ ε (cid:107) = 3. Likewise, we mainly focus on parameter configurationswhere the field is propagating in the higher-permittivity layers and in the effective medium,and evanescent in the lower-permittivity ones. This corresponds to an angular incidencerange θ ac ≡ arcsin (cid:18)(cid:114) ε a ε e (cid:19) < θ i (cid:46) arcsin (cid:18)(cid:114) ¯ ε (cid:107) ε e (cid:19) ≡ ¯ θ c . (14)As for the electrical thickness, we explore the range 0 . < d/λ < .
1, which guaranteesthat the layers remain deeply subwavelength. Our parametric studies below consider ThMmultilayers at various stages of growth n , which correspond to N = 2 n layers.For the lossless scenarios, besides the trace χ n and antitrace υ n , we consider as the mainphysical observables the transmittance T n = | τ n | = 4 | χ n + iυ n | , (15)8nd the electric field (magnitude) distribution in the multilayer [computed by means of thetransfer-matrix chain in Eqs. (9)]. For scenarios featuring optical losses or gain, we alsoconsider the reflectance R n = | ρ n | , (16)computed from the reflection coefficient (see Appendix A for details) ρ n = E ( r ) y E ( i ) y (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) z =0 = τ n (cid:104) m ( n )11 − im ( n )12 (cid:105) − , (17)with the superscript ( r ) tagging the reflected field. From Eq. (17), we observe that, unlikethe transmittance, the reflectance does not depend solely on the trace and antitrace. Westress that, in principle, it is possible to derive evolution maps [formally similar to those inEqs. (13)] for any of the transfer-matrix elements [43]. These, however, are not reportedhere for brevity.For lossy scenarios, we also compute the absorbance, which follows directly from powerconservation: A n = 1 − T n − R n . (18)To ascertain the possible enhancement of nonlocal effects, and the role played by aperiodicorder, we also study the two reference configurations considered in Ref. [26], namely, ahomogeneous slab with relative permittivity ¯ ε (cid:107) given by the EMT model in Eq. (4a) andthickness D , and a periodic multilayer with same type and total number of layers (and hencethickness D ). In both cases, the observables above can be computed analytically. For thehomogeneous EMT slab, they can be computed from a single transfer matrix as in Eq. (7)(by assuming ε ν = ¯ ε (cid:107) , d ν = D ), while for the periodic case they readily follow from theclosed-form solutions of the trace and antitrace maps [26, Eqs. (14)]. B. Lossless Case
We start considering the lossless scenario ( ε a = 1, ε b = 5). Figure 2 compares the trans-mittance responses of three representative ThM configurations at various stages of growth( n = 8 , ,
10) with the corresponding EMT and periodic benchmarks, as a function of theelectrical thickness d/λ of the ab -type bilayer (henceforth, simply referred to as “electricalthickness” for compactness) and the incidence direction θ i . This latter, according to Eq.914), varies within the range θ ac = 30 o < θ i (cid:46) ¯ θ c = 60 o . At a qualitative glance, we observea generally good agreement between the EMT [Figs. 2(d)–(f)] and periodic [Figs. 2(g)–(i)] responses, which exhibit the expectable small-to-moderate Fabry-P´erot-type oscillationsof the transmittance, with the possible exception of the region nearby the critical angle¯ θ c = 60 o , where the field undergoes a transition from propagating to evanescent in the effec-tive medium. Conversely, the ThM responses [Figs. 2(a)–(c)] exhibit a markedly differentbehavior also far away from the critical angle, with much more pronounced (bandgap-like)oscillations. In what follows, we examine in more detail two distinctive mechanisms under-lying these strong departures.
1. Near-Critical Incidence
In the periodic-multilayer case [20–26], significant differences between the exact opticalresponse and its EMT prediction were observed in the vicinity of the critical angle for whichthe field becomes evanescent in the effective medium ( θ i (cid:46) ¯ θ c ). In particular, we showed inRef. [26] that the trace and antitrace maps pertaining to the multilayer and a homogeneousEMT slab periodically depart according to a two-scale oscillatory law, whose maximumamplitudes may diverge asymptotically in the antitrace case (together with the slow scale)as the incidence direction approaches the critical angle ¯ θ c .For the ThM case of interest here, this regime remains critical, and other interestingeffects appear, which have no counterpart in the periodic scenario. Figure 3(a) shows arepresentative transmittance cut from Fig. 2(a) (ThM multilayer at stage of growth n = 8,i.e., N = 256 layers), at a fixed incidence angle θ i = 61 . o (cid:38) ¯ θ c , for which the field isevanescent in the effective medium. As it can be observed, the transmittance is very lowwithin most of the electrical-thickness range, but some high-transmittance resonant peaksappear for d/λ (cid:38) .
09. Conversely, the transmittance for the corresponding EMT andperiodic reference configurations remains always negligibly small ( < − ). Associated withthe high-transmittance peaks are some Fabry-P´erot-type states, as shown for in Fig. 3(b),which can exhibit strong field enhancements. Another representative example of such statesis shown in Fig. 3(c), for a higher stage of growth ( n = 12, i.e., N = 4096 layers).To gain some insight in this EMT-breakdown mechanism, which has no counterpartin the periodic case, it is instructive to look at the trace and antitrace maps. For fixed10lectrical thickness and incidence direction [corresponding to the Fabry-P´erot-type state inFig. 3(b)], Fig. 4 compares the evolution of the ThM trace, antitrace and transmittance,as a function of the stage of growth n , with those pertaining to the EMT and periodicconfigurations. We highlight that, for all three cases, trace and antitrace start from verysimilar values at the initial stage of growth, namely χ = 2 . υ = − .
550 for ThM andperiodic cases, and ¯ χ = 2 . υ = − .
482 in the EMT case, as an effect of the veryweak nonlocality. These values are only slightly beyond the “band-edge” condition χ = 2( k z = 0) for Bloch-type terminations [43], which indicates the occurrence of a bandgap. Asthe structure size increases, the maps pertaining to the EMT and periodic cases exhibitsimilar exponentially increasing behaviors, in stark contrast with those pertaining to theThM case, which remain bounded. Interestingly, the corresponding transmittance initiallyfollows the rapid, monotonic decrease of the EMT and periodic cases, but then abruptlyexhibits a “revival” at higher stages of growth ( n ≥ | χ | > | χ | = 2, thereby allowing a revival of the transmittance at higher stages ofgrowth. Although the ThM trace and antitrace maps in Eqs. (13) may actually exhibitperiodic orbits [43], this is not the case for the parameter configuration in Fig. 4, in spiteof the seeming periodicity (with n ) of ThM transmittance. For instance, no revivals areobserved for values 13 ≤ n ≤
20 (not shown).The above mechanism constitutes a first example of how negligibly weak nonlocality canbe enhanced by boundary effects so as to yield strong departures of the optical responsefrom the EMT prediction. We stress that, although these effects are still manifested inthe near-critical-incidence regime θ i ≈ ¯ θ c , they differ fundamentally from those observed inthe periodic-multilayer case [20–26], and can be genuinely attributed to the ThM aperiodicorder. 11 . Fractal Gaps and Quasi-Localized States Away from the near-critical incidence above, there are other distinctive features thatemerge in the ThM optical response. For a more quantitative assessment of the visualimpression from Figs. 2(a)–(c), Fig. 5(a) shows three representative cuts at fixed incidenceangle ( θ i = 35 . o ) and different stages of growth. For increasing size of the multilayer, weobserve the formation of gaps with growing complexity that resemble fractal-type structures.This is very different from the behavior of the EMT and periodic reference responses [cf. Fig.5(b)], which are in good agreement and exhibit only small oscillations around a near-unittransmittance.Fractal gaps are actually a well-known hallmark of ThM-based structures. Previousstudies on photonic crystals [37, 46] have explained the underlying mechanism in terms ofdistinctive interface correlation, and have derived the condition for a fractal gap to occur interms of a minimal bilayer electrical thickness √ ε a d a + √ ε b d b λ = 13 , (19)which, for our assumed parameters, corresponds to d/λ = 0 . d/λ ≈ .
06, i.e., by a factor ∼ . b -type layers and evanescent in the a -type ones. As a consequence,the phase accumulation is dominated by discrete jumps at the interfaces, rather than thepropagation across the layers [20].Figure 6 illustrates an interesting feature that is typically associated with fractal gapsin ThM optical structures, i.e., the appearance (at the lower or upper gap-edges) of stateswith hyperexponential localization properties lying somewhere in between the exponentialdecay of localized states and the extended character of Bloch-like gap-edge states in periodicstructures [37, 38, 46]. More specifically, Figs. 6(a) and 6(b) show two such states for theparameters as in Fig. 5, at two representative stages of growth, whereas Fig. 6(c) illustratesan example at a different incidence angle. Similar to what observed in the photonic-crystalregime [37, 38, 46], for increasing stages of growth, these “quasi-localized” states tend to12xhibit cluster-periodic distributions with large magnitude fluctuations and strong field en-hancement, and can attain very large quality factors, of potential interest for applicationsto optical cavities.To give an idea, the EMT prediction for the maximum field enhancement inside themultilayer can be expressed as (see Appendix B for details)¯ γ ≡ max 26) enhancement, in stark contrast with the actual values observed (4.3 and116, respectively).It is also instructive to look at the trace and antitrace maps. Figure 7 compares theevolutions of trace, antitrace and transmittance for the parameter configuration pertainingto the quasi-localized state in Fig. 6(b). As it can be observed, the EMT and periodicvalues maintain a generally good agreement, with only moderate oscillations in the traceand antitrace and near-unit transmittance, whereas the ThM ones exhibit markedly differentbehaviors for intermediate stages of growth 7 ≤ n ≤ 10. Qualitatively similar results can beobserved in connection with other quasi-localized states.The above results are a clear manifestation of a fundamentally different type of boundaryeffects, which can occur far away from the critical incidence, but are still genuinely inducedby the ThM aperiodic order.In what follows, with a view towards possible applications to sensing, absorbers andlasing, we study the effects of small losses and gain. C. Small Losses or Gain We add a small imaginary part to the permittivity of the b -type material, i.e., by assuming ε b = 5 + iδ , with | δ | (cid:28) 1; for the assumed time-harmonic convention, positive and negativevalues of δ correspond to optical losses and gain, respectively.Figure 8 shows some representative absorbance responses for the stage of growth n = 8( N = 256 layers), near-critical incidence ( θ i = 60 . o ), and δ = 10 − and 10 − (i.e., losses).13s it can be observed, there are a series of resonant peaks, even for electrical thicknessesas small as d/λ = 0 . δ ), since the field is evanescent inside the structure and gets almost completely reflected( ¯ R n ∼ . d/λ = 0 . δ = 10 − anddifferent stages of growth. Once again, several sharp peaks are observed, with absorbanceas high as 0.8. For three representative peaks, the comparison with the EMT and peri-odic counterparts is illustrated in Fig. 10 in terms of bar diagrams. In all three examples,the EMT and periodic look comparable, and substantially different from the ThM coun-terparts. For instance, the absorbance in the ThM case is significantly higher (by a factor20–50), and also the transmittance is quite different. As it can be observed in Fig. 11, thefield distributions at the resonant peaks exhibit the quasi-localized characteristics typical offractal-gap-edge states (cf. Fig. 6).We highlight that the results above pertain to specific parameter values, and in general theabsorbances exhibited by the three configurations (ThM, periodic, EMT) are comparable.Next, we consider a scenario featuring small optical gain, namely, δ = − − . Figure 12shows some representative transmittance responses in the vicinity of fractal gaps, at differentstages of growth, characterized by the presence of sharp peaks with very strong amplitudes(up to values of ∼ ), which are indicative of lasing conditions. Also in these cases, asshown in Fig. 13, the corresponding field distributions resemble quasi-localized states. Sim-ilar behaviors are also observed for the reflectance responses. Conversely, transmittance andreflectance for the EMT and periodic counterparts remain near-unit and very small, respec-tively. Qualitatively similar results (not shown for brevity) are also observed in connectionwith Fabry-P´erot-type resonant modes excited nearby the critical incidence.To better illustrate the difference between the observed response and the EMT prediction,we consider the lasing condition derived by enforcing a pole in the transmission (or reflection)coefficient for the EMT case (see Appendix B for details)tan (cid:0) ¯ k z D (cid:1) = 2 i ¯ k z k ze ¯ k z + k ze , (21)14hich can admit real-frequency solutions in the presence of gain. To give an idea, forthe parameters corresponding to the resonant peak in Fig. 12(b) ( n = 10, θ i = 31 . o , d/λ = 0 . D ∼ λ , i.e.,more than 60 time thicker than the ThM case. Alternatively, for the same overall thickness D = 31 . λ as for the ThM case, the EMT prediction yields a gain coefficient δ = − . D. Some Remarks A few remarks are in order on the assumptions and restrictions of our study. First, onemay argue that the material parameters considered in the multilayers (especially ε a = 1) arenot realistic for an experimental validation. As previously mentioned, the main motivationbehind this parameter choice was to facilitate direct comparison with the results in previousstudies on periodic and random scenarios [20, 24, 26, 27]. As also demonstrated by theexperimental studies in Refs. [21, 28], the phenomena of interest remain visible when realisticmaterials (e.g., silica and titania) are instead utilized. However, the lower the materialcontrast, the more difficult the observation. For the periodic multilayer case, in Ref. [26],we were able to derive analytically the relationship between the material contrast and thecritical size of the multilayer for which the breakdown phenomena could be observed. Forthe ThM case of interest here, an analytic study is not possible, but similar qualitativeconclusions are expected to hold. Within this framework, another potentially critical aspectis the exponential increase of the multilayer size with the stage of growth n . In the examplesshown, for the assumed parameters, the buildup effects leading to enhanced nonlocality turnout to occur for stages of growth n (cid:38) 8. While particularly high values of n would clearlylead to technologically unfeasible structures, we remark that stages of growth around thethreshold values n = 8 and n = 9 (i.e., hundreds of layers) are within reach for currentnanofabrication technologies, as demonstrated in recent experimental studies [28].Moreover, in connection with our assumption of TE polarization, once again to facilitate15irect comparison with the results from previous related studies [20, 24, 26, 27], we notethat enhanced nonlocality generally occurs for the transverse-magnetic polarization as well[20], although its visibility may be less pronounced [21].Finally, we remark that, unlike the periodic-multilayer case [26], it was not possible hereto identify some closed-form parameters relating the optical response with the ThM aperiodicgeometry, due to the more complex character of the arising trace and antitrace maps (notsolvable analytically). Nevertheless, our parametric studies elucidate some representativemechanisms and effects that are distinctive of the ThM geometry. V. CONCLUSIONS AND PERSPECTIVES To sum up, we have shown that, in aperiodically ordered, fully dielectric multilayeredmetamaterials based on the ThM geometry, the inherently weak nonlocality exhibited in thedeeply subwavelength regime can be substantially enhanced via the buildup of boundaryeffects that are fundamentally different from those observed in the periodic case [20–26].These effects are manifested as strong departures of the optical response (reflectance andtransmittance, as well as absorbance or lasing in the presence of small loss or gain, respec-tively) from the EMT prediction and periodic counterpart, with distinctive footprints suchas fractal gaps and quasi-localized states.We stress, once again, that the comparison with the periodic case is particularly mean-ingful, since the two geometries (ThM and periodic) contain exactly the same amounts ofeach of the material constituents, the only difference being the spatial order. This providesfurther evidence that, even at deeply subwavelength scales, spatial order may strongly affectthe optical response.Our outcomes constitute a first step toward extending the previous studies on periodic[20–26] and randomly disordered [27, 28] geometries to the intermediate realm of “orderlydisorder” and, albeit focused on a specific geometry, provide some generally applicable tools.As shown in Refs. [47, 48], the trace and antitrace map formalism can in principle be appliedto generic aperiodic sequences based on two-symbol substitution rules. Therefore, amongthe possible follow-up studies, it looks very intriguing to explore different aperiodicallyordered geometries. For instance, it would be very interesting to explore to what extentsome distinctive properties of the optical response of Fibonacci-type photonic quasycrystals16e.g., self-similarity in the spectrum, critical states of multifractal nature, etc.) [49] are alsoobservable in the deeply subwavelength regime. Also of great interest are more application-oriented studies on the promising potentials that have emerged in connection with opticalsensors, absorbers and lasers. Within this framework, we are currently pursuing a systematicstudy of the effects of enhanced nonlocality on the (bulk and surface) optical sensitivityresponse of ThM multilayers, which will be the subject of a forthcoming paper. Appendix A: Details on Eqs. (11) and (17) In view of Eq. (2) and the definitions of the transmission and reflection coefficients inEqs. (10) and (17), respectively, the total electric fields at the input ( z = 0) and output( z = D ) interfaces can be written as E y ( x, z = 0) = E (1 + ρ n ) exp ( ik x x ) , (A1a) E y ( x, z = D ) = E τ n exp ( ik x x ) . (A1b)By calculating the corresponding tangential magnetic fields from the relevant Maxwell’sequation, the matrix equation in Eq. (9a) can be rewritten as ρ n − i (1 − ρ n ) = M ( n ) · τ n − iτ n , (A2)from which Eqs. (11) and (17) readily follow by solving the linear system of equations, andexploiting the unimodular character of the matrix. Appendix B: Details on Eqs. (20) and (21) To calculate the EMT predictions, we consider a homogeneous slab of thickness D andrelative permittivity ¯ ε (cid:107) . For given incidence conditions, the electric field inside the structurehas the form of a standing wave¯ E y ( z ) = ¯ E + exp (cid:2) i ¯ k z ( z − D ) (cid:3) (cid:8) (cid:2) − i ¯ k z ( z − D ) (cid:3)(cid:9) , (B1)where an irrelevant exp ( ik x x ) term is omitted, ¯ E + is a complex amplitude to be determined,¯Γ = Z e − ¯ ZZ e + ¯ Z (B2)17s the partial reflection coefficient between the exterior and effective media, and¯ Z = ωµ ¯ k z (B3)is the TE wave impedance of the effective medium. In the exterior region z < 0, the totalfield is obtained by summing the incident [see Eq. (2)] and reflected contributions, viz.,¯ E y ( z ) = E [1 + ¯ ρ exp ( − ik ze z )] , (B4)where ¯ ρ can be obtained from Eq. (17) by assuming a single layer of thickness D and relativepermittivity ¯ ε (cid:107) . By enforcing the continuity of the electric field at the interface z = 0, aftersome algebra, we obtain¯ E + = 2 E ¯ Z (cid:0) Z e + ¯ Z (cid:1) exp (cid:0) − i ¯ k z D (cid:1) exp (cid:0) − i ¯ k z D (cid:1) (cid:0) Z e + ¯ Z (cid:1) − (cid:0) Z e − ¯ Z (cid:1) . (B5)By recalling that, in view of the assumed parameters, ε e > ¯ ε (cid:107) and hence Z e < ¯ Z , we observefrom Eq. (B2) that ¯Γ < 0. Therefore, by assuming the slab thicker than half a wavelength,it follows from Eq. (B1) that¯ γ = (cid:12)(cid:12) ¯ E + (cid:12)(cid:12) (cid:0) − ¯Γ (cid:1) E = 4 ¯ Z (cid:0) Z e + ¯ Z (cid:1) − (cid:0) Z e − ¯ Z (cid:1) = ¯ ZZ e = k ze ¯ k z , (B6)which corresponds to the result in Eq. (20).Likewise, from Eq. (11), the EMT prediction of the transmission coefficient is¯ τ = 1cos (cid:0) ¯ k z D (cid:1) − i (cid:18) ¯ k z k ze + k ze ¯ k z (cid:19) sin (cid:0) ¯ k z D (cid:1) , (B7)from which the lasing condition in Eq. (21) directly follows by zeroing the denominator. [1] F. Capolino, Theory and Phenomena of Metamaterials (CRC Press, Boca Raton, FL, 2009).[2] W. Cai and V. M. Shalaev, Optical Metamaterials: Fundamentals and Applications (Springer,New York, 2010).[3] A. M. Urbas, Z. Jacob, L. Dal Negro, N. Engheta, A. D. Boardman, P. Egan, A. B. Khanikaev,V. Menon, M. Ferrera, N. Kinsey, C. DeVault, J. Kim, V. Shalaev, A. Boltasseva, J. Valentine,C. Pfeiffer, A. Grbic, E. Narimanov, L. Zhu, S. Fan, A. Al`u, E. Poutrina, N. M. Litchinitser, . A. Noginov, K. F. MacDonald, E. Plum, X. Liu, P. F. Nealey, C. R. Kagan, C. B. Murray,D. A. Pawlak, I. I. Smolyaninov, V. N. Smolyaninova, and D. Chanda, J. Opt. , 093005(2016).[4] J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Moldingthe Flow of Light , 2nd ed. (Princeton University Press, Princeton, NJ, USA, 2008).[5] D. R. Smith and J. B. Pendry, J. Opt. Soc. Am. B , 391 (2006).[6] D. R. Smith, D. C. Vier, T. Koschny, and C. M. Soukoulis, Phys. Rev. E , 036617 (2005).[7] S. Arslanagi´c, T. V. Hansen, N. A. Mortensen, A. H. Gregersen, O. Sigmund, R. W. Ziolkowski,and O. Breinbjerg, IEEE Antennas Propagation Mag. , 91 (2013).[8] A. Sihvola, Electromagnetic Mixing Formulas and Applications , Electromagnetics and RadarSeries (IET, London, UK, 1999).[9] D. Landau and E. M. Lifshitz, Electrodynamics of Continuous Media (Pergamon Press, NewYork, 1960).[10] V. Agranovich and V. Ginzburg, Crystal Optics with Spatial Dispersion, and Excitons ,Springer Series in Solid-State Sciences (Springer-Vergal, Berlin/Heidelberg, 2013).[11] M. G. Silveirinha, Phys. Rev. B , 115104 (2007).[12] A. Al`u, Phys. Rev. B , 075153 (2011).[13] A. Ciattoni and C. Rizza, Phys. Rev. B , 184207 (2015).[14] J. Elser, V. A. Podolskiy, I. Salakhutdinov, and I. Avrutsky, Appl. Phys. Lett. , 191109(2007).[15] A. V. Chebykin, A. A. Orlov, A. V. Vozianova, S. I. Maslovski, Y. S. Kivshar, and P. A.Belov, Phys. Rev. B , 115438 (2011).[16] A. V. Chebykin, A. A. Orlov, C. R. Simovski, Y. S. Kivshar, and P. A. Belov, Phys. Rev. B , 115420 (2012).[17] R.-L. Chern, Opt. Express , 16514 (2013).[18] S. Maier, Plasmonics: Fundamentals and Applications (Springer-Verlag, New York, NY, USA,2007).[19] A. A. Orlov, P. M. Voroshilov, P. A. Belov, and Y. S. Kivshar, Phys. Rev. B , 045424(2011).[20] H. Herzig Sheinfux, I. Kaminer, Y. Plotnik, G. Bartal, and M. Segev, Phys. Rev. Lett. ,243901 (2014). 21] S. V. Zhukovsky, A. Andryieuski, O. Takayama, E. Shkondin, R. Malureanu, F. Jensen, andA. V. Lavrinenko, Phys. Rev. Lett. , 177402 (2015).[22] A. Andryieuski, A. V. Lavrinenko, and S. V. Zhukovsky, Nanotechnology , 184001 (2015).[23] V. Popov, A. V. Lavrinenko, and A. Novitsky, Phys. Rev. B , 085428 (2016).[24] X. Lei, L. Mao, Y. Lu, and P. Wang, Phys. Rev. B , 035439 (2017).[25] A. Maurel and J.-J. Marigo, Phys. Rev. B , 024306 (2018).[26] G. Castaldi, A. Al`u, and V. Galdi, Phys. Rev. Appl. , 034060 (2018).[27] H. Herzig Sheinfux, I. Kaminer, A. Z. Genack, and M. Segev, Nat. Commun. , 12927 (2016).[28] H. Herzig Sheinfux, Y. Lumer, G. Ankonina, A. Z. Genack, G. Bartal, and M. Segev, Science , 953 (2017).[29] D. Shechtman, I. Blech, D. Gratias, and J. W. Cahn, Phys. Rev. Lett. , 1951 (1984).[30] D. Levine and P. J. Steinhardt, Phys. Rev. Lett. , 2477 (1984).[31] E. Maci´a, Rep. Progr. Phys. , 397 (2006).[32] L. Dal Negro and S. Boriskina, Laser Photonics Rev. , 178 (2011).[33] M. Queff´elec, Substitution Dynamical Systems - Spectral Analysis , Lecture Notes in Mathe-matics (Springer-Verlag, Berlin/Heidelberg, 2010).[34] N.-H. Liu, Phys. Rev. B , 3543 (1997).[35] F. Qiu, R. W. Peng, X. Q. Huang, Y. M. Liu, M. Wang, A. Hu, and S. S. Jiang, Europhys.Lett. , 853 (2003).[36] L. Dal Negro, M. Stolfi, Y. Yi, J. Michel, X. Duan, L. C. Kimerling, J. LeBlanc, and J. Haav-isto, Appl. Phys. Lett. , 5186 (2004).[37] X. Jiang, Y. Zhang, S. Feng, K. C. Huang, Y. Yi, and J. D. Joannopoulos, Appl. Phys. Lett. , 201110 (2005).[38] J. M. Marianne Hiltunen and L. Dal Negro, in Proc. SPIE , Vol. 6989 (2008) pp. 6989–7.[39] V. Grigoriev and F. Biancalana, New J. Phys. , 053041 (2010).[40] V. V. Grigoriev and F. Biancalana, Photonics Nanostruct. Fundam. Appl. , 285 (2010).[41] W. J. Hsueh, S. J. Wun, Z. J. Lin, and Y. H. Cheng, J. Opt. Soc. Am. B , 2584 (2011).[42] S. Savoia, G. Castaldi, and V. Galdi, Phys. Rev. B , 235116 (2013).[43] X. Wang, U. Grimm, and M. Schreiber, Phys. Rev. B , 14020 (2000).[44] M. Born and E. Wolf, Principles of Optics , 7th ed. (Cambridge University Press, 1999).[45] S. Lang, Linear Algebra , 3rd ed. (Springer, Berlin, 1987). 46] H. Lei, J. Chen, G. Nouet, S. Feng, Q. Gong, and X. Jiang, Phys. Rev. B , 205109 (2007).[47] M. Kol´aˇr and M. K. Ali, Phys. Rev. A , 7112 (1990).[48] M. Kol´aˇr and F. Nori, Phys. Rev. B , 1062 (1990).[49] Z. V. Vardeny, A. Nahata, and A. Agrawal, Nat. Photon. , 177 (2013). e yzx EH k e d a d ba b a a a a a a a ab b b b b b b b i(i)(i) D FIG. 1. Problem geometry illustrating the ThM multilayered metamaterial and illumination con-ditions (details in the text). r an s m i tt an c e θ i ( deg ) n=8 Th M (a) n=9(b) n=10(c)(d) E M T θ i ( deg ) (e) (f) P e r i od i c (g) θ i ( deg ) d/λ (h) d/λ (i) d/λ FIG. 2. (a), (b), (c) Transmittance [see Eq. (15)] responses pertaining to ThM multilayeredmetamaterials at stages of growth n = 8 ( N = 256 layers), n = 9 ( N = 512 layers), and n = 10( N = 1024 layers), respectively, for ε a = 1, ε b = 5, f a = f b = 0 . d a = d b = d/ ε e = 4,as a function of the electrical thickness d/λ and incidence angle θ i . (d), (e), (f) Same as above,but EMT predictions (¯ ε (cid:107) = 3). (g), (h), (i) Same as above, but for periodic arrangements. r an s m i tt an c e −4 −3 d/λ (a) N o r m a li z ed m agn i t ude Layer No. (b) N o r m a li z ed m agn i t ude Layer No. (c) FIG. 3. Parameters as in Fig. 2. (a) Transmittance-cut from Fig. 2(a) [ThM multilayer at stageof growth n = 8 ( N = 256 layers)], at θ i = 61 . o . The corresponding EMT and periodic referenceresponse (not shown), are negligibly small ( < − ). (b) Field (magnitude) distribution inside themultilayer (normalized by the incident-field amplitude E ) at d/λ = 0 . n = 12 ( N = 4096 layers). hMEMTPeriodic No. of layers, N=2 n T r a c e −4−20246810 (a) A n t i t r a c e −30−20−10010 (b) T r an s m i tt an c e Stage of growth, n (c) FIG. 4. Parameters as in Fig. 2, with d/λ = 0 . 092 and θ i = 61 . o . (a), (b), (c) Comparisonsbetween the trace (blue circles), antitrace (green triangles) and transmittance (red squares) evolu-tions pertaining to ThM, EMT and periodic configurations, respectively, as a function of the stageof growth n . The corresponding number of layers is also shown on the top axis. Continuous curvesare guides to the eye only. =8n=9n=10 T r an s m i tt an c e −9 −6 −3 (a) EMTPeriodic T r an s m i tt an c e d/λ (b) FIG. 5. Parameters as in Fig. 2. (a) Transmittance-cuts at θ i = 35 . o from Fig. 2(a), forThM multilayers at stages of growth n = 8 ( N = 256 layers; blue-solid), n = 9 ( N = 512 layers;red-dashed), and n = 10 ( N = 1024 layers; green-dotted). Note the logarithmic scale and verylarge dynamics on the vertical axis. (b), (c) Same as panel (a) but for the corresponding EMT(purple-dashed) and periodic configurations (magenta-solid) at stage of growth n = 8 [i.e., cutsfrom Figs. 2(b) and 2(c), respectively]. Note the linear scale and much smaller dynamics on thevertical axis. o r m a li z ed m agn i t ude Layer No. (a) N o r m a li z ed m agn i t ude Layer No. (b) N o r m a li z ed m agn i t ude Layer No. (c) FIG. 6. Parameters as in Fig. 2. Field (magnitude) distributions inside the ThM multilayer(normalized by the incident-field amplitude E ) for representative quasi-localized states. (a) d/λ =0 . θ i = 35 . o , n = 8 ( N = 256). (b) d/λ = 0 . θ i = 35 . o , n = 12 ( N = 4096). (c) d/λ = 0 . θ i = 49 . o , n = 12 ( N = 4096). hMEMTPeriodic No. of layers, N=2 n T r a c e −100−80−60−40−200 (a) A n t i t r a c e −25−20−15−10−505 (b) T r an s m i tt an c e Stage of growth, n (c) FIG. 7. Parameters as in Fig. 6(b). (a), (b), (c) Comparisons between the trace (blue circles),antitrace (green triangles) and transmittance (red squares) evolutions pertaining to ThM, EMTand periodic configurations, respectively, as a function of the stage of growth n . The correspondingnumber of layers is also shown on the top axis. Continuous curves are guides to the eye only. b s o r ban c e d/λ δ =10 -4 δ =10 -3 N o r m a li z ed m agn i t ude Layer No. FIG. 8. Parameters as in Fig. 2, but with ε b = 5 + iδ (with δ > 0, i.e., losses), n = 8 (i.e., N = 256layers), and near-critical incidence θ i = 60 . o . Absorbance as a function of the electrical thickness,for δ = 10 − (blue-solid) and δ = 10 − (red-dashed). The inset shows the normalized (magnitude)field distributions of the Fabry-P´erot states corresponding to the peaks at d/λ = 0 . b s o r ban c e (a) A b s o r ban c e (b) A b s o r ban c e d/λ (c) FIG. 9. Parameters as in Fig. 2, but with ε b = 5 + i − (losses). Absorbance [from Eq. (18)] asa function of the ab -type bilayer electrical thickness, for (a) n = 9 ( N = 512 layers), θ i = 50 . o ,(b) n = 10 ( N = 1024 layers), θ i = 48 . o , (c) n = 12 ( N = 4096 layers), θ i = 40 . o . bsorbance Transmittance Reflectance T h M P e r i od i c E M T T h M P e r i od i c E M T T h M P e r i od i c E M T (c)(b)(a) FIG. 10. (a), (b), (c) Bar diagrams illustrating the absorbance, transmittance and reflectance atthree representative peaks ( d/λ = 0 . , . , . o r m a li z ed m agn i t ude Layer No. (a) N o r m a li z ed m agn i t ude Layer No. (b) N o r m a li z ed m agn i t ude Layer No. (c) FIG. 11. (a), (b), (c) Normalized (magnitude) field distributions of the quasi-localized statescorresponding to the absorbance peaks in Fig. 10. r an s m i tt an c e (a) T r an s m i tt an c e (b) T r an s m i tt an c e d/λ (c) FIG. 12. Parameters as in Fig. 2, but with ε b = 5 − i − (gain). Transmittance as a functionof the ab -type bilayer electrical thickness for (a) n = 9 ( N = 512 layers), θ i = 46 . o , (b) n = 10( N = 1024 layers), θ i = 31 . o , (c) n = 12 ( N = 4096 layers), θ i = 56 . o . o r m a li z ed m agn i t ude Layer No. (a) N o r m a li z ed m agn i t ude Layer No. (b) N o r m a li z ed m agn i t ude Layer No. (c) FIG. 13. (a), (b), (c) Normalized (magnitude) field distributions of the quasi-localized statesassociated with representative transmittance peaks ( d/λ = 0 . , . , .085, respectively) in thecorresponding panels of Fig. 12.