Circuit quantum electrodynamics of granular aluminum resonators

Granular aluminum (grAl) is a promising high kinetic inductance material for detectors, amplifiers, and qubits. Here we model the grAl structure, consisting of pure aluminum grains separated by thin aluminum oxide barriers, as a network of Josephson junctions, and we calculate the dispersion relation and nonlinearity (self-Kerr and cross-Kerr coefficients). To experimentally study the electrodynamics of grAl thin films, we measure microwave resonators with open-boundary conditions and test the theoretical predictions in two limits. For low frequencies, we use standard microwave reflection measurements in a low-loss environment. The measured low-frequency modes are in agreement with our dispersion relation model, and we observe self-Kerr coefficients within an order of magnitude from our calculation starting from the grAl microstructure. Using a high-frequency setup, we measure the plasma frequency of the film around 70 GHz, in agreement with the analytical prediction.


Current distribution
We perform all calculations in the limit of a one dimensional current distribution along the resonator. Here we present the results of a finite elements simulation of the current distribution for the five lowest stripline resonator modes. Supplementary Figure 1 shows the current distribution along (top row) and across (bottom row) the stripline resonator. The current along the resonator is at least two orders of magnitude higher than the current in the perpendicular direction. This validates the 1D current distribution assumption and allows us to exclude drum-like modes as a source of deviation from the linear dispersion relation. is shown in the bottom row. Jy is at least two orders of magnitude larger than Jx for all modes. Therefore, we conclude that up to the fifth mode the resonator follows the well known λ/2 current distribution, and that the measured nonlinear dispersion relation of the third mode, reported in the main text, is not caused by a drum mode like behavior.

Details on analytical model
In order to derive the equation of motion for a 1D JJ array we write the Kirchhoff laws and Josephson equations for two neighboring effective junctions as I n = I n+1 + C 0 dV n dt , where V n and χ n are the voltage and the phase on the n th node respectively, and I n is the current through the n th JJ. Combining these equations and introducing an excitation as an external current I ext cos(ωt) applied to the m th cell, where x m = l/2, we obtain where ϕ n = χ n+1 − χ n is the phase difference across the n th JJ. The JJ current is described as the following Substituting this expression to Supplementary Eq. (2) we obtain the equation of motion in the discrete limit In order to obtain the dispersion relation of the resonator we rewrite Supplementary Eq. (4) in the continuous limit We consider here the first resonance mode with sinusoidal current distribution I(x, t) = I(t) sin πx ; the corresponding phase difference is ϕ(x, t) = ϕ(t) sin πx . By substituting the phase difference ansatz in Supplementary Eq. (5), multiplying the equation by sin πx and integrating it along the resonator, we obtain the equation of motion of the resonator We approximate the Bessel function to first order J 1 [ϕ(t)] ∼ ϕ(t)/2, thus obtaining in the linear limit. Notably this equation is very similar to the motion equation of a current biased JJ. Solving Supplementary Eq. (7) we obtain the first resonance frequency of our system Performing the same calculations using the coordinate distribution of the higher resonance modes, we obtain the dispersion relation In order to derive the self-Kerr coefficient of the fundamental mode we solve the nonlinear equation of motion. The coupling of our resonator to the environment and internal losses of the resonator are introduced to Supplementary Eq. (6) as a damping term with a parameter γ = ω 1 /Q total (∼ 10 4 − 10 5 according to the experiment) whereĨ c = I c π 2 a 2 2 andC = C 0 + π 2 a 2 2 C J . Since we are working in resonance regime, we require the system to oscillate only with the driving frequency ω by assuming where ϕ a is the amplitude of the response for each JJ and δ is a phase delay due to losses in the system. Solving Supplementary Eq. (10) with ansatz (11) we obtain For small ϕ a we can expand the Bessel functions in series up to the third order, thus obtaining Since at resonance the response ϕ a reaches its highest value, we derive the resonance frequency of the nonlinear resonator by maximizing Supplementary Eq. (13).
One can see that in comparison to a single JJ with the resonance frequency, ω = ω 1 1 − ϕ 2 a 4 , the 1D array has similar, but lower first order nonlinearity. By relating the phase response to an average circulating photon numberN (see Appendix ), we obtain the self-Kerr coefficient for the fundamental mode In the following, we consider the cross-Kerr coupling between two different modes m and k with eigenfrequencies ω a and ω b . In order to obtain the cross-Kerr coefficients K mk one needs to solve the equation of motion with excitation terms aδ x − 2 (I m cos(ω a t) + I k cos(ω b t)), representing the drive of m th and k th modes We are mainly interested in the cross-Kerr coupling of the fundamental mode and we start with considering the coupling between the first and third modes. Similarly to the self-Kerr case, we look for a solution of the equation as a sum of the two driven modes Using the Jacobi-Anger identity we expand the nonlinear term in Supplementary Eq. (16) up to third order Here we limited the series only to the first and third modes, sin πx and sin 3πx , in which we are currently interested. We consider the case of strong pumping of the first mode and weak probing of the third mode. Therefore, Supplementary Eq. (16) can be simplified and splits into two separate equations for each mode Since the first mode drive is much stronger than the third mode drive, the equation where ϕ a is the amplitude of the third mode response for each JJ and δ is a phase delay due to losses in the system. Solving Supplementary Eq. (20) with ansatz (21) we obtain Again, we derive the resonance frequency by maximizing Supplementary Eq. (22) which gives the cross-Kerr coefficient between the first and third modes Performing the same procedure for first and all the other modes, we obtain the cross-Kerr coefficients

