1. Introduction

Moiré effects encompass a variety of physical phenomena resulting from superimposing two periodic structures with a slight geometric mismatch. A well-known setup exhibiting these effects are moiré metasurfaces, parallel stacks of metasurfaces with a relative twist angle [1]. The study of moiré effects has attracted significant attention in recent research, as it promises novel avenues to tune the properties of light and its interaction with matter [2–4]. This has already led to profound findings such as light localization in lattice structures written in photorefractive materials [5], the ability to engineer the plasmon dispersion relation in evanescently coupled metasurfaces [2] and the discovery of exotic optical states exhibiting characteristics of both moiré flat bands and (quasi-)bound states in the continuum [6]. However, a deeper understanding of moiré effects from first principles and their impact on photonic band structures remains a prime challenge for current research in nano-optics [1,7,8] and numerical and computational methods that can explore such systems are investigated more recently [9–14]. Moiré flat bands have also been explored in graphene nanoribbons [15], electric circuits with extremely high degrees of freedom [16] and 1D electronic geometries [17]. Since the understanding of moiré flat bands in bilayer photonic crystals has been deepened considerably by analytical means [18], we aim to provide a comparable treatment of 1D photonic structures composed of atomic dipoles. The lack of such a treatment is hindering a more thorough understanding of higher-dimensional geometries, built from e.g. individual meta-atoms, as well.

To respond to this challenge, this work explores one-dimensional linear arrays of atoms characterized by an electric dipolar response [19,20]. By stacking two of these arrays on top of each other, moiré effects can be investigated by introducing a relative lattice constant mismatch between the upper and the lower array, parameterized by a moiré parameter $\theta$ (see Fig. 1) [21,22]. In the limit of strong coupling, where the spatial separation between the two arrays is minimal, heavily localized collective excitations emerge that are approximately characterized by infinite lifetimes and a vanishing dispersion. In the band structure, these excitations appear as flat bands of a near-lossless system, while they manifest as sharp peaks in the density of optical states (DOS). The spectral location of these field-enhancing flat bands can be revealed to be a moiré effect that is controllable mainly via adjustment of the moiré parameter $\theta$. Using an analytical approach relying on the diagonalization of the system’s interaction matrix, the spectral location of these bands can be calculated from basic expressions depending on the distances between nearest-neighbor pairs in a unit cell. A comparison to numerical simulation results obtained for finite systems demonstrates that spectral control of system resonances is largely achievable by an appropriate choice of the moiré parameter $\theta$. In our treatment, we focus on rational values of $\theta$ since incommensurable configurations, corresponding to irrational values of $\theta$, do not exhibit a periodically repeated unit cell. This makes them inaccessible to our numerical method and precludes us from applying Bloch’s theorem, which our theoretical treatment crucially relies on. However, since any irrational value of $\theta$ might be approximated by a sufficiently close rational one, this presents no inherent limitation to the applicability of our results.

 

Fig. 1. Two moiré unit cells for a setup of two stacked arrays with lattice constants $\Lambda _2$, $\Lambda _1$. The lateral separation between the two arrays is given by $d$. The depicted geometry is described by a moiré parameter value of $\theta = \frac {3}{4}$. The periodicity of the overall setup is then given by $\Lambda = 3 \Lambda _1 = 4 \Lambda _2$. The $i$-th atom in the $n$-th unit cell carries a dipole moment $\textbf{p}_{i,n}$. By Bloch’s theorem, these dipole moments are identical for equivalent atoms in different unit cells up to a phase. This relationship is indicated by an identical color code.

Download Full Size | PDF

Our results provide a better understanding of photonic moiré effects and point to their potential use in future photonic devices such as optical sensors and light traps. The remainder of this paper is organized as follows: In Sec. 2, we briefly introduce the theoretical framework underlying the description of dipolar collective excitations in atomic structures of finite spatial extent. A particular focus is given to the finite system’s interaction matrix, a non-Hermitian Hamiltonian capturing atomic dipole-dipole interactions, which our numerical treatment is based on. This quantity delivers the spectral information needed to construct approximate band structures and compute field enhancement in terms of the DOS. In Sec. 3, we introduce the interaction matrix for infinite, periodic structures. The resulting matrix is a function of crystal momenta in the irreducible Brillouin zone and captures the effective lattice-mediated interaction between dipoles in a unit cell. It is intimately related to the appearance of effective polarizabilities, detailed, e.g., in Refs. [23] or [24]. This $k$-space interaction matrix enables us to derive analytical expressions for the band structure and decay rates of periodic arrays. These expressions take comparatively accessible form in the perturbative limit of small spatial inter-array separation. We identify this limit as an inter-array nearest-neighbor approximation that still retains the total intra-array lattice interactions. This lattice-mediated NN-approximation provides valuable insight into the underlying physics of the phenomena under study and thus guides our numerical investigations. In Sec. 4, we develop the perturbative approach described above starting from perfectly aligned geometries and finally apply it in Sec. 5 to study the effects of moiré geometries. It is revealed that moiré flat bands arise from collective excitations characterized by out-of-phase atomic dipole oscillations in different arrays. The spectral position of a flat band is given by a comparatively simple composition of harmonic functions and can be controlled by adjusting inter-array nearest-neighbor distances via appropriate choices for $\theta$. We then conclude our investigation with a summary of the results in Sec. 6. A brief derivation of the interaction matrix in the finite and infinite periodic case is presented in Appendix A. Details on the numerical simulation procedure and links to the source code are provided in Appendix B. In Appendix C, we illustrate the validity of the presented results in the strong-coupling regime, while in Appendix D, we discuss the impact of fabrication errors conceptualized as positional disorder on moiré flat bands.

