Bulk properties of honeycomb lattices of superconducting microwave resonators

We have realized different honeycomb lattices for microwave photons in the 4 to 8 GHz band using superconducting spiral resonators. Each lattice comprises a few hundred sites. Two designs have been studied, one leading to two bands touching at the Dirac points and one where a gap opens at the Dirac points. Using a scanning laser technique to image the eigenmodes of this new type of photonic lattices, we are able to reconstruct their band structure. The measured bands are in excellent agreement with ab initio models that combine numerical simulations of the electromagnetic properties of the spiral resonator and analytical calculations.

Superconducting photonic lattices for microwave photons hold the promise to simulate the behaviour of strongly interacting bosons in tailored 1D and 2D lattices [1][2][3]. In such systems, the excitation are microwave photons with a typical frequency around 6 GHz stored in superconducting resonators with high quality factor. The main interest of this platform is the possibility to reach the strongly interacting regime using Josephson junctions as a non-linear element. But, Josephson junctions are prone to disorder, and the realization of large disorder free non-linear lattices remains a very difficult task. Nevertheless, many-body effects have already been demonstrated in small uni-dimensional lattices [4][5][6][7][8].
Independently of controlling the interactions in the lattice, it seems also interesting to develop techniques allowing to probe lattices in the linear regime and to understand their properties [9,10]. This research direction, towards the realization of engineered band structures for photons, catches up with the rapidly growing field of photonic lattices and topological photonics [11][12][13]. In comparison to other photonic systems, the cryogenic environment and the low energy of the photons bring some technical difficulties to characterize the properties of such lattices.
In this manuscript, we use a laser scanning imaging technique in order to map the spatial distribution of the resonant modes across the lattice. This technique is similar to the laser scanning microscopy developed to observe the current density in superconductor films [14] and has already been used to characterize single a superconducting resonator [15,16] or coupled resonators [17]. It is an alternative to the scanning technique developed in [9]. From a Fourier analysis of the spatial distribution of the modes, we reconstruct the dispersion relation of the lattice bands. We apply this method to three honeycomb lattices with different designs, where the A and B sites in the unit cell correspond to identical or different resonators. In one case, we observe two bands touching at the Dirac point, while in the latter, a gap opens at the Dirac points. These properties of the bands are similar to the ones for electrons in graphene, even though the two systems are not described by the same wave equation. In order to understand our lattices and pinpoint the difference with electronic systems, we have developed two approaches to predict the properties of the lattice. The first approach consists in projecting the Maxwell equations on the basis of the resonator modes, which corresponds to a coupled mode theory (CMT), a technique widely used in photonics [18]. The band structure of the lattice is then obtained in terms of overlap integrals between the electric and magnetic fields of neighbouring resonators. The second approach relies on more intensive numerical simulations to obtain an equivalent circuit to the lattice, from which the band structure can be calculated. These two ab initio approaches allow us to reproduce the measured band structures with a very good accuracy.
The manuscript is organized as follows. In the first part, we present the properties of the spiral resonator that is used as the site of the different lattices. In the second part, we present the measurement of the band structure of the three lattice designs. And the last part details the two models that we have developed to predict the band and mode structures of the lattices.

