Playing relativistic billiards beyond graphene

The possibility of using hexagonal structures in general and graphene in particular to emulate the Dirac equation is the basis of our considerations. We show that Dirac oscillators with or without restmass can be emulated by distorting a tight binding model on a hexagonal structure. In a quest to make a toy model for such relativistic equations we first show that a hexagonal lattice of attractive potential wells would be a good candidate. First we consider the corresponding one-dimensional model giving rise to a one-dimensional Dirac oscillator, and then construct explicitly the deformations needed in the two-dimensional case. Finally we discuss, how such a model can be implemented as an electromagnetic billiard using arrays of dielectric resonators between two conducting plates that ensure evanescent modes outside the resonators for transversal electric modes, and describe an appropriate experimental setup.

Abstract. The possibility of using hexagonal structures in general, and graphene in particular, to emulate the Dirac equation is the topic under consideration here. We show that Dirac oscillators with or without rest mass can be emulated by distorting a tight-binding model on a hexagonal structure. In the quest to make a toy model for such relativistic equations, we first show that a hexagonal lattice of attractive potential wells would be a good candidate. Firstly, we consider the corresponding one-dimensional (1D) model giving rise to a 1D Dirac oscillator and then construct explicitly the deformations needed in the 2D case. Finally, we discuss how such a model can be implemented as an electromagnetic billiard using arrays of dielectric resonators between two conducting plates that ensure evanescent modes outside the resonators for transversal electric modes, and we describe a feasible experimental setup. 5 Author to whom any correspondence should be addressed.

Introduction
Dirac operators play a central role in relativistic quantum dynamics, from the early work of Dirac and the exact solutions for the hydrogen atom to later work including quantum chromodynamics and the Dirac oscillator [1]- [3]. Based on the early work of Wallace [4], there has been much work recently that is centered around the fact that the mean field theory of graphene is well represented by the free Dirac equation near the edges of the first Brillouin zone, i.e. at the center of the band with so-called Dirac points around which linear dispersion relations hold [5]- [7]. Simulations of this kind of situations are ongoing at various laboratories using hexagonal and occasionally triangular arrays in single-particle problems or classical waves, including acoustics [8], microwaves [9,10] and photonic crystals [11] as well as true quantum simulations in nanostructures [12]- [18]. We wish to recall that any lattice with coordination number three will yield points with approximate linear dispersion relations similar to Dirac points, but we need a tight-binding situation to guarantee isotropic cones and a hexagonal structure to introduce an additional discrete degree of freedom for the small and large components of the effective spinor. This becomes very transparent as we note that by using a hexagonal structure made up of two different triangular structures, we obtain a finite gap in the spectrum, i.e. a Dirac equation for a massive particle, a situation that would correspond to a B-N (boron nitride) lattice. It is important to note that emulations of these structures in microwave billiards are restricted to reproducing a single fermion propagating in a material, since Fermi-Dirac statistics cannot be emulated by electromagnetic waves. However, a wave can explore all the energy states analogous to the single-particle electronic states in a solid.
In this paper, we will address the task to find similar analogues for Dirac operators, which describe situations other than the free particle. Formally, we shall assume arrays of potential wells that can hold exactly one bound state. In between the wells, these states decay exponentially. As the overlap of the functions of two wells describes the coupling, this implies to a good approximation a tight-binding system. In this framework, in principle one-and two-dimensional (1D and 2D) Dirac-type problems can be formulated and in some cases exactly solved. We shall here focus on the Dirac oscillator as introduced by Moshinsky and Szczepaniak [3]; this system is frequently considered as a paradigm of integrability in relativistic quantum mechanics [19]- [25], [54] and has been used to model fermions in confining potentials as well as composite relativistic systems such as hadrons [26]- [28]. We shall see that it is particularly simple to implement. Other solvable problems, such as gyroscopes [29] and the Coulomb problem as well as systems with random potentials, will be touched upon in the conclusions.
In particular, the 1D and 2D Dirac oscillators have been considered recently in the context of both relativistic quantum mechanics and quantum optics [30]- [32]. Their importance as paradigmatic integrable models is well established and their realization in the context of tightbinding models is desirable. We want to stress that the dimensionality in our examples plays a crucial role in their algebraic properties and spectra. As we shall see, the exceptional infinite degeneracy of the 2D Dirac oscillator (as opposed to the finite degeneracy in 1D) is intimately related to its specific realization on the lattices described above. It is therefore interesting to consider both dimensionalities, as they exhibit clear differences. It remains an open question whether 3D relativistic wave equations can be emulated on a 3D lattice.
The description of such systems will depend on a massive distortion of the hexagonal lattice. This causes probably prohibitive obstacles to any implementation in terms of mean fields of graphene-related structures. For a recent theoretical work dealing with deformations on carbon sheets, see [33]. Other works dealing with deformed graphene in connection with integrable systems such as Landau electrons can be found in the litertature [34,35]. For classical wave models, on the other hand, such distortions can be implemented and we shall discuss specifically how this is done in what we usually call a microwave billiard. In two dimensions, such systems yield scalar wave equations and can readily be interpreted as single-particle quantum systems [36].
Microwave billiards have been used to simulate a wide variety of phenomena including quantum chaos [37]- [40], scattering from open billiards [41,42], transport phenomena [43,44], fidelity decay [45,46] and disordered systems [47]- [49], as well as recent work on Dirac points [9,10]. This wide range of success combined with recent experiments with arrays of small dielectric microwave resonators gives us a good reason to hope for interesting results in ordered structures such as the ones we need. The evanescence of the cavity modes that we intend to use guarantees exponential decay and thus tight-binding situations can be emulated.
In the next section, we study a one-dimensional crystal; this will serve to fix the notation and basic ideas in a simpler context, which is not without relevance in itself. This toy model reveals Dirac-like Hamiltonians in both periodic and deformed structures. In section 3, we analyze and reproduce similar results for the 2D case with and without deformations. We also consider corrections to nearest-neighbor interactions and obtain the form of energy surfaces beyond tight binding. In section 4, we access the experimental applications of our treatment by considering arrays of resonators in microwave cavities. The applicability of the tight-binding model in this context is discussed. Finally, we draw some conclusions and give an outlook on a selection of other Dirac systems that can be emulated. Useful results are included in the appendix.