2. Interaction matrices for finite systems

We consider a collection of identical atoms with an electric dipolar response and model the polarizability of a single atom by a Lorentzian [25,26]

The resonance behaviour of this geometry is captured by a $3N \times 3N$ non-Hermitian matrix $M$. It is called the interaction matrix and plays the role of an effective, non-Hermitian Hamiltonian describing the atomic dipole-dipole interactions. Denoting the position of the $i$-th atom by $\textbf{r}_i$, it can be expressed in block-form as follows

Here, we denote the $n$-dimensional unit matrix by $I_n$. The matrix $G(k, \textbf{r})$ is the dyadic Green’s tensor, whose components are given by [29,30]

In Eqn. (4), we have defined the resonance frequency of the quasi-mode relative to the single atom resonance frequency $\omega _0$ for convenience. The connection of these quantities to the field enhancement is given by the normalized density of optical states (DOS) of atomic excitations at frequency $\omega$, which is given by [31]

In this paper, we numerically diagonalize the interaction matrix for finite 1D geometries. We aim to identify co-occurring system resonances, i.e., a set of resonance indices $L$ such that for each index combination $l,l' \in L$, it holds that $\Delta _l \approx \Delta _{l'}$. From Eqn. (6), it can be seen that for $\omega - \omega _0$ close to these co-occurring resonances, the DOS strongly increases. This implies large field enhancement in the corresponding spectral region.

To gain a deeper understanding of the underlying physical effects, we additionally derive closed analytical expressions for quantities analogous to Eqns. (4) and (5) in infinite systems. Since the infinite case deals with system properties in reciprocal space, where the resonance and line width depend on a 1D crystal momentum $k$, we also need to construct approximate dispersion diagrams for finite systems.

This can be achieved by associating a crystal momentum to each resonance and decay rate in Eqns. (4) and (5) by a Fourier analysis of the corresponding eigenvector. This approach is outlined in greater detail in Appendix B. Represented as part of the band structure, co-occurring resonances manifest as flat bands in the dispersion. We then study how moiré tuning of the external geometry can influence the occurrence and spectral location of these flat bands.

3. Interaction matrices for infinite periodic systems

An infinite system of resonant dipoles with periodicity $\Lambda$ and $N$ dipoles in the unit cell can again be described by a $3N \times 3N$ interaction matrix $M_P(k)$, which takes the following block form

Following Eqns. (4) and (5), we obtain the band structure and collective decay rates of the infinite system from the $3N$ eigenvalues $m_{P,l}(k)$ of $M_P(k)$ as

Additionally, it is useful to express both the collective resonance shift and the decay rates in terms of the single-atom decay rate as follows

We do so for quantities relating to finite and periodic systems. However, in the case of periodic systems, we still indicate the dependence of the normalized quantities above on the crystal momentum with the notation $D(k), G(k)$.

4. Perfectly aligned geometries

Our initial focus is on perfectly aligned geometries without moiré effects to develop a perturbative ansatz concerning the exact solution for a single linear dipolar array.

Let us first consider a single infinite periodic array with lattice constant $\Lambda$, described by a single dipole in a unit cell. Such a setup has already been analyzed in the literature by analytical means slightly different from ours, e.g. in [32], but we briefly discuss it here as it serves as the basis of our perturbation theory and a convenient illustration of the overall methodology. We chose the direction of periodicity to be along the $z$-axis and compute the corresponding $3 \times 3$ interaction matrix $M_{0}(k)$ from Eqn. (7)

The eigenvalues of $M_0(k)$ are identical to its diagonal elements. To arrive at closed expressions, we split the series occurring in Eqns. (13) and (14) into real and imaginary parts, which yield well-known limits, which are given in Ref. [33]. The real part converges to a sum of the first three generalized Clausen functions $\operatorname {Cl}_{i}$. In contrast, the imaginary part converges to a sum of Bernoulli polynomials, which are restricted to be periodic on the interval $[0, 2\pi ]$. Collecting terms, and with the dimensionless wave numbers

We note that formally equivalent expressions have been derived elsewhere, e.g., in Refs. [34] and [32]. Figure 2 displays a comparison to approximate band structures for a finite system of 40 atoms. Of special interest in the above equations is the fact that modes with a crystal momentum $k$ greater than the single atom resonance wave number $k_0$, are perfectly guided and do not radiate to infinity as their decay rates vanish.

 

Fig. 2. Normalized band structure (left) and collective decay rates (right) for a single array with lattice constant $\Lambda = 0.2 \lambda _0$. The solid lines show the analytical expression for the infinite case, while the dots mark the simulation results for the approximate band structure of a finite system of 40 atoms. The dashed gray line indicates the light line, where the crystal momentum is identical to the atomic resonance wave number $k_0$.

Download Full Size | PDF