I. SPIRAL RESONATOR PROPERTIES
The three different lattices studied in this manuscript are honeycomb lattices where each site consists of a superconducting spiral resonator as shown in figure 1. The spiral is made of a 4.3 µm wide Nb wire that is winded in an hexagonal pattern with an overall size of 300 µm. The lattices are patterned through photo-lithography and reactive ion etching starting from a ∼ 300 nm thick sputtered Nb layer on top of a silicon wafer. The design of the spiral was adjusted through numerical simulations such that the resonance frequency ω 0 of the fundamental mode is close to 2π × 6 GHz. This frequency can be finely tuned by adjusting the length of the wire at the center of the spiral. For the experiments shown in this manuscript, we use the two designs shown in figure 1a,b that were chosen to obtain two sites with slightly different ω 0 spaced by 2π × 120 MHz. The figure 1c shows the simulated self-admittance of the two designs, from which we deduce ω 0 as well as the resonance frequencies of the higher order modes. The second mode is expected to resonate around 14 GHz and plays no role in the experiments presented here, where all measurements are performed between 5 and 7 GHz.
When two spirals are approached at a distance d as (c) Self-admittance of the two spirals calculated for a port located at the end of the wire at the center of the spiral. The resonances correspond to the zero-crossings. The inset shows a zoom close to the fundamental resonance.
shown in figure 2, the fundamental modes couple, giving rise to two resonances at frequencies ω − and ω + . If the coupling is not too strong, the coupled modes can be expressed as a linear combination of the uncoupled modes. This approach that consists in using a reduced set of well chosen modes as a basis is known as the Coupled Mode Theory (CMT) [19][20][21]. Here, we use a basic CMT with only one mode per resonator. The Maxwell equations projected in this basis lead to a linear system, whose eigenvalues correspond to the resonance frequencies ω − and ω + . If the two resonators are identical with equal ω 0 , one obtains: where κ e (κ m ) is the electric (magnetic) coupling constant. These coupling constants are proportional to the overlap of the electric (magnetic) fields of the two resonators. We define where E i (H i ) is the electric (magnetic) field spatial dependence of the mode associated with the i = 1, 2 resonator. We assume these functions to be real. We explain in Appendix A how we compute these integrals from the charge and current distribution in the spiral that is obtained from the numerical simulation of a single resonator. In particular, one has to take into account charge and current images due to the dielectric interface at the surface of the sample and to the metallic ground plane below the sample. The coupling constants are then given by κ e = D 12 /D 11 and κ m = G 12 /G 11 . Figure 2a shows the evolution of κ e and κ m as a function of d for two spirals shown in fig. 1a. Because κ e > 0 and κ m < 0, the magnetic and the electric couplings add up to increase the mode splitting, while the mean frequency (ω + −ω − )/2 remains close to ω 0 . At large distances (d > 100 µm), the magnetic coupling dominates, while at short distances both couplings are important.
We compare the CMT predictions to a full numerical simulation of the two coupled resonators. From the simulation, we obtain the 2 × 2 admittance matrix Y (ω), which corresponds to the admittance matrix for two ports located at the ends of the two spirals (see figure). The two resonances ω + and ω − are then obtained as the zeros of det Y (ω). The CMT predictions coincide with the ones of the full numerical simulation at large distances (d > 20 µm) but fails at short distances as shown in figure 2b. This is because, at short distances, the coupling is too strong and the coupled modes are not linear superposition of the uncoupled modes. In the following, we show results for lattices where d = 5 µm or d = 30 µm. We therefore expect lattices with d = 30 µm to be well described by the CMT method, while lattices with d = 5 µm will require a more advanced model including the simulation of coupled sites. The half splitting (ω + − ω − )/2 obtained from the admittance matrix method is 2π × 200 MHz for d = 5 µm and 2π × 120 MHz for d = 30 µm. This number gives an estimate of the nearest neighbour coupling amplitude in a tight-binding description of the lattice.

A. Density of states
We have built honeycomb lattices consisting of a few hundred spiral resonators that fit on a 20 × 10 mm 2 sample, which is then mounted on the 1 K stage of a dry dilution fridge. We have characterized three designs of honeycomb lattice labeled G,SI and SII, as shown in figure 3. The G design corresponds to the situation where all sites are identical, as in graphene, and are occupied by the spiral shown in figure 1a. The other two designs, SI and SII, correspond to the more general situation where the two nonequivalent A and B sites of the honeycomb lattice are occupied by the two different spirals shown in figure 1a and 1b. Because, to first approximation, the sites only differ by their onsite frequency, the SI and SII designs realize a so-called Semenoff insulator [22]. In the graphene case, the expected band structure consists of two bands touching at the two Dirac points, while in the Semenoff case a gap opens at the Dirac points. In the two Semenoff designs, a horizontal boundary divides the sample in two halves: The A (B) sites in the lower half are occupied by the resonators that occupy the B (A) sites in the upper half. The two halves correspond to the same infinite lattice and therefore have the same bulk properties. The role of this boundary is to create valley Hall boundary states as discussed in [23].
We probe the transmission through the lattice using four microwave ports that are connected to four coplanar waveguides. Each waveguide is coupled to a single site located on the edge of the lattice (see figure 3). We first characterize the bulk properties of the lattice by estimating the density of states (DOS) that we deduce from transmission measurements. Figure 5 shows a typical transmission spectrum obtained through a lattice of type SII. The lattice modes are visible as sharp resonant peaks. Depending on the sample, we identify between 43% (G design) and 80% (SI design) of the expected resonances. This number is limited by our signal to noise ratio and by the finite width of the peaks. Peaks that are too weakly coupled to the measurement ports are not identified as resonances and peaks too close in frequency are counted as a single peak. By counting the number of peaks in a frequency window of 15 MHz for the design G and 15 MHz for the design SI and SII, we can however obtain an estimate of the DOS as shown in figure 5. As expected, we observe two bands for all three samples, with a clear gap for the SI and SII samples. We can also identify the two van Hove singularities corresponding to the two maxima in the DOS of each band. The two ab initio models to which we compare our data in the figure 5 only have a global offset frequency and are detailed in section ??. The frequencies of the remarkable points observed in the DOS (band minima, maxima and Dirac point) are well reproduced, but the measured DOS is not in quantitative agreement with the model predictions because a significant fraction of the modes is missing from the measured DOS.