1D crystal
In this section, we start with two lattices. First, we consider the Schrödinger equation with square wells supporting a single bound state, located in a periodic 1D array. Here, traditional aspects of the existing theory (i.e. Bloch waves) are reviewed under an algebraic approach. The specific form of the localized wave functions is ignored, as it turns out to be irrelevant. An analogy between degeneracy points of the spectrum and a 1D Dirac point is revealed in this toy model. Near these points we shall cover a one-dimensional Dirac equation.
Once this is achieved, we proceed to deform the lattice from its periodic configuration by forcing a specific operator algebra, namely the one corresponding to Harmonic oscillator ladder operators. As an outcome, a 1D Dirac oscillator shall be realized.

The Dirac point in one dimension
We start by fixing some notation. In the following, we adopt natural unitsh = c = 1. A lattice consisting of two periodic sublattices is considered. Let λ be the distance between neighboring potential wells. They have the same period and are denoted as types A and B. Each sublattice point can be labeled by an integer n according to its position on the line, i.e. x n for type A and y n for type B (see figure 1). The energy of the single level to be considered in the well is denoted by α and β for types A and B, respectively. The state corresponding to a particle at site n of lattice A is denoted by |n A and the corresponding localized wave function is given by ξ A (x − x n ) = x|n A . For lattice B we define the wave function ξ B (x − y n ) = x|n B . Our Hamiltonian is well approximated by a tight-binding model if the overlap between localized wave functions can be neglected and the nearest-neighbor coupling matrix element is taken as the first off-diagonal element of the Hamiltonian in the localized basis.
For convenience we have split the lattice into two sublattices and we write the Hamiltonian in the form rather than in the usual tridiagonal form. The entries of each block are given by H n,m i j = n|H i j |m with n, m being integers. Thus, each block extends from −∞|H i j | − ∞ to ∞|H i j |∞ . In this basis, we write the Hamiltonian as where the elements outside the indicated diagonals are all zero. Expression (2) can be cast in terms of Pauli matrices σ 3 , σ + = σ 1 + iσ 2 , σ − = σ † + by defining = H AB ⊗ 1 and setting M = (α − β)/2, E 0 = (α + β)/2. We have displaying explicitly the Dirac-like structure of the Hamiltonian. To ensure the analogy between the Dirac Hamiltonian and (3), we consider the spectrum of H . We note that has the remarkable properties From these relations we may compute the spectrum by squaring H Note that Bloch's theorem manifests itself in the spectrum and eigenfunctions of † as Therefore, energies and eigenfunctions of H are given by where N is a normalization constant. When M = 0, i.e. when sites A and B are equal, we can find conical points of the form k 0 = ±1/(2λ). Expanding around these points we find the usual expressions for the relativistic energies of Dirac particles with momentum κ = k − k 0 , effective speed of light : This is valid also for nonzero rest energy M in the case of different lattices. Then we find a gap in the spectrum. The eigenfunctions satisfy the Dirac equation in momentum space In the next subsection, we shall proceed to show that we can obtain the Dirac oscillator by deforming the double lattice.