We use Eqn. (12) as the starting point of our perturbative treatment by considering two identical, perfectly aligned arrays separated by a distance $d$ in $y$-direction (this corresponds to $\Lambda _1 = \Lambda _2$ in Fig. 1). Since the unit cell of this geometry encompasses two atoms, the interaction matrix $M(k)$ is $6 \times 6$ and can again be computed from Eqn. (7)

We consider the case of strong geometric coupling $d << \lambda _0$ and take the limit $d \rightarrow 0$, effectively collapsing the two arrays into one. To do this, we isolate the singular contributions in the lattice summations that always originate from nearest neighbor interactions with $n=0$ and retain them as explicit functions of $d$ while dropping all other powers of $d$. Physically, this makes our treatment equivalent to approximating the total inter-array interaction by its nearest-neighbor contribution, while still taking the full intra-array interactions into account. Inspecting the occurring sums, one can see that the off-diagonal elements in Eqn. (26) vanish in this limit because the summand $n=0$ does not contribute, and the remaining terms decrease linearly with $d$. Similarly, the first term in Eqn. (24) may be neglected, leaving us with two distinct singularities in Eqns. (23)–(25), which need to be retained. Taking this together, the inter-array coupling matrix takes the following diagonal form, valid in the perturbative limit

Thus, the spectrum separates into a $k$-independent and a $k$-dependent part, a phenomenon that generalizes to moiré setups. We stress that the limit of strong spatial coupling $d \rightarrow 0$, which lies at the heart of our perturbative approach, is necessary for the interaction matrix to obtain the easily diagonalizable form above. From the real and imaginary part of these eigenvalues, the associated band structure and decay rates can again be computed using Eqns. (10) and (9). To this end, we define

 

Fig. 3. Approximate band structure (left) and associated DOS (right) for a system of perfectly aligned arrays with $\Lambda = 0.3 \lambda _0$, separated by a distance $d = 0.1 \lambda _0$ containing 60 atoms. The colorbar marks the associated decay rates. The solid lines indicate the band structure according to the nearest-neighbor approximation given by Eqns. (36)–(38). As can be seen, flat bands are nearly lossless and correspond to a strongly enhanced density of optical states. The dashed grey line in the left plot marks the light line given by $k=k_0$ as in Fig. 2. We note that, in accordance with Eqns. (41)–(43), the flat bands exhibit strongly suppressed losses. Additionally, we remark on the degeneracy of two of the flat bands $D_{1,-}(k)$ and $D_{3,-}(k)$ obtained from Eqns. (36) and (38).

Download Full Size | PDF

To obtain the decay rates, we substitute $x = d k_0$ and compute

The singularity of these expressions in the limit $x \rightarrow 0$ is removable. Application of L’Hôpital’s rule then yields $\lim _{x \rightarrow 0} \text {Im}(c_1(x)) = \lim _{x \rightarrow 0} \text {Im}(c_2(x)) = \frac {k_0}{6 \pi } = -(\alpha _0 k_0^2)^{-1}$. Together with Eqns. (32) and (33), this ultimately gives

Of note is the appearance of lossless flat bands, encoded by $F_1(d), F_2(d)$, which are compositions of harmonic functions. We note that the flat bands obtained from Eqns. (36) and (38) are degenerate. In Fig. 3, the connection of the band structure to the DOS of a finite system with 60 atoms is illustrated. The physical meaning behind the occuring flat bands is the out-of-phase oscillation of atomic dipole moments in different arrays resulting in highly suppressed radiation. A mathematical justification for this explanation will be given in the following discussion of general moiré setups.

5. Moiré effects

We now focus on general moiré geometries, illustrated in Fig. 1. Again, we consider the direction of periodicity to be along the $z$-axis and the arrays to be separated by a distance $d$ in $y$-direction. In contrast to the previous discussion, we now allow the array periodicities to differ. The periodicity of the array at $y=0$ is denoted by $\Lambda _1$. It is related to the periodicity $\Lambda _2$ of the array at $y=d$ by the moiré parameter $\theta$ according to

A commensurable configuration is present if the entire stack can be obtained by translating an appropriate unit cell with an overall periodicity $\Lambda$. Commensurable configurations correspond to values of $\theta$ in the rational numbers, such that there exists a pair of natural numbers $(n,m)$ fulfilling $\theta = m/n$. The periodicity of the composite system $\Lambda$ is then given by $m \Lambda _1 = n \Lambda _2$. Without loss of generality, we assume that $n > m$.

The interaction matrix is thus $3N \times 3N$, where $N = n + m$ is the number of atoms in the unit cell. By Eqns. (7) and (12), we can decompose it into a diagonal contribution resulting from $N$ uncoupled arrays and a coupling matrix $C(k)$, analogous to Eqn. (22) as follows

In the moiré unit cell, we identify $m$ inter-array nearest-neighbor pairs, leaving $n-m$ unpaired atoms. Analogous to our treatment of perfectly aligned arrays, we make the inter-array nearest-neighbor approximation. This implies taking the limit of $d \rightarrow 0$, which diagonalizes all coupling matrices by the equations above. Additionally, we take the limit $d_z \rightarrow 0$ required to define inter-array nearest-neighbors in our approximation reliably. The procedure of isolating singular terms and dropping all non-singular powers of $d$ and $d_z$ in Eqns. (46)–(50) is carried out analogously to the previous section. The asymptotic form of the NN-coupling matrices follows