B. Mode imaging and dispersion relation
In order to further characterize the lattice, we use a laser scanning technique to map the spatial variation of the modes identified in the transmission spectra. This imaging technique allows us to obtain a partial information on the dispersion relation of the lattice modes. The measurement consists in monitoring the transmission of one (or many) mode while scanning a laser spot across the lattice. The experimental setup is shown in figure 6. This method has been previously used to map the spatial profile of the resonant modes of a single resonator [15] or a chain of resonator. The optical setup was designed to obtain a laser waist on the sample of 60 µm, which is much larger than the size of the spiral wire width and spacing but smaller than the overall resonator size. This allows us to average the mode distribution over each site, while keeping sufficient resolution to resolve adjacent sites. We observe that the main effect of the laser is to induce dissipation on the illuminated site as observed in [15,24]. Because the modes are under-coupled to the probe ports, the laser induced increase of the loss results in a decrease of the mode transmission. To first order, the transmission drop for a given mode is proportional to the mode squared amplitude averaged over the illuminated area. In order to get rid of slow drifts and improve the signal to noise ratio, the laser intensity is modulated at a few kHz. We then digitally demodulate the transmitted signal measured with a VNA and record the amplitude of the in phase signal as a function of the laser position. Figure 6d shows a fine scan of one lattice site, which appears as a blurred hexagon with a central dip. This is consistent with the fact that the laser induced loss is maximal where the current density is large [24] as shown by the simulation in figure 5d. We attribute a single value for the mode intensity per site by averaging over a few measurement points well inside the hexagon surrounding one site. We have checked that the final result is rather insensitive to the details of the averaging procedure. In order to optimize the acquisition time, different scanning techniques have been tested, continuous (as shown in figure 5e) or raster. We obtain best results with a raster scan consisting of six measurement points per site, while we monitor the transmission change of all the peaks in a frequency window of about 1 GHz. This scan method allows us to image tens of modes in a single scan through the lattice. Figures 7 and 8 show the results of the mode imaging technique applied to some modes of G, SI and SII samples. In order to attribute a wavevector to each measured mode, we have developed a reconstruction technique to deduce the signed mode amplitude from the measured intensity. We do so by supposing that the mode is a linear combination of three known basis modes. As a basis, we choose the modes expected for a lattice having the same finite size geometry as the measured and described by a tight-binding model with only nearest neighbour coupling. We then look for the combination that matches best the measured mode intensity and use its sign to at-tribute a sign to the data. Details of the method are given in Appendix B. We then take the Fourier transform of the reconstructed signed mode amplitude. In the reciprocal space, several peaks lying on a circle are observed. This allows us to attribute a single value k corresponding to the norm of a wavevector to each mode. Figure 9 shows the dispersion relation ω(k) that we obtain through this analysis for the G and SI samples. The measured dispersion for the SII sample is almost identical to the one of the SI sample and is not shown here. We recover the two asymmetric bands observed in the DOS estimation. The mode dispersion is quadratic at small |k| leading to an effective mass for the photons in the lattice which is on the order of () for the lower (upper) band. At larger k, the dispersion relation clearly deviates from a quadratic behaviour. For the G sample, we expect a linear dispersion around the Dirac points but we are not able to image a sufficiently large number of modes in this region to clearly reproduce this behaviour. This is due to the finite size of the sample and also to the radial averaging over the direction of the wavevector. For the SI sample, we observe that bands curve again in the opposite direction close to the Dirac points.

III. LATTICE MODELS
In comparison to the DOS, the measured dispersion shown in is less affected by the fact that we are not able to probe all the modes. It allows us to precisely compare our data with two ab initio models. Following the analysis of the coupling of two spiral resonators, we extend the CMT and the admittance matrix model to the case of an infinite lattice. The two models have no free parameters and differ in the following way: the CMT model solely relies on the simulation of the charge and current distribution of a single spiral at resonance, while the admittance matrix model simulates the admittance matrix of a cluster of a few coupled sites. As explained in I, the CMT model gives a clear physical picture of the coupling in the lattice in terms of overlap integrals between the different sites but fails when the coupling is too strong. As a consequence, the band structure of the S lattice is well predicted by the CMT but observe a discrepancy for the G lattice, which is more pronounced for the lower band. The Y matrix calculation gives a more exact description of the lattice at any coupling but requires a more intensive numerical computation.