Circuit quantization
The total Q total and coupling quality factor Q c can be extracted from the measurement, allowing the average number of photons in the resonator to be calculated. In the case of one port waveguide for the first resonance mode it can be written asN where P in is the input power at the cavity port. The average number of photons relates to the amplitude of the current circulating in the resonator as where L J = 2eIc is the Josephson inductance of one junction 1 . At resonance, ϕ a reaches its maximal value, which is where Φ 0 = h/2e is the (superconducting) magnetic flux quantum. The self-Kerr nonlinearity is proportional to the second order of the phase response Substituting Supplementary Eq. (29) in Supplementary Eq. (14) we obtain

Details on samples
We measure the fundamental frequency of all resonators directly with the VNA. Due to its symmetry, the second mode is decoupled from the waveguide mode. For the two longest resonators made of grAl#1 film (4000 µΩ cm), the third mode, which is outside the frequency range of our VNA, can be excited by a second tone, generated by an RF generator, and detected via its cross-Kerr interaction with the first mode. Knowing the first and third resonance frequencies we use Eq. (3) of the main text to derive the plasma frequency which for f 1 = 6.287 ± 0.001 GHz and f 3 = 18.255 ± 0.001 GHz gives ω p = 68 ± 0.1 GHz. Figure 3a and b depict the measured S 11 response in the single photon regime for the stripline resonator of 0.6 mm length and 0.04 mm width made of grAl#1 film. A circle fit routine 2 is used to extract the internal and coupling quality factors Q i = 10 5 and Q c = 10 4 respectively 3 , and the resonance frequency f 1 ≈ 6.3 GHz. An overview of samples holders and sample geometries is provided in Supplementary Fig. 2 (in this appendix). The source radiation is two combined beams from two blackbodies at 77 K (liquid Nitrogen) and 300 K. The rotating polarizer P1 combines and polarizes the beam, which is then divided by two partial beams by the beam splitter BS. The wire grids of the beam splitter BS are oriented at 45 • to the normal of the drawing so that the polarization component perpendicular to the grid is transmitted and the component parallel to the grid is reflected. Two roof mirrors (one of each is movable) bring two components back to BS with a 90 • rotation of the polarization. If the two roof mirrors are equally spaced from BS, the input and output beams of the BS are polarized identically. After polarizer P2 only one orthogonal polarization is transmitted to the cryostat. The cryostat optical system consists of a lens at room temperature, and two aperture and lens pairs, at 4 K, and at 100 mK, in front of the sample.
We use the Martin-Puplett interferometer (MPI) as a broad-band illumination source, with a resolution up to 1 GHz. A schematic drawing of the MPI is shown in Supplementary Fig. 3. Three wire grids are used as polarizer P1, beam splitter BS and polarizer P2. If the polarization of an incident wave is parallel to the wires, the wire grid behaves like the surface of a metal and the wave is reflected, whereas it would be a perfectly transparent element for a wave that is polarized orthogonally to the wires. The source radiation, emitted from two black bodies, at room temperature (red) and liquid nitrogen temperature (blue), is combined on the polarizer P1, which is a wire grid inclined by 45 • with respect to the plane of the drawing. The polarizer P1 is rotating with frequency ω 0 about its axis, which creates the output beam containing two orthogonal polarizations that are swept over all possible orientations by a full rotation of P1. In the following discussion we'll consider P1 at a fixed moment in time. The beam is divided into two partial beams by the beam splitter BS, which is a fixed wire grid at a 45 • angle to the incoming beam. The component of the incident beam with polarization parallel to the wires of BS is transmitted towards the fixed roof mirror, and the orthogonal component is reflected towards the movable roof mirror. The roof mirrors reflect the incident beams and flip their polarizations by 90 • . The movable roof mirror can be displaced by an amount dx. The two beams recombine at BS with an accumulated path difference ∆ = 2dx, resulting in a phase difference 2π∆/λ which gives rise to interference. After the reflection both polarizations are flipped with respect to their first encounter with BS, therefore the transmission/reflection routine is inverted. The polarizer P2 provides a reference for the ω 0 modulation of the polarization produced by P1. The single polarization output beam enters the cryostat through the room temperature lens. Here, the incoming beam is collected and sent to the cold optics, consisting on two additional focusing lenses, an in-focus aperture to reduce the intensity, and an out-of-focus aperture to crop the image. The intensity of the on-sample radiation for a single wavelength is 4 where I 0 is defined among others by the diameter of the apertures and position of the sample with respect to the optical exes of cryostat optical system. I(∆) is an even, periodic function of the roof mirror displacement, with the first maximum appearing at the origin, which is result of the alternatingly constructive and destructive interference. The interferogram is the modulated term of Supplementary Eq. (32) integrated over all wavelengths. The parity of the integrand allows us to recast the cosine modulation as an exponential having the same argument, thus showing that the interferogram and the spectrum are Fourier transform (FT) pairs. The interferogram has a global maximum at the origin since different wavelengths will interfere in a fully constructive fashion only in the case of zero path difference. Furthermore, the interferogram is a stronger signal than the monochromatic intensity, since it encodes contributions from all wavelengths. Interferograms and their corresponding spectra (MPI responce versus illumination frequency) are shown in Supplementary Fig. 4.
The interferogram is generated by recording the on-sample irradiation at roof mirror steps T x ∼ 10 µm. This is equivalent to multiplying the a priori continous signal with a comb of Dirac deltas with spacing T x . The FT of this product is a convolution of the spectrum with a 1/T x Dirac comb in impulse space, i.e. an array composed by images of the spectrum spaced 1/T x apart. The images are symmetrical and bounded by some ±f max . The highest frequency that can be attained before image overlapping and subsequent aliasing is given by step T x , for T x ≈ 50 µm the highest frequency f max = c/2T x ≈ 3000 GHz. For all measurements presented in Supplementary Fig. 4 lowpass filters are used, for Al samples the cutoff frequency is 180 GHz, for grAl#2 and grAl#3 the cutoff frequency is 300 GHz

Critical current
The geometry and the switching current measurements for a DC-SQUID made with granular aluminum are presented in Fig 5. The measurements over 4 samples with resestivities 250, 1520, 3200 and 5550 µΩ·cm are summerized in Supplementary Fig. 6.