As before, $\tilde {G}$ depends only on the nearest-neighbor distance. Notice that this asymptotic form only applies to the $2m$ matrices coupling nearest neighbors and not to any of the other matrices entering $C(k)$.

We now employ a partial diagonalization procedure for the entire interaction matrix $M(k)$. First, we subdivide the Cartesian basis of $\mathbb {C}^{3N}$ into $N$ vectors $\textbf{v}^i, i \in \{1, \hdots, N \}$, defined as

The physical intuition behind this notation is that $\textbf{v}^i_a$ corresponds to the $a$-th dipolar excitation component of the $i$-th atom in the moiré unit cell. Then, we can represent the interaction matrix as

The summation over Cartesian components is redundant because all matrices are diagonal, so we omit it in the following. Consider now an arbitrary pair of inter-array nearest neighbors identified by the indices $i_1,i_2$. We perform a change of basis

We recognize this basis as the optical analogon of bonding and antibonding states given for electrons, e.g., in Ref. [17], which is a phenomenon well observed in plasmonics [35,36]. The vector $\textbf{v}_+$ describes in-phase oscillations of the nearest-neighbor pair, while the vector $\textbf{v}_-$ describes out-of-phase oscillations [37,38]. In the limit $d \rightarrow 0, d_{i_1 i_2,z} \rightarrow 0$, the coupling to any dipole in the unit cell indexed by $j$ fulfills $\tilde {G}(k, \textbf{d}_{i_1j}) = \tilde {G}(k, \textbf{d}_{i_2j})$ as well as $\tilde {G}(k, \textbf{d}_{ji_1}) = \tilde {G}(k, \textbf{d}_{ji_2})$ since the nearest neighbors effectively collapse into a single point. It follows that

Despite its unwieldy appearance, this expression has a simple physical interpretation. First, Eqn. (55) is the sum of bonding and antibonding self-interaction matrices. Equation (56) implies that only the bonding state couples to the excitation manifold supported by the remaining $N-2$ atoms in the moiré unit cell. Finally, Eqn. (57) describes the interaction among these remaining $N-2$ atomic dipoles. We conclude that in the inter-array NN-limit a system of $N$ interacting lattice modes is equivalent to a system composed of $N-1$ interacting modes, an additional, completely decoupled mode.

We can thus repeat the procedure detailed above $m$ times for every pair of inter-array nearest neighbors, each time performing the change of basis to bonding and antibonding states. We ultimately arrive at the following partially diagonal form of the system interaction matrix

Here, $M_{\text {flat}}$ is a diagonal matrix, while, in contrast to the matrix in Eqn. (22) for perfectly aligned arrays, the matrix $M _{\text {curved}}(k)$ is not diagonal in general as it describes coupling among bonding states and untransformed lattice modes. The matrix $M_{\text {flat}}$ gives rise to flat bands with infinite lifetimes, thus corresponding to lossless and localized modes. From Eqn. (54), these flat bands can be given a physical meaning in terms of out-of-phase dipolar oscillations of nearest-neighbors in a unit cell. We might thus conceptualize flat band modes as interfering destructively over the small spatial separation between the arrays to produce non-radiating, trapped states.

 

Fig. 4. Approximate low-loss ($G(k) < 0.1$) band structure for strongly coupled $d = 0.1 \lambda _0$ arrays with $\Lambda _1 = 0.3 \lambda _0$. Each finite moiré setup is composed of 30 unit cells to ensure approximate applicability of Bloch’s theorem. Theoretically predicted low-loss flat bands are shown as dashed lines and approximate most of the near-lossless approximately constant bands seen in the numerically determined band structure. As the value of the moiré parameter approaches 1, the number of nearest neighbors in the unit cell is increased and, consequently, the number of flat bands introduced in the spectrum grows.

Download Full Size | PDF

Using Eqns. (9) and (55), we observe that every nearest neighbor pair $i,j$ separated by distance $d_{ij}$ introduces two flat bands given by

The exact value and number of distinct distances $d_{ij}$ depend on the choice of $\theta$. It is worthwhile to contrast these expressions to the ones obtained in the perfectly aligned case given by Eqns. (36)–(38). In the perfectly aligned case, only two distinct spectral regions of enhanced DOS are expected, which is due to two-fold degeneracy in the flat band spectrum. For a moiré setup, however, the occurrence of a large number of flat bands is possible, depending on the number of distinct distances $d_{ij}$ between dipoles in the unit cell. For smaller differences between array periodicities, corresponding to values of $\theta$ closer to $1$, we observe an increasing number of flat bands as the number of nearest-neighbor dipoles in the unit cell increases. We thus recognize the occurrence of these flat bands as a moiré effect enabled by inter-array NN-interactions.

Since these bands are ideally lossless, Fig. 4 highlights the low-loss spectrum of three different moiré setups, characterized by different values of $\theta$. Indeed, the occurrence of most flat bands is correctly predicted by Eqns. (59) and (60). In Fig. 5, we present the ability of these expressions to approximately identify most of the critical spectral regions where the DOS is strongly amplified. Large field enhancement in a desired, albeit narrow, frequency domain may thus be achieved by strong spatial coupling and an appropriate choice of the moiré parameter, where the expressions in Eqns. (59) and (60) can serve as a guiding principle.

 