A. Coupled mode theory
We follow the derivation of [21] and adapt it to the specific case of coupled superconducting resonators. We only consider one mode per resonator and look for a solution to the Maxwell equations at a frequency ω as where E i (r) (H i (r)) is the electric (magnetic) field of the mode associated to the resonator at site i. The two fields verify where ω i is the resonance frequency of the mode of the resonator at site i. In addition to the overlap integrals defined in (3), we consider the integral where the surface S corresponds to all the metallic boundaries in the circuit, which consists of the resonators and the walls of the box enclosing the sample. This integral can be rewritten in terms of the overlap integrals D ij and G ij assuming that  where S i is the surface of the resonator at site i. The integral over S i is null because E i (r) is normal to the surface of the resonator. Using that the current density in the resonator at site j is given by J j (r) = n × H j (r), where n is the normal to the sample surface, we obtain Using this approximation and considering the volume integral of ∇ · (E i × H j ), one obtains In matrix form, the last equation writes  7). The frequency of the mode is indicated above each image.
where Ω is the diagonal matrix with elements Ω ii = ω i . In the same way, the volume integrals of ∇ · (H i ×E) and ∇ · (E i × H) lead to the following two equations Using the approximation (11) for M and eliminating b, we look for periodic solutions at frequency ω and obtain the following eigenvalue problem In the case of two resonators only, we recover the result given in (1). The matrices G, Ω and D are real symmetric matrices, but, in general, ΩG −1 ΩD is not symmetric. But its eigenvalues are real and positive with non-orthogonal eigenvectors. Equation (14) can be used to find the coupled mode frequencies of any ensemble of coupled resonators. In order to take advantage of the lattice periodicity, we relabel the a i and b i amplitudes in a given lattice cell as a µ (R) and b µ (R), where R is Bravais lattice vector identifying the position of the cell in the lattice and the µ index distinguishes the A and B sites. Using the same notation, we define the following overlap integrals between neigh-bouring sites in the lattice We now look for periodic solutions over the lattice as and define the Fourier transform of the D and G overlap matrices We suppose that the lattice contains N cells and we look for solutions with periodic boundary conditions such that R e jk · R = N δ k,0 , where k can take N values in the first Brillouin zone. Equations (12) and (13) then lead to ΩG −1 (k)ΩD(k)a(k) = ω 2 (k)a(k) (20) where Ω is now the diagonal matrix with the resonance frequencies of the A and B sites. The resulting eigenvalue FIG. 9. Band structure of the G and SI samples. For a given resonance frequency, the corresponding norm of the wavevector associated to the mode is obtained from the Fourier analysis of the measured spatial dependence as shown in figures 7 and 8. The measured dispersion relation is compared to the predictions of the two ab initio models detailed in section III (see also figure 5). The theoretical predictions are plotted as shaded areas, because, for a given |k|, the models predict different resonance frequencies depending on the orientation of the wavevector.
problem (20) gives the dispersion of the two bands as a function of k.
We obtain the band structure of our lattices by including all couplings up to the third neighbour coupling (see figure 4) in equations (18) and (19) and solving (20). The results are shown in figure . In order to compare the obtained band structure with the one of a tight-binding model with the same range of coupling, we consider the case of the G lattice where the A and B sites have the same resonance frequency. In order to obtain a simple analytical formula, we neglect the dependence of the overlap integrals with the direction of separation (e.g. we suppose that E A E B a1 = E A E B a ). With this approximation, a straightforward calculation shows that D(k) and G(k) are diagonal in the same basis, leading to the following dispersion relation ω ± (k) where the electric coupling constants are given by The magnetic couplings κ (i) m are given by the same expressions in terms of magnetic field overlaps and the link functions f i (k) correspond to f 1 (k) = 1 + e ik · a1 + e ik · a2 f 2 (k) = 2[ cos k · a 1 + cos k · a 2 + cos k · (a 1 − a 2 ) ] f 3 (k) = e ik · (a1+a2) + 2 cos k · (a 1 − a 2 ) Expression (21) can be compared to the tight-binding dispersion relation which comes from the diagonalisation of the tight-binding Hamiltonian The coupling t i corresponds to the i-th nearest neighbour coupling. The dispersion relation (21) and (22) match if one keeps only nearest neighbour terms, performs a first order expansion in κ (1) e and κ (1) m and identifies t 1 to ω 0 (κ e,1 −κ m,1 )/2. Including 2nd and 3rd nearest neighbour overlap integrals lead to a dispersion relation that cannot be identified to its tight-binding counterpart with the same coupling range.
As mentioned in the introduction, a strong motivation to build lattices of superconducting resonators is to include non-linear elements in the lattice. For example, in the case considered here, one could design spirals with an empty area in the center where a transmon might be inserted and coupled to the spiral. In this context, it can be interesting to obtain an effective tight-binding Hamiltonian that describes the linear behaviour of the lattice. Using the CMT approach, we obtained the two equations of motionȧ = Ωb (24) And the total energy in the system is given by In order to identify a and b with the two quadratures of the resonator modes, we normalize the basis fields to With this definition, and in the absence of coupling between the resonators (G = D = Ω), the equations of motion coincide with the Hamilton equations derived from H assuming standard commutation rules between the mode quadratures, [a m , a n ] = 0, [b m , b n ] = 0 and [a m , b n ] = i δ mn . In the following, we drop the factors and set = 1. When the overlap between adjacent sites are non-zero, the equations of motion do not coincide anymore with the naive Hamilton equations considering that the a's and b's commute. This comes from the fact that the basis used to project the Maxwell equations is not orthonormal, which modifies the commutation relations. This problem is well known in the calculation of electronic band structure using the LCAO method. In order to find an orthonormal basis, we adapt the Löwdin procedure used for electrons [25]. We first apply the canonical transformation a → Ω 1/2 a and b → Ω −1/2 b to remove the Ω dependence in the equation of motion, which transforms the Hamiltonian to withD = Ω 1/2 DΩ 1/2 andG = Ω −1/2 GΩ −1/2 . This step is not necessary if Ω is proportional to the identity. The equations of motion becomeȧ = b and Gḃ = −Da. We now apply the Löwdin transformation a →G −1/2 a and b →G −1/2 b in order to restore canonical commutation relations. The Hamiltonian becomes The equations of motion are modified toȧ = b anḋ b = −G −1/2DG−1/2 a, which now correspond to the Hamilton equations obtained with H assuming canonical commutation relations. We finally apply the reverse canonical transformation a → Ω −1/2 a and b → Ω 1/2 b to explicitly restore the Ω dependence. The final Hamiltonian is We can now introduce the bosonic operators c = (a + ib)/ √ 2 and rewrite H as with 2J = Ω −1/2G−1/2DG−1/2 Ω −1/2 − Ω. A rotating wave approximation can then be performed to transform the last term to c T (J/2)c † + hc.