The 1D Dirac oscillator and lattice deformations
Using the notation of the previous section, we modify the energy operator by deforming the lattice of potential wells. This implies abandoning the periodic structure and leads to a site dependence of the couplings in the corresponding tight-binding model. We denote by n,n+1 the coupling between sites at y n and x n+1 , while n,n denotes the coupling between sites x n and y n . These couplings are approximately proportional to the overlap between neighboring sites. They decay exponentially as a function of the separation distance between the potential wells, i.e.
where d n and d n are the deviations from the periodic configuration, i.e. d n + λ = |y n+1 − x n |, d n + λ = |y n − x n |. When d n = d n = 0, the periodic configuration is recovered (see figure 1). The length is the penetration depth into the classically forbidden region for the wave function of the single well. We expect a modification of the operators , † caused by the change in the position of the wells. One should keep in mind that such deformations have the effect of breaking the periodic symmetry of the system and Bloch wave functions cease to be useful. However, if one finds an exactly solvable model for a deformation, which is continuously connected to the periodic configuration, then the corresponding solutions could be considered to constitute a generalization of Bloch's waves.
The algebraic properties of observables in Hamiltonians have a clear connection with integrability and exact solvability. Therefore, the simplest way to extend our Hamiltonian is by replacing the operators , † by a, a † , such that their algebraic properties are those of known solvable systems. In particular, we propose the Harmonic oscillator algebra, since it is a paradigmatic example [51]. The Hamiltonian (3) becomes 7 We require that such extensions reduce to the periodic case in some free limit. We propose then and impose the condition where ω stands for a frequency. If ω = 0, we recover the algebra of , † (the Bloch limit).
Computing the commutator in (14), one finds the conditions where the first equality is a recurrence relation solved by the constant , whereas the second equality is a recurrence relation solved by and thus Imposing the conditions (11) in this solution, we find The first condition implies d n = 0 (a chain of dimers). The second condition relates the displacement between resonators d n with the corresponding couplings. For convenience we impose a third condition, namely the inequality This condition ensures that the displacement d n increases monotonically as a function of the distance from an arbitrarily chosen origin. Combining it with the second condition, this implies a finite grid, which in turn leads to a finite Hilbert space of dimension n max + 1 with The operator a takes the finite matrix form In principle, other non-monotonic choices of scaling d n are possible and lead to other finite or infinite arrays. For our modeling purposes finite arrays are advantageous, since they can be easily implemented in a real situation. Moreover, a localization effect of wave functions around the origin of the array is possible due to an increasing distance of resonators. The analogy with a particle trapped by a potential is thus realized and the explicit form of wave functions can be 8 found in the appendix. In any case our Hamiltonian will emulate a Dirac oscillator only in the vicinity of zero energy (energies near E 0 ) as we will see later.
We now have to check the commutation relations between a and a † [a , Note that the correction term is of order 1/n max as the principal term is accompanied by the identity, while the correction acts only in the first component of this basis. The prefactor is of Now we can replace a by a and return to the Hamiltonian (12) for the finite system. Such a Hamiltonian allows analytical solutions for energies and eigenfunctions. Our approximate construction of (12) allows us to compute the spectrum and eigenfunctions up to a shift smaller than unity. As we shall see, they correspond to those of the Dirac oscillator [3]. We consider eigenvectors φ n of the number operator such that a † a φ n = (ω )nφ n , a † φ n = √ ω (n + 1)φ n+1 and aφ 0 = 0. The Hamiltonian has an integral of the motion given by the operator These blocks can be diagonalized, leading to the energies An additional 1 × 1 block is due to the singlet The eigenfunctions corresponding to the doublets are obtained in the form where N is a normalization constant. The specific form of φ n is given in the appendix. We would like to emphasize that the validity of these solutions is only restricted by our approximate construction of the modeling Hamiltonian. Once the tight-binding approximation is imposed and the O(1/n max ) corrections are neglected, the resulting energies and wave functions are analytical solutions of the stationary problem. The usual representation of a Dirac oscillator in terms of momentum and a confining potential is obtained by the transformation P = i(a + a † ) − k 0 , X = (a + a † ). In the Bloch limit, ω = 0 and P reduces to the momentum operator around the conical point.
In the next section, we shall see that a similar construction occurs very naturally if we deform hexagonal lattices in two dimensions, which in the undeformed case, are closely related to the mean field theory of graphene and B-N sheets.

The 2D crystal
The concepts given in the last section are now extended to two dimensions. We shall consider hexagonal structures that emulate 2D Dirac equations and can emulate the mean field of graphene or B-N near the gap at the center of the usual band. We shall use the same algebraic strategy to derive spectra and a possible extension through deformations, namely the 2D Dirac oscillator.

Hexagonal lattice
Let us fix the notation for this system. The honeycomb lattice is divided in to two triangular sublattices, one of them generated by the set of vectors a 1 = (3 √ 2, 0), a 2 = (−3/ √ 2, 3/2) and a 3 = (3/ √ 2, −3/2) (type A) while the other sublattice is obtained by adding the vectors b 1 = (0, 1), b 2 = (−3/ √ 2, −1/2) and b 3 = (3/ √ 2, −1/2). These vectors are given in arbitrary units (see figure 2). We denote the linear combinations of a i by A, where A is a vector parameterizing the points of sublattice A. For sublattice B we use the vector parameter A + b 1 . The position vectors r A , r B of the periodic lattices are obtained by introducing the factor λ. For periodic arrays this means r A = λA, r B = λB. In further considerations this notation will be useful, since deformed lattices admit a parameterization by vectors A, B, but the corresponding position vectors r A , r B become more complicated functions of a i , b i . The state vectors for individual potential wells on grid A shall be denoted by |A , giving wave functions of individual wells as ξ A (r − r A ) = r|A . For grid B we use |A + b 1 . As before, we consider different energies for the wells on grids A and B.
Similarly as in the 1D case, the Hamiltonian in the tight-binding approximation is constructed as where the first two terms indicate the total on-site energy on grids A and B, respectively, while the last sum indicates the nearest-neighbor interaction with coupling strength . We shall analyze this system by considering again a subdivision of the Hilbert space according to sublattices A and B. Due to the coordination number in this lattice, the matrix representation in section 2 is no longer feasible. However, we may construct the usual Pauli operators through the definitions while the operators , and † are defined through They have the algebraic properties With M and E 0 given as in section 2, we obtain the Hamiltonian and Using Bloch waves, we have eigenvectors of the form [50] φ a k = for grids A and B, respectively. They satisfy † φ a,b The spectrum and the eigenfunctions are then The degeneracy points of the spectrum for the massless case are k 0 = ± 1 2λ (1, − √ 3). Expanding around such points one finds as expected. One can verify that the Dirac equation is again satisfied.

The importance of the tight-binding approximation
We have formulated a theory describing the propagation of waves in hexagonal arrays of resonators with and without deformations. The treatment has been successful in predicting the existence of Dirac points, in compliance with the common knowledge about this system when periodic symmetry is present. However, one may ask whether the tight-binding approximation is an essential ingredient for a Dirac-like behavior of the propagating wave.
To answer this, let us review some of the assumptions we made and the corresponding properties obtained for the 2D Dirac wave function. The spin emerges as the probability of the electrons to be located at sites of the triangular sublattice A (spin up) or B (spin down). The momentum (or wave vector) as a conserved quantity came directly from Bloch's theorem, i.e. the periodicity of the system. Linearity around degeneracy points (which is essential to produce a Dirac Hamiltonian) is a consequence of the hexagonal structure, while the existence of such degeneracy points came from the symmetry under the interchange of the two sublattices. We have seen that the appearance of mass corresponds to the lifting of such degeneracy.
In addition to all these properties, one should consider isotropy as a fundamental requirement for the emulation of a free Dirac particle. We shall see that rotational symmetry around degeneracy points is a direct consequence of the tight-binding approximation. It is well known that rotational symmetry in the Dirac equation demands a transformation of both orbital and spinorial degrees of freedom. It is on the orbital part that we shall concentrate by studying the energy surfaces around degeneracy points beyond the tight-binding model.

Let us recall that the transition amplitudes between nearest neighbors (hopping transition) gave rise to the Hamiltonian
with T b i being the one-site translation operator in the direction of b i . Since degeneracy points are located using the condition H ψ 0 = 0 for a Dirac state ψ 0 , and given the fact that [ , † ] = 0, it suffices to impose ψ 0 = 0. The operators T b i are unitary and commute with each other, implying that their eigenvalues can be found simultaneously as e iλk·b i , where k is any real wave vector. Small deviations from degeneracy points (denoted by k 0 ) in the form k = k 0 + κ give the energy which is rotationally invariant in κ. However, one may try to introduce interactions between sites separated by more than one step under lattice translations. It is clear that a second-neighbor interaction of strength modifies the kinetic operator as where the vectors a i have now appeared, connecting a point with its six second neighbors. The energy equation becomes We expect a deviation of degeneracy points k 0 , for which k = k 0 + κ. Upon linearization of the exponentials in κ, we find the energy where the vectors are given by Thus, the presence of yields the energy surfaces (42) as cones with elliptic sections whenever κ is inside the first Brillouin zone. Regardless of how we complete the energy contours to recover periodicity, it is evident that the resulting surfaces are not invariant under rotations around degeneracy points. The circular case is recovered only when = 0, leading to k 0 = k 0 . In this case, the vectors reduce to v = (1, 0), u = (0, 1) when k 0 is the degeneracy point at (1/2λ, 0). In summary, extending the interactions to second neighbors has the effect of breaking the isotropy of space around degeneracy points, which is an essential property of the free Dirac theory.

The 2D Dirac oscillator
Before analyzing lattice deformations, let us recall [24,30] that the two-dimensional Dirac oscillator Hamiltonian is quite similar to the 1D case, except for the replacement a → a R , where a R = a x + ia y is the chiral (right) annihilation operator in terms of Cartesian annihilation operators a x , a y . The corresponding number operator is a conserved quantity and is given by N R = N − L, i.e. the difference between the total number operator and the orbital angular momentum. Since the chiral left operator is absent in the expression, the spectrum is infinitely degenerate. Eigenfunctions are constructed as a combination of two usual Harmonic oscillator functions with defined orbital angular momentum. Our aim is to produce the spectrum for this problem in the hexagonal array, together with a deformation that allows localization of the wave functions around some center.
We proceed to deform the lattice through an extension of the kinetic operators, in analogy to the 1D case. Let us consider site-dependent couplings (A, A + b i ) connecting the sites labeled by A, A + b i . Again, these are related to distances d(A, A + b i ) between potential wells as (A, A + . We define the ladder operator and impose [a R , a † R ] = ω . After some algebra, one can prove that this leads to the conditions The vector b 1 in the first equation was chosen arbitrarily, but due to symmetry a choice of b 2 or b 3 would be equivalent. Due to the coordination number three, we obtain three relations rather than the two of the 1D case. As complicated as the recursion relations may seem, one can easily construct a lattice reproducing them consistently. Relations (46) and (47) establish an equality between the lengths of opposite sides of a given hexagon. Relation (48) containing ω gives the deformation and can be split in to two parts where θ is an arbitrary angle. Having chosen previously the privileged direction b 1 , the angle θ determines the relative stretching between directions b 2 and b 3 on the seminal cell. The choice θ = 0 produces deformations only in one direction of the lattice (b 3 ), resembling the 1D case discussed above. This leads to a logarithmic law for the deviation distance similar to (11). Logarithmic stretching will hold for an arbitrary angle θ . To construct the grid in the general case, one starts with a regular hexagon as a seed and completes the scheme in figure 3 by extending lines of equal length in the direction of b 1 . Then, one completes the hexagons by drawing parallel lines for the opposite sides as shown also in figure 3. Hexagonal cells satisfy the recursion relations above trivially. A restriction to a finite-dimensional space occurs in a way similar to the one-dimensional case. The maximum number of levels is N max = [| ω |]. Considering a R as the chiral operator restricted to this finite-dimensional space, the resulting Hamiltonian of this problem is As the Hamiltonian (51) is formally identical to that in the 1D case, we find the eigenvalues The shape of the eigenfunctions is obtained by solving a R φ 0 = 0 and applying the raising operators similar to those in the appendix. The Hamiltonian does not depend on the left operators a L , a † L , where a L = (a R ) * . In the full space this would imply a Landau electron-like infinite degeneracy. The degeneracy does not occur for a fixed array of resonators, but it can be interpreted to reflect the arbitrary choice of θ if we consider the hypothetical use of an ensemble of arrays for all angles θ .
Note that there is a physical limitation to the allowed degree of distortion, which results from the fact that the nearest neighbors can change. Before this happens, the coupling of potential wells can no longer be dominated by the original three nearest neighbors.
In figures 4-6 we give some realizations of lattice deformations following the procedure indicated previously. The examples are related to the choices of θ = 0, π/4 and π/2, which determine the vectors to be deformed with a logarithmic law near the x-axis of the graphs.

Experimental implementation in microwave cavities
We shall now underline the feasibility of such an emulation of the Dirac equation by making a specific suggestion for an implementation, although by no means the only possible one. To emulate Dirac-like equations on 2D arrays of resonators, we have to generate a scalar field. For electromagnetic fields, this can be achieved in a 2D metallic cavity that supports two types of independent modes: transverse magnetic (TM) mode, ψ(r) = E z (r), and transverse electric (TE) mode, ψ(r) = H z (r), r lying in the plane of the cavity. For a cavity of height h (in the z-direction) and for frequencies ν < c/2h, the scalar field obeys the Helmholtz equation The experimental setups using chaotic microwave cavities adopt various configurations depending on the specific studies they are intended for, but the technique to feed microwaves into the cavity and to collect the signal of interest is ubiquitous. The central device is the network analyzer, which performs emission and lock-in detection of microwaves from a few tens of MHz to a few tens of GHz. The cavity is linked to the network analyzer through flexible coaxial 50 cables connected to monopolar antennas whose central conductor penetrates into the cavity. Usually, several antennas are dispatched over the top and/or bottom plates. For a measurement, only one antenna is used at a time as a microwave emitter and another (in transmission) or the same (in reflection) as a receiver. The other unused antennas are terminated by 50 loads so that all antennas behave the same way regarding the losses they imply. The measurements are given in terms of scattering coefficients, which form the complex S-matrix where S 11 (respectively S 22 ) measures the reflection on port 1 (respectively 2) and S 12 (respectively S 21 ) measures the transmission from port 2 (respectively 1) to port 1 (respectively 2). All the measurements are performed after a proper calibration to get rid of any influence from the cables and connectors and even from the analyzer itself. In a microwave cavity, a varying potential can be obtained by introducing substances of varying permittivity (r). The wave equation then reads Note that the effective potentialṼ (r) = (1 − (r))k 2 is energy dependent. This does not preclude the quantum-classical analogy.
Recently, one of us developed experiments implementing equation (55) in a disordered microwave cavity [47]. The physical phenomenon put under scrutiny was Anderson (or strong) localization (see [52] for a recent survey of this prolific domain). A network analyzer Rohde and Schwarz ZVA-24 was used in a frequency range from 1 to 10 GHz. The disordered potential was introduced through 200 dielectric cylinders (Temex-Ceramics, E2000 series) of high dielectric permitivity ( = 37) and low loss (quality factor Q = 7000 at 7 GHz). Their height fitted that of the cavity, 5 mm, and several diameters ranging from 6 to 8 mm were used. The disorder was numerically generated and the scatterers were precisely positioned. In this experiment, the central conductor of the antennas was perpendicular to the plane of the cavity. Thus, the TM polarization of the electromagnetic field was selected [47,53].
The same dielectric cylinders, used for the localization experiments, can be arranged in periodic or other ordered patterns. In the domain of optics, periodic structures in semiconductor materials are widely used to obtain the particular transport properties of guided light. Thanks to photonic crystals (also called photonic band gap materials) the technology of photons conspires to supplant the one of electrons in domains such as communication and information technologies, computing and sensing [11]. Microwave cavity experiments, which are definitely intended for more fundamental issues, bear the benefits of their versatility. They allow, for example, distortions destroying the periodicity of the potential such as the ones described in sections 2 and 3 for emulating Dirac oscillators.
As emphasized in section 3.2, we have to apply the tight-binding condition. In the localization experiment, the field filled all the cavity, TM polarization being supported inside and outside the resonators. For the experiments we proposed with ordered structures, the requirements are quite different. In order to transmit the energy efficiently into the cylinders, we use in-plane antennas and consequently TE polarization. Due to the wavelength reduction by a factor of √ 6 inside the dielectric material, TE modes can be excited into the scatterers above 5 GHz, while 30 GHz is required in air between the resonators. In the proposed frequency range, resonators support modes that decay evanescently in the surrounding space. This implies exponential decay outside the resonators. For a dielectric cylinder of 8 mm diameter, the first TE mode appears at 6.66 GHz. We experimentally checked that this mode is isotropic: ψ(r) ∼ J 0 (r), thus enabling isotropic coupling between resonators. We studied the range of coupling and its spatial dependence. Measurements were done with 2 and 3 resonators equally spaced on a line, and 6 resonators placed at the vertices of a hexagon (a benzene-like structure). For all configurations, we observed an exponentially decreasing coupling with a characteristic length of 400 m −1 . The three-cylinder and benzene measurements clearly established that the second-neighbor coupling is negligible for distances between the centers of the scatterers above 10 mm. This gives a large range of the coupling constant for which the tight-binding model is fulfilled. All these experimental results will be published in a forthcoming paper focused on transport properties in graphene-like structures [9].
As an advance of how the Dirac oscillator shall be realized, we consider resonators placed at the vertices of the arrays given in figures 4-6, corresponding to different states with the same theoretical energy levels (extraordinary degeneracy of the 2D Dirac oscillator). Two types of resonators are needed to produce a mass. In order to introduce a second resonance coming from type B cylinders, we must ensure that there is no overlap with other resonances of type A cylinders, located nearly at 5 and 8 GHz according to our experimental results. This suggests a mass M proportional to a frequency difference from 10 to 100 Mhz. For separation distances of 3 mm we expect a coupling ∼ 70 Mhz. The frequency ω is fixed by the decrease of the coupling as a function of the deformation distance. This can be obtained from the second formula in (19) for n = 1 or from the corresponding relation for the 2D case. A deformation distance of 1 mm between nearest resonators would produce ω ∼ 0.1 ∼ 7 MHz and a number of levels n max ∼ 10 in a window from 6.4 to 7 GHz.
We have thus established that we can meet the conditions for the emulation of the Dirac oscillator and related problems can be met with arrays of dielectric microwave oscillators between two conducting plates. These conditions will also allow us to emulate other Dirac operators corresponding to gyroscopes, disordered systems, etc. Note though that it will be difficult to emulate a Dirac hydrogen atom as the distances between resonators should become extremely small, causing conflict with the diameter of the discs.

Conclusions and outlook
We have proposed that a wide class of relativistic equations of the Dirac type can be implemented by arrays of potential wells. Such quantum systems can be well approximated by tight-binding Hamiltonians if sufficiently fast exponential decay of the wave functions in the classically prohibited region is ensured. This idea was specifically implemented for the Dirac oscillator. Furthermore, we have shown by way of example that we can implement such a situation experimentally with arrays of dielectric resonators. The use of TE modes and the conducting plates below and above the resonators ensures at appropriate frequencies wellisolated resonances in the resonators and exponential decay outside. We have thus revealed a practical way to emulate Dirac-like equations for both massive and massless particles using classical wave systems.
While we presented an emulation for the Dirac oscillator in one and two dimensions, it is clear that our algebraic treatment allows other possibilities. The most promising example is the emulation of Dirac gyroscopes [29] as the number of states in this case is finite to begin with owing to the conservation of total angular momentum. This forces its realization on finite grids without approximations in that respect. The relativistic hydrogen atom, on the other hand, presents obvious difficulties regarding the steep potential needed near its singularity and the nearly flat potential at large distances requires separations too small for the radii of the resonators used. Some intermediate region of Rydberg states could possibly be emulated. Leaving the realm of integrable systems, it might be interesting to introduce random small perturbations in the positions of the wells. This might mimic some properties of random matrix Dirac operators [55]. The space-dependent wave functions φ n (x) are obtained by defining a specific form of the localized functions ξ A and ξ B . For definiteness, we choose these functions to be the ground state of a deep square well. The ground state φ 0 (x) and the density |φ 0 (x)| 2 are shown in figure A.1.