Fig. 5. DOS for strongly coupled $d = 0.1 \lambda _0$ arrays with $\Lambda _1 = 0.3 \lambda _0$. Each finite moiré setup is composed of 30 unit cells. Lossless flat bands from a theoretical inter-array nearest neighbor approximation are shown as dashed lines and can approximate the location of most of the peaks in the DOS.

Download Full Size | PDF

6. Conclusion

We have studied a 1D moiré geometry of atomic dipoles distributed among two coupled linear arrays in a modified NN-approximation. Our approximation scheme fully considers lattice interactions within a single array, allowing us to capture interactions between different arrays via their nearest-neighbor contributions. Similar to how lossless curved bands arise in an NN-approximation of uncoupled geometries [39], our treatment revealed the occurrence of lossless flat bands in coupled geometries.

The spectral location of these flat bands is given by the distances between inter-array nearest neighbors in a unit cell. Since a moiré unit cell contains multiple inequivalent nearest neighbors, a multitude of flat bands emerge. These flat bands lead to strong amplification of the DOS around the corresponding frequencies. We have numerically confirmed this effect of potentially high field enhancement in the regime of strong geometric coupling, where we expect the NN-contributions to dominate.

We have thus unveiled photonic moiré geometries as a promising platform for highly tunable and spectrally selective field enhancement that might be useful in applications geared towards sensing or tunable light-matter interaction. From a theoretical perspective, our findings provide further insight into the general nature of photonic moiré effects and might aid a better understanding of these effects in more complicated setups regarding scatterer type and geometry.

Further research could focus on the generalization of the above results to 2D atomic dipolar moiré structures, where we expect similar effects to occur and investigate the potentially rich interplay between moiré geometry and other complex phenomena such as the study of topological properties [40], investigating, e.g., topological (non-) triviality of moiré flat bands. Another interesting avenue to explore could be the incorporation of higher-order multipole contributions based on insights from the corresponding single-array geometries, e.g., in [41].

A. Interaction matrix

Here, we closely follow [31] in deriving the interaction matrix for finite geometries, which we then generalize to the infinite periodic case, similar to [32].

As in Sec. 2, assuming a collection of $N$ atoms at resonance, the total electric field at the position $\textbf{r}_i$ of the $i$-th atom is

(61)$$\textbf{E}(\textbf{r}_i) = \frac{k_0^2}{\varepsilon_0} \sum_{j \neq i}^{N} G(k_0, \textbf{r}_i, \textbf{r}_j) \textbf{p}_{j} + \textbf{E}_{\text{ext}}(\textbf{r}_i),$$

where $\textbf{E}_{\text {ext}}(\textbf{r}_i)$ denotes the external field at $\textbf{r}_i$, $\textbf{p}_i$ refers to the $i$-th atomic electric dipole moment, and $G$ is the dyadic Green’s tensor introduced in Eqn. (3). Bringing the sum in Eqn. (61) to the opposite side and using $\textbf{p}_i = \alpha (\omega _0) \textbf{E}(\textbf{r}_i)$ leads to the equation

(62)$$\frac{1}{\alpha(\omega_0)}\textbf{p}_i - \frac{k^2}{\varepsilon_0} \sum_{j \neq i}^{N} G(k_0, \textbf{r}_i - \textbf{r}_j) \textbf{p}_{j} = \textbf{E}_{\text{ext}}(\textbf{r}_i) .$$

Defining $\textbf{p} = \bigoplus\limits _{i=1}^N \textbf{p}_i, \textbf{E}_{\text {ext}} = \bigoplus\limits _{i=1}^N \textbf{E}_{\text {ext}}(\textbf{r}_i)$, this can be written compactly in terms of a $3N \times 3N$ matrix $M$

(63)$$\frac{1}{\epsilon_0 \alpha_0} M \textbf{p} = \textbf{E}_{\text{ext}},$$

where $M$ is the interaction matrix, which can be written in block-form as in Eqn. (2). Notice that, as in Ref. [31], we define the matrix $M$ as a dimensionless quantity and work in units of the constant scaling factor $(\epsilon _0 \alpha _0)^{-1}$ in Eqn. (63), implicitly adopting units where the atomic polarizability in Eqn. (1) is dimensionless. The real and imaginary parts of the eigenvalues of this matrix determine the system resonances and collective decay rates as in Eqns. (4) and (5). In the case of a single dipolar array, discussed in the introductory part of Sec. 4, the interaction matrix takes diagonal form and its eigenvalues are given by its diagonal elements.

The analytical treatment of an infinite system with periodicity $\Lambda$ is largely analogous to the previous discussion. As in Sec. 3, we describe the coupling between $N$ dipoles in a unit cell, including lattice interactions. To this end, we again consider the resonant periodic scattering problem as expressed in Eqn. (61) and divide the periodic setup into unit cells, each containing $N$ particles. One can then apply Bloch’s theorem to connect the $i$-th atomic dipole moment in the $n$-th unit cell to the $i$-th atomic dipole moment in the $0$-th unit cell

(64)$$\textbf{p}_{i,n} = \textbf{p}_{i,0} \mathrm{e}^{\mathrm{i} k n \Lambda},$$