B. Admittance matrix
In order to model more accurately the G lattice, we extend the admittance matrix method introduced in section I to an infinite lattice of resonators and we also take into account the possibility to have N ≥ 1 microwave ports per resonator. We introduce the voltage amplitude V (r ) that is a complex vector with 2N components corresponding to the voltage at the 2N ports describing the voltage in the lattice cell at position r . At a given frequency ω, the Kirchhoff's circuit laws describing the lattice can be written r Y r (ω)V (r + r) = 0 (32) where the sum over r is over the points of the Bravais lattice. The matrix Y 0 (ω) is the admittance matrix between the ports belonging to the same cell, while Y r (ω) is the admittance matrix between the ports corresponding to two cells separated by r, it has the property Y −r (ω) = Y T r (ω). We look for a periodic solution V (r) = V (k)e ikr and obtain r Y r (ω)e ikr V (k) = 0 (33) where V (k) is a 2N element vector. The dispersion relation is then obtained by solving det r Y r (ω)e ikr = 0 (34) The Y r (ω) matrix is obtained from the numerical simulation of a finite size lattice. The CMT calculations indicate that one should simulate a lattice with enough sites such that the central site is surrounded by all neighbouring sites up to the 3rd nearest neighbour. We therefore choose the geometry shown in figure with 16 sites. In order for the sum over r to converge rapidly when the distance between sites increases, the number of ports N and their location on the spiral must be well chosen. We observe from the simulations that N = 1 is not sufficient, while results obtained for N > 2 and various port locations give the same results. The results of the final simulations shown in are obtained with N = 3 ports located as shown in .
Compared to the CMT model, the admittance matrix is numerically more demanding because it requires the simulation of a much larger (here 16) number of sites. The dispersion obtained for the S lattice confirms the prediction of the CMT model, which were already in good agreement with our experimental data. For the G lattice, we observe a deviation compared to the CMT model leading to a much better agreement with the measured dispersion. We can thus conclude that our lattice are well under control and that ab initio electromagnetic simulation accurately describe their properties.