where $k \in [0, \frac {\pi }{\Lambda }]$ is the crystal momentum of the Bloch mode in the irreducible Brillouin zone. The electric field at the position of the $i$-th dipole in the $0$-th unit cell is then

(65)$$\textbf{E}(\textbf{r}_{i,0}) = \frac{k^2}{\varepsilon_0} \sum_{n ={-}\infty}^{\infty} \sum_{j \neq i}^{N} G(k_0, \textbf{r}_{i,0} - \textbf{r}_{j,n}) \textbf{p}_{j,0} \mathrm{e}^{\mathrm{i} k n \Lambda} + \frac{k^2}{\varepsilon_0} \sum_{n \neq 0} G(k_0, \textbf{r}_{i,0} - \textbf{r}_{i,n}) \textbf{p}_{i,0} \mathrm{e}^{\mathrm{i} k n \Lambda} + \textbf{E}_{\text{ext}}(\textbf{r}_{i,0}),$$

Analogous to how we arrived at Eqn. (62), this equation can be cast in matrix form by defining $\textbf{p}_0 = \bigoplus\limits _{i=1}^N \textbf{p}_{i,0}, \textbf{E}_{\text {ext}} = \bigoplus\limits _{i=1}^N \textbf{E}_{\text {ext}}(\textbf{r}_{i,0})$ and arriving at a $3N \times 3N$ interaction matrix for a periodic system $M_P(k)$

(66)$$\frac{1}{\epsilon_0 \alpha_0} M_P(k) \textbf{p}_0 = \textbf{E}_{\text{ext}}.$$

Expressing this matrix in block-form, we obtain Eqn. (7).

B. Simulation method

We adopt an FFT-based approach similar to [42] to construct approximate band structures for sufficiently large finite systems. First, we create a geometry containing $m$ unit cells with $N$ atoms. Denoting by $\textbf{p}^{i,j}$ the $i$-th atomic dipole moment in the $j$-th unit cell, we construct the interaction matrix such that it acts on vectors of the form

(67)$$\textbf{p} = \bigoplus_{i=1}^m \bigoplus_{j=1}^N \textbf{p}^{i,j}.$$

We then subdivide each right eigenvector $\textbf{v}$ of the interaction matrix into $3N$ vectors of length $m$ according to

(68)$$\textbf{w}_{i,a} = \bigoplus_{j=1}^{m} \textbf{v}_a^{i,j},$$

where $i \in \{1, \hdots, N\}$ and $a = 1,2,3$ is the Cartesian component such that each vector $\textbf{w}$ is identified by a unique combination $i,a$. Out of these vectors, we select the critically contributing one identified by

(69)$$(i_c, a_c) = \text{arg max}_{(i,a)} \sum_{j=1}^{m} |\textbf{w}_{i,a;j}|.$$

For this vector, we perform an FFT. From the resulting frequencies, we select and normalize the dominant contribution to compute the approximate crystal momentum $k$.

Since the single atom resonance wavelength $\lambda _0$ enters our treatment merely as a scaling factor without influencing the physics, we set $\lambda _0 = 1$ throughout our simulations. The simulations were conducted using the library JAX for the Python programming language.

C. Validity regime of the perturbative treatment

The perturbative approach detailed in this work is formally similar to expanding the interaction matrix in inverse powers of the inter-array distance $d$ and dropping all terms non-singular in the limit $d \rightarrow 0$. As such, we expect our results to be affected by errors in the order of $\mathcal {O}(d)$. For numerical confirmation, we performed simulations with a moiré setup characterized by $\theta = 3/4$. The lattice constant of the array at $y=0$ is given by $\Lambda _1 = 0.3 \lambda _0$ and we varied the inter-array distance $d$. The results are displayed in Fig. 6. We conclude that our perturbative results provide a mostly valid approximation in the near-field limit.

 

Fig. 6. Approximate low-loss ($G(k) < 0.1$) band structure for strongly coupled $d = 0.1 \lambda _0$ arrays with $\Lambda _1 = 0.3 \lambda _0$ and $\theta = 3/4$. The finite moiré setup is composed of 30 unit cells. The prediction of flat bands according to Eqns. (59) and (60) ceases to be valid once the strong coupling regime is left.

Download Full Size | PDF

D. Robustness against fabrication errors

To investigate the effect of fabrication errors on the optical properties discussed in Sec. 5, simulations involving positional disorder in the sample are conducted. These simulations are performed with a moiré setup characterized by $\theta = 3/4$ and an inter-array distance of $d = 0.1 \lambda _0$. The lattice constant of the array at $y=0$ is given by $\Lambda _1 = 0.3 \lambda _0$. Positional disorder is modeled by randomly displacing each dipole in the geometry. The new random positions are sampled according to a 3D Gaussian centered at the unperturbed dipole position with a standard deviation of $D \cdot d$, where the parameter $D$ is used to quantify the degree of fabrication errors relative to the inter-array separation. The results of these simulations are presented in Fig. 7. The simulations indicate robustness of the predicted moiré flat bands against the effect of positional disorder up to a critical threshold of the disorder parameter $D = 10\%$. Our simulations indicate that for fabrication errors within this tolerance, the desired effect of strong field enhancement will largely persist. Outside of this regime, we observe the broadening and ultimate disappearance of moiré flat bands as the structure deviates too strongly from perfect regularity.

 

Fig. 7. Approximate low-loss ($G(k) < 0.1$) band structure for strongly coupled $d = 0.1 \lambda _0$ arrays with $\Lambda _1 = 0.3 \lambda _0$ and $\theta = 3/4$ under the influence of positional disorder (described in the main text). The finite, perfectly ordered moiré setup is composed of 30 unit cells. The prediction of flat bands according to Eqns. (59) and (60) and the corresponding strong enhancement of the DOS is expected to become invalid for positional disorder parameters of $D >= 10\%$.

Download Full Size | PDF

Funding

Alexander von Humboldt-Stiftung; Helmholtz Association; Carl-Zeiss-Stiftung; Deutsche Forschungsgemeinschaft (390761711); KIT-Publication Fund of the Karlsruhe Institute of Technology.

Disclosures

The authors declare no conflicts of interest.

Data availability

The code to generate the results used in this paper is made publicly available under [43].

References

1. G. Hu, M. Wang, Y. Mazor, C.-W. Qiu, and A. Alù, “Tailoring light with layered and moiré metasurfaces,” Trends Chem. 3(5), 342–358 (2021). [CrossRef]  

2. G. Hu, A. Krasnok, Y. Mazor, C.-W. Qiu, and A. Alù, “Moiré hyperbolic metasurfaces,” Nano Lett. 20(5), 3217–3224 (2020). [CrossRef]  

3. A. Krasnok and A. Alù, “Low-symmetry nanophotonics,” ACS Photonics 9(1), 2–24 (2022). [CrossRef]  

4. R. Röhrich and A. F. Koenderink, “Double moiré localized plasmon structured illumination microscopy,” Nanophotonics 10(3), 1107–1121 (2021). [CrossRef]  

5. P. Wang, Y. Zheng, X. Chen, C. Huang, Y. V. Kartashov, L. Torner, V. V. Konotop, and F. Ye, “Localization and delocalization of light in photonic moiré lattices,” Nature 577(7788), 42–46 (2020). [CrossRef]  

6. L. Huang, W. Zhang, and X. Zhang, “Moiré quasibound states in the continuum,” Phys. Rev. Lett. 128(25), 253901 (2022). [CrossRef]  

7. Z. Wu and Y. Zheng, “Moiré metamaterials and metasurfaces,” Adv. Opt. Mater. 6(3), 1701057 (2018). [CrossRef]  

8. S. Liu, S. Ma, R. Shao, L. Zhang, T. Yan, Q. Ma, S. Zhang, and T. J. Cui, “Moiré metasurfaces for dynamic beamforming,” Sci. Adv. 8(33), eabo1511 (2022). [CrossRef]  

9. N. S. Salakhova, I. M. Fradkin, S. A. Dyakov, and N. A. Gippius, “Fourier modal method for moiré lattices,” Phys. Rev. B 104(8), 085424 (2021). [CrossRef]  

10. J. Sperrhake, M. Falkner, S. Fasold, T. Kaiser, and T. Pertsch, “Equivalence of reflection paths of light and Feynman paths in stacked metasurfaces,” Phys. Rev. B 102(24), 245108 (2020). [CrossRef]  

11. C. Menzel, J. Sperrhake, and T. Pertsch, “Efficient treatment of stacked metasurfaces for optimizing and enhancing the range of accessible optical functionalities,” Phys. Rev. A 93(6), 063832 (2016). [CrossRef]  

12. X. Romain, F. I. Baida, and P. Boyer, “Spectrally tunable linear polarization rotation using stacked metallic metamaterials,” J. Opt. 19(8), 085102 (2017). [CrossRef]  

13. B. Lou, N. Zhao, M. Minkov, C. Guo, M. Orenstein, and S. Fan, “Theory for twisted bilayer photonic crystal slabs,” Phys. Rev. Lett. 126(13), 136101 (2021). [CrossRef]  

14. H. Tang, F. Du, S. Carr, C. DeVault, O. Mello, and E. Mazur, “Modeling the optical properties of twisted bilayer photonic crystals,” Light: Sci. Appl. 10(1), 157 (2021). [CrossRef]  

15. H. Z. Zhang, H. Y. Qin, W. X. Zhang, L. Huang, and X. D. Zhang, “Moiré graphene nanoribbons: Nearly perfect absorptions and highly efficient reflections with wide angles,” Opt. Express 30(2), 2219 (2022). [CrossRef]  

16. W. Zhang, D. Zou, Q. Pei, W. He, H. Sun, and X. Zhang, “Moiré circuits: Engineering magic-angle behavior,” Phys. Rev. B 104(20), L201408 (2021). [CrossRef]  

17. S. Carr, D. Massatt, M. Luskin, and E. Kaxiras, “Duality between atomic configurations and bloch states in twistronic materials,” Phys. Rev. Res. 2(3), 033162 (2020). [CrossRef]  

18. D. X. Nguyen, X. Letartre, E. Drouard, P. Viktorovitch, H. C. Nguyen, and H. S. Nguyen, “Magic configurations in moiré superlattice of bilayer photonic crystals: Almost-perfect flatbands and unconventional localization,” Phys. Rev. Res. 4(3), L032031 (2022). [CrossRef]  

19. R. Alaee, B. Gurlek, M. Albooyeh, D. Martín-Cano, and V. Sandoghdar, “Quantum metamaterials with magnetic response at optical frequencies,” Phys. Rev. Lett. 125(6), 063601 (2020). [CrossRef]  

20. K. Ballantine and J. Ruostekoski, “Optical magnetism and huygens’ surfaces in arrays of atoms induced by cooperative responses,” Phys. Rev. Lett. 125(14), 143604 (2020). [CrossRef]  

21. M. Martí-Sabaté and D. Torrent, “Dipolar localization of waves in twisted phononic crystal plates,” Phys. Rev. Appl. 15(1), L011001 (2021). [CrossRef]  

22. J. Zeng, Y. Hu, X. Zhang, S. Fu, H. Yin, Z. Li, and Z. Chen, “Localization-to-delocalization transition of light in frequency-tuned photonic moiré lattices,” Opt. Express 29(16), 25388–25398 (2021). [CrossRef]  

23. V. E. Babicheva and A. B. Evlyukhin, “Analytical model of resonant electromagnetic dipole-quadrupole coupling in nanoparticle arrays,” Phys. Rev. B 99(19), 195444 (2019). [CrossRef]  

24. E. Shahmoon, D. S. Wild, M. D. Lukin, and S. F. Yelin, “Cooperative resonances in light scattering from two-dimensional atomic arrays,” Phys. Rev. Lett. 118(11), 113601 (2017). [CrossRef]  

25. A. Lagendijk and B. A. Van Tiggelen, “Resonant multiple scattering of light,” Phys. Rep. 270(3), 143–215 (1996). [CrossRef]  

26. K. D. Bonin and V. V. Kresin, Electric-dipole Polarizabilities of Atoms, Molecules, and Clusters (World Scientific, 1997).

27. A. Rahimzadegan, R. Alaee, C. Rockstuhl, and R. W. Boyd, “Minimalist Mie coefficient model,” Opt. Express 28(11), 16511 (2020). [CrossRef]  

28. J. Rui, D. Wei, A. Rubio-Abadal, S. Hollerith, J. Zeiher, D. M. Stamper-Kurn, C. Gross, and I. Bloch, “A subradiant optical mirror formed by a single structured atomic layer,” Nature 583(7816), 369–374 (2020). [CrossRef]  

29. L. Novotny and B. Hecht, Principles of Nano-Optics (Cambridge University Press, 2012), 2nd ed.

30. V. E. Babicheva and A. B. Evlyukhin, “Metasurfaces with electric quadrupole and magnetic dipole resonant coupling,” ACS Photonics 5(5), 2022–2033 (2018). [CrossRef]  

31. S. E. Skipetrov, “Finite-size scaling of the density of states inside band gaps of ideal and disordered photonic crystals,” Eur. Phys. J. B 93(4), 70 (2020). [CrossRef]  

32. A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht, H. J. Kimble, and D. E. Chang, “Exponential improvement in photon storage fidelities using subradiance and “selective radiance” in atomic arrays,” Phys. Rev. X 7(3), 031024 (2017). [CrossRef]  

33. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (Dover, 1964).

34. A. Alù and N. Engheta, “Theory of linear chains of metamaterial/plasmonic particles as subdiffraction optical nanotransmission lines,” Phys. Rev. B 74(20), 205436 (2006). [CrossRef]  

35. F. Kretschmer, S. Muehlig, S. Hoeppener, A. Winter, M. D. Hager, C. Rockstuhl, T. Pertsch, and U. S. Schubert, “Survey of plasmonic nanoparticles: from synthesis to application,” Part. Part. Syst. Charact. 31(7), 721–744 (2014). [CrossRef]  

36. H. Wang, D. W. Brandl, P. Nordlander, and N. J. Halas, “Plasmonic nanostructures: artificial molecules,” Acc. Chem. Res. 40(1), 53–62 (2007). [CrossRef]  

37. N. Harris, M. D. Arnold, M. G. Blaber, and M. J. Ford, “Plasmonic resonances of closely coupled gold nanosphere chains,” J. Phys. Chem. C 113(7), 2784–2791 (2009). [CrossRef]  

38. W. A. Murray and W. L. Barnes, “Plasmonic materials,” Adv. Mater. 19(22), 3771–3782 (2007). [CrossRef]  

39. R. P. H. Wu and H. C. Ong, “Small mode volume topological photonic states in one-dimensional lattices with dipole-quadrupole interactions,” Phys. Rev. Res. 4(2), 023233 (2022). [CrossRef]  

40. R. J. Bettles, J. Minář, C. S. Adams, I. Lesanovsky, and B. Olmos, “Topological properties of a dense atomic lattice gas,” Phys. Rev. A 96(4), 041603 (2017). [CrossRef]  

41. A. Alù and N. Engheta, “Guided propagation along quadrupolar chains of plasmonic nanoparticles,” Phys. Rev. B 79(23), 235412 (2009). [CrossRef]  

42. Y. Chen, M. Kadic, and M. Wegener, “Roton-like acoustical dispersion relations in 3D metamaterials,” Nat. Commun. 12(1), 3278 (2021). [CrossRef]  

43. D. Dams, “Lorentzian band structures,” Github, 2023, https://github.com/tfp-photonics/dipole_band_structures.