<figure><img src="image_1.png" /></figure>

HAL
open science

# Dynamical Coulomb blockade: An all electrical probe of the ultrastrong light-matter coupling regime

Usama Iqbal, Yanko Todorov, Christophe Mora

► To cite this version:

Usama Iqbal, Yanko Todorov, Christophe Mora. Dynamical Coulomb blockade: An all electrical probe of the ultrastrong light-matter coupling regime. Physical Review Research, 2024, 6 (3), pp.033097. 10.1103/PhysRevResearch.6.033097. [https://arxiv.org/abs/2009.03097] [https://arxiv.org/pdf/2009.03097.pdf] [https://arxiv.org/abs/2009.03097]

HAL Id: hal-04696603

https://hal.science/hal-04696603v1

Submitted on 16 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamical Coulomb blockade: An all electrical probe of the ultrastrong light-matter coupling regime

Usama Iqbal and Yanko Todorov*

Laboratoire de Physique de l'École Normale Supérieure, ENS, Université, PSL, CNRS, Sorbonne Université,
Université de Paris Cité, F-75005 Paris, France

Christophe Mora

Université Paris Cité, CNRS, Laboratoire Matériaux et Phénomènes Quantiques, F-75013 Paris, France

<figure><img src="image_3.png" /></figure>

(Received 12 April 2024; accepted 3 July 2024; published 23 July 2024)

The regime of ultrastrong light-matter coupling (USC) has attracted considerable interest recently, as it predicts fascinating quantum phenomena such as the presence of virtual photons in the ground state. Here we propose theoretically an approach to probe the quantum correlations in USC systems based on an electrical transport experiment. Namely, we consider the case in which the light-matter coupled system constitutes the electromagnetic environment of an ultrasmall tunnel junction operating in the regime of a dynamical Coulomb blockade (DCB). Unlike optical spectroscopy, such a probe can be seen as a nondemolition measurement of the USC system ground state. We have developed two approaches in order to evaluate the effect of the light-matter coupled states on the tunneling process. The first is a generalized, multimode Peierls substitution for the tunneling Hamiltonian. The second is based on the P(E) theory describing a DCB, where the USC system is treated as an equivalent circuit. The two approaches are consistent, whereas the second one enables us explicitly to take into account the presence of dissipation. This study has been performed for the zero-temperature case, which is pertinent for the terahertz spectral range where excitation energies are typically much higher compared to the microwave domain, where DCB was first studied.

DOI: 10.1103/PhysRevResearch.6.033097

## I. INTRODUCTION

Since the seminal work of Purcell [1], it has been well known that light-matter interaction can be modified substantially by coupling the quantum transition with a single-mode electromagnetic resonator. When the light-matter coupling energy $\hbar\Omega_R$ exceeds the dissipation rates of the system, energy is periodically exchanged with the resonator at a frequency $\Omega_R$, also known as the vacuum Rabi frequency [2]. The system enters the strong-coupling regime where the cavity mode is split into two light-matter coupled (polariton) states separated by an energy $2\hbar\Omega_R$ [2]. The past decade has seen the emergence of an even stronger interaction regime, where the coupling constant $\Omega_R$ becomes comparable to the frequency of the matter excitation $\omega_m$. This regime with $\Omega_R/\omega_m \approx 1$ is known as ultrastrong light-matter coupling (USC) [3], and it establishes new frontiers for cavity quantum electrodynamics [4,5]. Indeed, coupling with light can be so strong that it leads theoretically to important changes of the material properties, such as the emergence of new phases [6–8] and the

possibility to trigger superconductivity [9,10]. The USC regime also contains fascinating quantum phenomena. As the coupling constant $\Omega_R$ grows, the polariton states become squeezed radiation states [3,11,12]. The ground state of the system is also completely changed, and it acquires a population of virtual photons [3]. Under certain conditions, these can be extracted from the system [3,13,14]. The effect of converting such virtual photons into real photons is known as the dynamical Casimir effect [15], which is analogous to Hawking radiation from black holes [16].

Up to now, USC has been observed in many spectroscopic experiments [4,5], however such experiments do not provide direct access to the quantum optical properties of the light-matter coupled states. Indeed, in quantum optics, one can access the nonclassical nature of radiation states by measuring high-order correlation functions [17], which requires designing polariton emitters and the exploitation of very sensitive detectors. A very appealing alternative, inspired from mesoscopic physics, is to access the quantum fluctuations of “photon” states in electronic transport experiments [18]. Here “photon” has been set in quotes, as the majority of mesoscopic devices employ transmission lines and inductor-capacitor (LC) resonators rather than optical microcavities [19]. Recently, such quantum devices featuring USC and electronic transport have been studied in the terahertz range [20–23]. In addition to fundamental theoretical studies of light-matter coupled states [24–27], the electrical transport in semiconductor and organic devices is currently an actively

<footnote>*Contact author: yanko.todorov@phys.ens.fr</footnote>
investigated experimental topic [28–31]. The USC regime between superconducting qubits and transmission lines has also been studied [4].

A basic system that allows us to probe the quantum fluctuations of an electromagnetic resonator is an ultrasmall tunnel junction with a capacitance $C_T$, such as the charging energy of a single electron $e^2/2C_T$, is comparable to the oscillation quantum $\hbar\omega_{LC}$ [32–34]. When the ratio $\rho_0 = (e^2/2C_T)/\hbar\omega_{LC}$ is sufficiently large, the tunneling of a single electron is conditioned by the quantum fluctuations of the electromagnetic environment, leading to the suppression of low-voltage conductance. This effect is known as a dynamical Coulomb blockade (DCB). Recently, it was shown theoretically that DCB can be used to probe nonclassical states of the electromagnetic resonator [35]. On the experimental side, DCB has been observed in the GHz range with quantum circuits [36–41]. At higher energies, DCB has been observed with graphene phonons in the terahertz (THz) range [42], and in scanning tunneling microscopy (STM) experiments for higher frequencies [43,44].

Here, we investigate a THz quantum-metamaterial architecture that takes advantage of the DCB in order to probe the quantum-optical properties of the polariton quasiparticles. This can be achieved by using a THz metamaterial resonator with highly sub-wavelength capacitive parts that interacts with highly doped quantum wells in the USC regime. The lumped element nature of the THz meta-atom allows us to implement a tunnel junction in one of the capacitors, which experiences the polaritonic environment in the DCB regime. In our metamaterial architecture, there is an intriguing link between the DCB phenomenon and the USC regime with a single electron in the quantum wells. We show that transport measurements through the junction allow one to access the squeezing of the light-matter states that appear in the USC [3,11], in realistic systems, opening the possibility to perform quantum optical experiments “on chip” with THz photons.

This paper is organized as follows. In Sec. II we present the system without dissipation where we recall the main quantum-optical features of the USC. We derive an equivalent impedance for the polaritonic environment of the junction, enabling the polaritonic system to be treated on the same footing as mesoscopic circuits where DCB has been observed so far. We discuss a derivation of the DCB in that system based on a generalized Peierls substitution. Since there is no absorption in the system, and the tunnel junction is treated perturbatively, we can see the DCB as a nondemolition probe of the polariton states. In Sec. III we present a more realistic model, including dissipation, within the framework of the circuit $P(E)$ theory, where the DCB effect is evaluated from the real part of the impedance seen by the tunnel junction [34]. We show that the Peierls substitution approach yields the same results as the impedance approach of the $P(E)$ theory. In Sec. IV we discuss our results and the possibility to measure quantum-statistical properties for realistic devices operating in the THz range.

Our studies are performed here for the zero-temperature case, which is justified for the THz frequency range where excitation energies are typically much higher compared to the microwave domain where DCB was first studied. Indeed, the frequency of 1 THz corresponds to a temperature of

50 K, which means that in principle the zero-temperature approximation can be valid for experiments performed with liquid helium. Furthermore, as shown in Ref. [34], the $P(E)$ theory can be generalized for a finite temperature by numerical methods.

## II. HAMILTONIAN MODEL

### A. Main features of USC

In the THz range, the USC regime between the electronic transitions in highly doped quantum wells and electromagnetic resonators has been explored both with microcavity resonators [45] [Fig. 1(a)] and circuitlike architectures [46–48] [Fig. 1(b)]. Both cases can be described by the same Hopfield-like Hamiltonian, which allows us to study the features of USC, as further recalled. Here, we will focus on the latter case, as it naturally allows implementing the DCB. The first ingredient of the coupled system is an LC resonator of frequency $\omega_{LC}$ described by a quantum Hamiltonian $\hat{H}_{LC} = \hat{Q}^2/2C + \hat{\Phi}^2/2L$, where $\hat{Q}$ and $\hat{\Phi}$ are the charge and flux operators which satisfy the commutation relation $[\hat{Q}, \hat{\Phi}] = i\hbar$. We consider for now the simple geometry with a single capacitor schematized in Fig. 1(b), but our results remain valid for more complex circuits comprising multiple capacitive parts as described further. The charge and flux operators can then be expressed from the ladder operators $a$ and $a^\dagger$ as

$$ \hat{Q} = \sqrt{\frac{C\hbar\omega_{LC}}{2}}(a + a^{\dagger}), \quad \hat{\Phi} = \sqrt{\frac{L\hbar\omega_{LC}}{2}}i(a - a^{\dagger}) \quad (1) $$

such that the circuit Hamiltonian becomes $\hat{H}_{LC} = \hbar\omega_{LC}(a^\dagger a + 1/2)$.

The electronic excitation that we consider arises from the transition between two confined subbands 1 and 2 in a two-dimensional electron gas hosted in the quantum wells depicted in Figs. 1(a) and 1(b). The wells are highly doped in order to ensure a high density of dipoles [45,49]. The Fermi energy $E_F$ is then located close to the edge of the second subband. Electronic interactions give rise to a collective electronic excitation, the intersubband plasmon [45,49] with a frequency $\tilde{\omega}_{21} = \sqrt{\omega_{21}^2 + \omega_p^2}$, where $\omega_{21}$ is the $1 \to 2$ transition frequency, $\omega_p^2 = N_e e^2 / (m^* \epsilon \epsilon_0 S L_{eff})$ is the plasma frequency of the electron gas with $N_e$ the total number of electrons per quantum well, $m^*$ is the effective electron mass in the semiconductor material constituting the quantum wells, $\epsilon$ and $\epsilon_0$ are, respectively, the material and vacuum electric permittivity, $L_{eff}$ is an effective thickness of the quantum well [49], and $S$ is its surface. We will consider the linear regime of light-matter interaction, where we can associate a bosonic Hamiltonian with the plasmon excitation $\hat{H}_{pl} = \hbar\tilde{\omega}_{21}(p^\dagger p + 1/2)$ [49].

The electronic polarization of the quantum wells interacts with the quasistatic electric field of the capacitor [50], which is proportional to the charge operator from Eq. (1). We express the interaction in the framework of the Power-Zienau-Wooley (PZW) representation of quantum electrodynamics [51], which is a full representation of the light-matter interacting Hamiltonian that contains all features of the USC. The
<figure><img src="image_4.png" /><figcaption>FIG. 1. Two alternative arrangements for probing ultrastrong light-matter coupling: the quantum wells are embedded in a microwave cavity (a) [45], or coupled to a lumped element LC circuit in (b) [46]. For both devices, we represent the following: (c) Frequencies of the normal modes $\omega_{\pm}$ of the coupled system as a function of the frequency of the electromagnetic resonator $\omega_{LC}$. The other parameters for this plot are $\omega_{P} = 1.2$ THz, $\omega_{21} = 2$ THz ($\tilde{\omega}_{21} = 2.33$ THz), and $F = 0.7$. The gray area is the polariton gap. (d) Normal modes as a function of the coupling parameter $\omega_{P}$ for $\omega_{LC} = \omega_{21} = 2$ THz and $F = 0.7$. (e) Ground-state photon population $\langle G|a^{\dagger}a|G \rangle$ as a function of $\omega_{P}$ using the same parameters as in (d). (f) Quantum fluctuation enhancement factors for charge, Eq. (8) (red curve), and flux, Eq. (9) (blue curve). The dashed line is the product of the two curves. Note that panel (f) pertains only for the case of a circuit [panel (b)], whereas for the microcavity case (a), the two field quadratures should be considered instead of charge and flux.</figcaption></figure>

complete Hamiltonian of the system now becomes [50]

$$ \hat{H} = \hbar\tilde{\omega}_{21}(p^{\dagger}p + 1/2) + \hbar\omega_{LC}(a^{\dagger}a + 1/2) + \frac{\hbar\omega_{P}}{2}\sqrt{F\frac{\omega_{LC}}{\tilde{\omega}_{21}}}(p + p^{\dagger})(a + a^{\dagger}). \quad (2) $$

This is a Hopfield-like Hamiltonian [52] which describes two coupled harmonic oscillators. The coefficient $F$ introduced above is an overlap coefficient that satisfies $0 < F < 1$. It quantifies the geometric overlap between the electronic polarization and the electric field of the capacitor $C$; it eventually takes into account the presence of other capacitive elements in the circuit. For the case of a circuit with a single capacitor,

we can write $F = f_{21}N_{qw}L_{eff}/d$, where $N_{qw}$ is the number of quantum wells, $d$ is the thickness of the capacitor [50], and $f_{21} = 2m^*\omega_{21}z_{21}^2/\hbar$ is the oscillator strength of the QW transition, with $z_{21}$ the length of the microscopic dipole, which is determined by the envelope wave functions of the heterostructure potential [53]. In this expression, the capacitor is treated as a perfect planar capacitor where fringing field effects are neglected. Note that in our representation, the light-matter term is proportional to the plasma frequency $\omega_{P}$, which will be further referred to as a “coupling parameter.”

The Hamiltonian from Eq. (2) is solved by employing the Hopfield-Bogoliubov procedure [45,52], which converts it into a sum of two independent harmonic-oscillator
contributions with normal mode frequencies $\omega_+$ and $\omega_-$:

$$ \hat{H} = \hbar\omega_{+}(\Pi_{+}^{\dagger}\Pi_{+} + 1/2) + \hbar\omega_{-}(\Pi_{-}^{\dagger}\Pi_{-} + 1/2). \quad (3) $$

The normal mode frequencies satisfy the bi-quadratic equation [49]

$$ (\omega^2 - \omega_{LC}^2)(\omega^2 - \tilde{\omega}_{21}^2) = F\omega_p^2\omega_{LC}^2, \quad (4) $$

and the new polariton operators $\Pi_{\pm}$ are expressed as linear combinations from the old ones:

$$ \Pi_{\pm} = x_{\pm}a + y_{\pm}a^{\dagger} + m_{\pm}p + t_{\pm}p^{\dagger}. \quad (5) $$

Note that we have adopted a slightly different form of the light-matter interaction Hamiltonian where the ladder operators for the circuit and matter excitation are perfectly symmetric. As a result, the Hopfield coefficient in the decomposition Eq. (5) are all real numbers. As shown in Appendix A, we can further express them analytically as a function of the frequencies $\omega_{\pm}$, $\omega_{LC}$, $\tilde{\omega}_{21}$. The inverse Hopfield transformation can also be provided explicitly, which is very useful for the subsequent results.

In Figs. 1(c) and 1(d) we recall the well-known behavior of the normal mode frequencies $\omega_{\pm}$. In Fig. 1(c) the frequency of the material excitation is fixed, with $\omega_p = 1.2$ THz, $\omega_{21} = 2$ THz ($\tilde{\omega}_{21} = 2.33$ THz), and $F = 0.7$, while the frequency of the LC circuit is varied. We observe the typical anticrossing behavior of the two polariton modes, with an opening of a gap that arises from the antiresonant terms of the light-matter Hamiltonian Eq. (2) [45,49]. The gap, indicated by the gray area in Fig. 1, is comprised between the frequencies $\tilde{\omega}_{21}$ and $\tilde{\omega}_{21} = \sqrt{\omega_{21}^2 + (1-F)\omega_p^2}$ [45,49].

In Fig. 1(d) we set $\omega_{LC} = \omega_{21} = 2$ THz and we vary the coupling parameter $\omega_p$. We recover the nonlinear behavior of the polariton branches $\omega_{\pm}$ as a function of $\omega_p$, which is due to the progressive detuning of the matter frequency $\tilde{\omega}_{21}$ from the LC resonance. We recall that the renormalization of the matter excitation by plasmonic effects arises from the quadratic term of the PZW Hamiltonian [45].

We now review the main features of the USC regime in the spirit of the original work by Ciuti et al. [3]. First of all, let us note that the ground state of the Hamiltonian Eq. (3) is the tensor product $|G\rangle = |G\rangle_+|G\rangle_-$, where $|G\rangle_{\pm}$ are the ground states of each normal mode that satisfy $\Pi_{\pm}|G\rangle_{\pm} = 0$. As pointed out in Ref. [3], this ground state cannot be expressed as a tensor product of the photon and electron vacua, $|0\rangle_{LC}|0\rangle_{mat}$, for the uncoupled systems, as would be the case if the rotating-wave approximation were applied to the Eq. (2). Instead, the ground state of the coupled system contains a finite number of circuit or matter excitations. Indeed, by using the inverse Hopfield transform, we can show that

$$ \langle G | a^{\dagger} a | G \rangle = y_{+}^{2} + y_{-}^{2} = \gamma_{+} \frac{(\omega_{+} - \omega_{LC})^{2}}{4\omega_{+}\omega_{LC}} + \gamma_{-} \frac{(\omega_{-} - \omega_{LC})^{2}}{4\omega_{-}\omega_{LC}}. \quad (6) $$

Here we introduced the notations

$$ \gamma_{+} = \frac{\omega_{+}^{2} - \tilde{\omega}_{21}^{2}}{\omega_{+}^{2} - \omega_{-}^{2}}, \quad \gamma_{-} = \frac{\tilde{\omega}_{21}^{2} - \omega_{-}^{2}}{\omega_{+}^{2} - \omega_{-}^{2}}, \quad (7) $$

which will be systematically used in the following. Equation (6) expresses the ground-state population of circuit photons through the antiresonant Hopfield coefficients $y_{\pm}$, as explained in Ref. [3]. In Fig. 1(e) we plot $\langle G|a^{\dagger}a|G\rangle$ as a function of the coupling constant $\omega_p$ using the same parameters as those for Fig. 1(d). The ground-state population of photons increases quadratically as $\approx \omega_p^2/\omega_{LC}^2$, and it can reach values as high as 0.02 in the USC, where the coupling constant is an important fraction (70%) of the bare resonator frequency.

In the context of quantum circuits, it is also interesting to examine the value of the charge and flux quantum fluctuations in the ground state. By using Eq. (1) as well as the inverse Hopfield transform, we arrive at the following results:

$$ \langle G | \hat{Q}^2 | G \rangle = \langle 0 | \hat{Q}^2 | 0 \rangle \left( \frac{\gamma_{+}}{\omega_{+}} + \frac{\gamma_{-}}{\omega_{-}} \right) \omega_{LC}, \quad (8) $$

$$ \langle G | \hat{\Phi}^2 | G \rangle = \langle 0 | \hat{\Phi}^2 | 0 \rangle \frac{\gamma_{+} \omega_{+} + \gamma_{-} \omega_{-}}{\omega_{LC}}. \quad (9) $$

In Fig. 1(f) we plot the values of $\langle G|\hat{Q}^2|G\rangle^{1/2}/\langle 0|\hat{Q}^2|0\rangle^{1/2}$ (red curve) and $\langle G|\hat{\Phi}^2|G\rangle^{1/2}/\langle 0|\hat{\Phi}^2|0\rangle^{1/2}$ (blue curve) as a function of $\omega_p$. Here $|0\rangle = |0\rangle_{LC}$ is ground state of the uncoupled LC circuit. These quantities represent the enhancement of quantum fluctuations of charge and flux in the new ground state of the system. As the coupling constant $\omega_p$ increases, we observe a strong antisqueezing of the charge fluctuations while the flux fluctuations are squeezed. These tendencies are correlated to the increase of the ground-state photon occupancy [Fig. 1(e)]. The black dotted line is the product of the two enhancement factors; since it is larger than unity, we have $\langle G|\hat{Q}^2|G\rangle^{1/2}\langle G|\hat{\Phi}^2|G\rangle^{1/2} > \langle 0|\hat{Q}^2|0\rangle^{1/2}\langle 0|\hat{\Phi}^2|0\rangle^{1/2} = \hbar/2$ and the Heisenberg inequality is well-satisfied. We also note that the relative importance of the charge antisqueezing is almost twice that of the virtual photon occupancy of the USC ground state. The antisqueezing of charge fluctuations can thus be considered as a more tangible manifestation of the USC regime with circuit resonators.

For the case of photonic resonators, similar results to those described in Fig. 1(f) hold for the two quadratures of the microcavity electric field. For the circuit case, however, charge and flux can be coupled with the electronic transport in mesoscopic circuits, which is an intriguing possibility as described further.

## B. Link between USC and DCB

The minimum splitting $\Omega$ between the polariton states in Fig. 1(a) can be expressed as $\Omega = \omega_p\sqrt{F}$ [45,49]. The ratio $\Omega/\tilde{\omega}_{21}$ is often used as a figure of merit, with the onset of the USC considered to occur for $\Omega/\tilde{\omega}_{21} \gtrsim 0.1$ [3]. It is interesting to examine this figure of merit in the single-electron limit. Theoretically, this limit is achieved in a circuit with a very small capacitor [50] and a single quantum well $N_{qw} = 1$. In this limit, collective electronic effects can be discarded and $\tilde{\omega}_{21} \to \omega_{21}$. We will consider instead the square of the expression $(\Omega/\omega_{21})^2$, which becomes

$$ \frac{\Omega^2}{\omega_{21}^2} (N_e = 1) = F \frac{e^2}{m^* \omega_{21}^2 \varepsilon_0 S L_{\text{eff}}} \quad (10) $$
Since the overlap factor is $F = f_{21}L_{\text{eff}}/d$, Eq. (10) can be rewritten as

$$ \frac{\Omega^2}{\omega_{21}^2} (N_e = 1) = \frac{2e^2 z_{21}^2}{\hbar \omega_{21} \varepsilon \varepsilon_0 S d} = \left(2 \frac{z_{21}}{d}\right)^2 \frac{e^2}{\hbar \omega_{21} 2C}. \quad (11) $$

We used here the geometrical capacitance of the circuit $C = \varepsilon \varepsilon_0 S/d$. The ratio in the parenthesis, $2z_{21}/d$, must be typically on the order of unity in order to achieve good spatial overlap between the electronic polarization and the electric field of the capacitor. The quantity $e^2/2C$ is the Coulomb energy of a single electron deposited on the capacitor plates, whereas $\hbar\omega_{21}$ is the oscillation quantum of the circuit. The ratio $e^2/2C\hbar\omega_{21}$ is precisely the figure of merit introduced by Devoret et al. to quantify the feasibility of a dynamical Coulomb blockade in a single-mode circuit [32], as recalled further: the latter occurs if this ratio is close to unity. Thus, Eq. (11) shows that there is an intriguing link between the USC in a circuit and the DCB. Clearly, mesoscopic systems that can feature strong DCB are also very suitable to observe USC at the single-electron limit [21].

Let us recall the basics of the DCB phenomenon. To this end, we consider an ultrasmall tunnel junction with a capacitance $C_T$ connected to an inductor $L$. This arrangement forms an LC circuit, with a total Hamiltonian $\hat{H}_{\text{LC}} + \hat{H}_{\text{tr}} = \hbar\omega_{\text{LC}}(a^\dagger a + 1/2) + \hat{H}_{\text{tr}}$, where we introduced the Hamiltonian describing the electronic transport through the tunnel junction [34,54]:

$$ \hat{H}_{\text{tr}} = \sum_{l} E_{l} c_{l}^{\dagger} c_{l} + \sum_{r} E_{r} c_{r}^{\dagger} c_{r} + \hat{T} + \hat{T}^{\dagger}, \quad (12) $$

$$ \hat{T} = \sum_{l,r} \tau_{l,r} c_l^\dagger c_r e^{i e \hat{\Phi} / \hbar}. \quad (13) $$

The first two terms of the right-hand side of Eq. (12) describe the quasiparticle energies of the left/right lead of the junction. The extra term $e^{i e \hat{\Phi} / \hbar}$ is known as a Peierls phase factor [55]. Thanks to its presence, the tunnel coupling operator $\hat{T}$ effectively achieves the simultaneous transfer of an electron across the junction while displacing the influence charge of the capacitor from $\hat{Q}$ to $\hat{Q} - e$ [54]. Assuming zero temperature, and following the standard theory of DCB [32,34], we need to compute the autocorrelation function for the flux operator:

$$ \langle G | \hat{\Phi}(t) \hat{\Phi}(0) | G \rangle = \frac{L}{2} \hbar \omega_{\text{LC}} e^{-i \omega_{\text{LC}} t}. \quad (14) $$

We proceed by introducing the function $J(t) = (e^2/\hbar^2)\langle G | [\hat{\Phi}(t) - \hat{\Phi}(0)] \hat{\Phi}(0) | G \rangle$ and thus define the $P(E)$ function,

$$ P(E) = \int_{-\infty}^{\infty} \frac{dt}{2\pi\hbar} e^{J(t)} e^{iEt/\hbar}, \quad (15) $$

central in the theory of DCB. The function $P(E)$ represents the probability density for the environment to absorb the energy $E$ during a tunneling event. Using the above expressions, this function can be determined in a closed form:

$$ P(E) = e^{-\rho_0} \sum_{n}^{\infty} \frac{\rho_0^n}{n!} \delta(E - n\hbar\omega_{\text{LC}}). \quad (16) $$

<figure><img src="image_5.png" /><figcaption>FIG. 2. A terahertz meta-atom resonator coupled both with highly doped quantum wells and a tunnel junction. The current-voltage characteristics of the system is probed between the leads a and b.</figcaption></figure>

Respectively, the zero-temperature differential conductance is derived by integrating the $P(E)$ function:

$$ G(U)/G_K = e^{-\rho_0} \sum_{n}^{\infty} \frac{\rho_0^n}{n!} \Theta(U - n\hbar\omega_{\text{LC}}/e), \quad (17) $$

with $\Theta(x)$ being the Heaviside function, and $U$ is the bias applied on the junction. Note that we can directly compare energies with values of the bias by taking all energy scales in eV or meV. The following quantity has been defined above:

$$ \rho_0 = \frac{Z_{LC}}{\pi R_K} \quad (18) $$

with $Z_{LC} = \sqrt{L/C_T}$ the characteristic impedance of the circuit, and $R_K = 1/G_K = \hbar/e^2$ the resistance quantum. We can also reexpress the quantity $\rho_0$ as

$$ \rho_0 = \frac{e^2/2C_T}{\hbar\omega_{\text{LC}}}. \quad (19) $$

DCB can be observed when $\rho_0$ approaches unity, indicating that the electrostatic charging energy of a single electron must be comparable to the oscillation quantum of the circuit. Alternatively, we can state that the typical size of the quantum fluctuations of the circuit, $\langle 0|\hat{Q}^2|0 \rangle^{1/2} = (C_T \hbar \omega_{\text{LC}}/2)^{1/2}$, must be comparable to the value of the elementary charge $e$. Remarkably, precisely the same quantity appears for the USC figure of merit with a single electron, Eq. (11).

The DCB phenomenon is thus related to the quantum fluctuations of the ground state of the circuit. From our discussion in the previous paragraph, it is clear that the modifications of the charge and flux quantum fluctuations are yet another feature of the USC regime. Guided by these considerations, we now consider a more complex device where a tunnel junction is coupled to a polaritonic system. The resulting quantum metamaterial resonator is described in Fig. 2, and it is inspired from a three-dimensional architecture that has been recently introduced to study the USC regime with THz intersubband excitations [46] as well as for quantum-well detectors [56]. We shall not discuss here the possible fabrication steps that can lead to the realization of the device concept from Fig. 2, but we will explore it instead from a theoretical point of view.
As shown in Fig. 2, the device proposed can be described as a high-frequency circuit made from two double-metal capacitors, C and $C_T$, which are connected with an inductive element to form an equivalent LC resonator. In this configuration, only the capacitor C is filled with highly doped quantum wells, while a tunnel junction is inserted in the capacitor $C_T$. The circuit, also called a “meta-atom” in the context of THz spectroscopy [46], is supplied with electrical leads (see also Ref. [56]) that allow us to probe the current-voltage characteristics of the junction. The light-matter coupled system is thus still described by the Hamiltonian [Eq. (2)], but now the resonator frequency is provided by $\omega_{LC} = (LC_{eq})^{-1/2}$ with $C_{eq} = C_T C/(C_T + C)$. The new expression of the filling factor is $F = (f_{21}N_{qw}L_{eff}/d)C_{eq}/C$; this expression takes into account the redistribution of the vacuum electric field of the circuit between both capacitors C and $C_T$. The light-matter Hamiltonian Eq. (2) is supplemented with the transport contribution, Eq. (12). However, the flux that intervenes in the Peierls phase factor is no longer the flux $\Phi$ of the inductor of the uncoupled LC circuit alone. Rather, we must introduce a more general quantity $\mathcal{F}$, which takes into account the increased complexity of the circuit. This is shown further in Sec. II E. Indeed, since the capacitors C and $C_T$ are in series, a charge tunneling in the junction will feel the influence charge of both capacitors. The quantum fluctuations of the matter harmonic oscillator, $\hbar\tilde{\omega}_{21}(p^\dagger p + 1/2)$, should also influence the tunneling process. To take correctly into account the overall quantum fluctuations of the polaritonic system, we must first establish an equivalent circuit mapping of our light-matter coupling Hamiltonian, as described in Sec. II D. As a first step, we establish the mapping of the quantum-well system into an equivalent circuit, as explained in the next section.

Before proceeding with these ideas, let us briefly discuss the experimental feasibility of this system. Recently, there have been experimental demonstrations of LC metamaterial resonators with typical sizes of the capacitor parts on the order of 100 nm [57]. The corresponding capacitance can be evaluated in order of magnitude to be $C \approx \varepsilon_0\varepsilon_{100} \text{ nm} \approx 12 \text{ aF}$, where we used the dielectric function of the material (GaAs) filling the capacitance, $\varepsilon \approx 14$. This leads to typical charging energy $e^2/2C \approx 6 \text{ meV}$. For a structure resonating at 2 THz we have $\hbar\omega_{LC} \approx 8.3 \text{ meV}$, which provides $\rho_0 \approx 0.7$. Very high values of $\rho_0$ can thus be obtained in principle with such structures. Further estimations of $\rho_0$ from electromagnetic simulations are provided in Appendix B.

Since we consider here tunnel junctions with very small electrical areas, another important issue to address is the feasibility of electrical measurements from such structures. In particular, the tunnel junction should be sufficiently transparent as to allow measurable electrical current. Such junctions can be designed from semiconductor heterostructures [58] where both the height and the width of the tunneling barrier can be designed. They rely on the same growth technology that is used to fabricate the quantum wells interacting with the meta-atom. In Appendix B, we also show that such junctions can be reliably designed for the observation of the DCB effect in the THz frequency range. A very appealing perspective is that electrical measurement can thus be performed on a single metamaterial resonator, whereas optical spectroscopy typically requires dense arrays of identical resonators in order to

<figure><img src="image_6.png" /><figcaption>FIG. 3. Mapping of the quantum-well excitation into an equivalent circuit.</figcaption></figure>

overcome the diffraction limit [59,60]. Thus, an all-electrical probe of the USC regime offers more degrees of freedom to explore meta-atoms with various geometries.

### C. Equivalent circuit for the quantum-well excitation

The use of the $P(E)$ theory is justified by the linearity of the equivalent circuit describing the excitations of the environment of the tunnel junction. The description of mesoscopic systems in terms of circuit elements is extremely natural as these systems are mainly concerned with the electronic transport across nanostructures [18]. In most cases, resistances, capacitances, and inductances can be inferred in a straightforward way from the geometry of the system. For the case of coherent electronic transport, the notion of resistance is revisited by quantum mechanics [61], but it still involves the direct current transport of electrons from one contact to another.

An electric circuit description is far less obvious for the case of optical excitations in quantum-confined systems such as semiconductor quantum wells. In that case, electrons oscillate at a very high frequency (1 THz and above) between two quantum states, and one is rather concerned with the magnitude of the resulting dipole that couples to the electromagnetic field, as expressed by Eq. (11). The majority of USC systems with quantum wells published in the literature up to now thus do not rely on any circuit description.

In the present case, where we consider the DCB phenomenon in the context of polariton excitations, it is important to provide a mapping of the optical intersubband electronic excitation into an equivalent electrical circuit. A circuit approach is also needed for the description of light-matter coupling in such systems. Such mapping, which is absent in the current literature, allows us to bridge the languages employed by two distinct scientific areas: mesoscopic physics and semiconductor optoelectronics.

Since the bosonized Hamiltonian $\hbar\tilde{\omega}_{21}(p^\dagger p + 1/2)$ describes a quantum harmonic oscillator, one can always find a mapping into an equivalent circuit. This mapping is illustrated in Fig. 3; we will show that it is not a formal mapping, but the elements of the circuit actually have a clear physical meaning. The circuit consists of two capacitors $C_{12}$, $C_P$ and a kinetic inductance $L_{kin}$ which will be defined later. Anticipating the
discussion in Sec. III, we also introduce a resistive element $R_{kin}$. To establish the contributions $C_{12}$ and $C_p$, we are guided by the expression of the matter excitation, $\tilde{\omega}_{21}^2 = \omega_{12}^2 + \omega_p^2$, which suggests splitting the total capacitance of the circuit into a quantum contribution $C_{12}$ in series with a geometric or “plasmon” part $C_p$:

$$ \tilde{\omega}_{21}^2 = \omega_{12}^2 + \omega_p^2 = \frac{1}{L_{\text{kin}}C_{12}} + \frac{1}{L_{\text{kin}}C_p}. \quad (20) $$

The parameters $L_{kin}$, $C_p$, and $C_{12}$ can be chosen in a unique way by simple physical considerations. Let us first consider the case of a single quantum well, $N_{qw} = 1$, whereas the results will be immediately generalized for any number of identical quantum wells. Indeed, multiple quantum wells are used in order to obtain better geometrical overlap with the electric field of the LC resonator. The optical polarization corresponds to charge oscillations over a typical thickness $L_{eff}$, where positive/negative charges accumulate next to the edges of the well (see also Fig. 5(b) from Ref. [49]). We can thus identify the quantum well with a parallel plate capacitor of a thickness $L_{eff}$, such as $C_p = \varepsilon\varepsilon_0S/L_{eff}$. This identification is also compatible with the expression of the light-matter interaction energy as explained further. The kinetic inductance is then $L_{kin} = 1/C_p\omega_p^2 = m^*L_{eff}^2/N_e e^2$. We can introduce an effective three-dimensional (3D) electronic density $n_{3D} = N_e/(S \times L_{eff})$, such that the expression $L_{kin} = m^*L_{eff}/n_{3D}Se^2$ has a clear physical meaning: this is the kinetic inductance of a rectangular conductor $S \times L_{eff}$, where the current flows along the thickness $L_{eff}$. Indeed, this is the direction of the polarization microcurrent within the quantum wells [49].

Having established the expression for the kinetic inductance $L_{kin}$, we can now determine simply $C_{12}$ from the definition $C_{12} = 1/L_{kin}\omega_{12}^2$. To grasp the physical meaning of that quantity, let us consider as an example the transition between the first and second level of an infinite quantum well of a thickness $L_{qw}$. In that case, the transition energy is $\hbar\omega_{21} = 3\pi^2\hbar^2/2m^*L_{qw}^2$ and the effective thickness of the plasma oscillations is $L_{eff} = (3/5)L_{qw}$ [62]. We obtain the expression $C_{12} = S(m^*e^2/\pi\hbar^2)\beta(E_F/\hbar\omega_{12})$ with $\beta = 1.18$ a numerical factor. Here we expressed the electronic density through the Fermi energy $E_F$. We recognize $C_{12}$ to be in the form of a quantum capacitance [63] with a degeneracy factor $g = \beta(E_F/\hbar\omega_{12})$. The latter is on the order of unity, $g \approx 1$, for highly doped quantum wells where $E_F \le \hbar\omega_{21}$. A clearly similar expression can be derived for any confining potential, with only the numerical parameter $\beta$ changing value as a function of the shape of the electronic wave functions.

In the case of a stack with $N_{qw}$ identical quantum wells, the quantity $L_{eff}$ must be replaced by $N_{qw} \times L_{eff}$. On the other hand, neither of the frequencies $\omega_{12}$ and $\omega_p$ depends on $N_{qw}$. To keep Eq. (20) valid, the following expressions for the equivalent circuit lumped elements must be adopted:

$$ \begin{aligned} C_p &= \frac{\varepsilon\varepsilon_0 S}{N_{qw} L_{eff}}, & L_{\text{kin}} &= \frac{m^* N_{qw} L_{eff}}{n_{3D} S e^2}, \\ C_{12} &= \frac{Sm^* e^2}{\pi \hbar^2} \beta \frac{E_F}{\hbar \omega_{21}} \frac{1}{N_{qw}}. & & \end{aligned} \quad (21) $$

The dependence of $C_p$, $C_{12}$, and $L_{kin}$ on $N_{qw}$ is definitely in line with our physical interpretation, where the stack of quantum wells can be considered as capacitors and inductors in series. With the help of these quantities, we can now define quadrature operators for the matter equivalent circuit: $\hat{\phi} = i\sqrt{\hbar\tilde{\omega}_{21}}L_{kin}/2(p^\dagger - p)$ and $\hat{q} = \sqrt{\hbar\tilde{\omega}_{21}}\tilde{C}_{12}/2(p^\dagger + p)$. Here $\tilde{C}_{12} = C_{12}C_p(C_{12} + C_p)$, and the electronic Hamiltonian is re-expressed as

$$ \hbar\tilde{\omega}_{21} \left( p^\dagger p + \frac{1}{2} \right) = \frac{\hat{\phi}^2}{2L_{\text{kin}}} + \frac{\hat{q}^2}{2\tilde{C}_{12}}. \quad (22) $$

Using the relation $\tilde{\omega}_{21}^2\tilde{C}_{12} = \omega_p^2 C_p$, the charge variable can also be expressed as $\hat{q} = \sqrt{\hbar\omega_p^2 C_p/\tilde{\omega}_{21}}(p^\dagger + p)$, which will be used in the next section.

### D. Full circuit of the USC system

The mapping of the fundamental electromagnetic mode of the meta-atom structure into an equivalent circuit is straightforward. In the simplest approach, we can consider the structure from Fig. 2 as two capacitors $C$ and $C_T$ in series connected with an inductance $L$. This model can be refined if necessary (see Ref. [64] and Appendix B), but this configuration is already enough to understand the physics of the system. The presence of loss in the circuit is accounted for by a resistance $R$ in series with the inductor $L$; for the moment, we shall consider a lossless circuit with $R = 0$. The quantum Hamiltonian of the circuit alone is then, as discussed previously, $\hat{H}_{LC} = \hat{Q}^2/2C_{eq} + \hat{\Phi}^2/2L$, with $C_{eq} = C_T C/(C_T + C)$.

The full equivalent circuit of the polaritonic system is shown in Fig. 4(a). Following the scheme in Fig. 2, as well as experimental realizations [46], the plasmon capacitor $C_p$ is embedded in the capacitor $C$ of the meta-atom resonator. This picture is compatible with the fact that the light-matter interaction in these systems is mediated by the collective electronic polarization that is quantified by the capacitor $C_p$. On the other hand, the quantum contribution $C_{12}$ is taken outside the interacting region. To verify that this configuration indeed accurately reproduces the polariton system, we will momentarily adopt a classical description, defining $q$ and $Q$ as the charges induced on the capacitor plates $C$ and $C_p$, respectively, as illustrated in Fig. 4(b). The total electric field in the bulk of the capacitor $C$ is then a superposition of the fields created from both capacitors $C$ and $C_p$. This is shown in Fig. 4(c), where we treat the two capacitors as parallel plate capacitors, and the electric fields are considered to be homogeneous. The electric field in the capacitor $C$ is then $Q/S\varepsilon\varepsilon_0$ spanning over a distance $d$, and the electric field from $C_p$ is $q/S\varepsilon\varepsilon_0$ spanning over a distance $N_{qw}L_{eff}$. The voltage drop $V_{12}$ between the terminals 1 and 2 defined in Fig. 4(b) is then the area below the red curve in Fig. 4(c). Using the graph in Fig. 4(c), we obtain $V_{12} = Qd/S\varepsilon\varepsilon_0 + qN_{qw}L_{eff}/S\varepsilon\varepsilon_0$, which is also equal to $V_{12} = Q/C + q/C_p$. The electrostatic energy of the two nested capacitors is then evaluated at $W_C = \int_0^Q V_{12}dQ = Q^2/2C + qQ/C_p$. Considering the electrostatic energy stored in the other capacitors, $C_T$ and $C_{12}$, we obtain the following expression for the total electrostatic energy:

$$ W_{\text{elec}} = \frac{q^2}{2C_{12}} + \frac{q^2}{2C_p} + \frac{Q^2}{2C} + \frac{Q^2}{2C_T} + \frac{qQ}{C_p}. \quad (23) $$
<figure><img src="image_7.png" /><figcaption>FIG. 4. (a) Full equivalent circuit, which is now a combination between the equivalent circuit of the matter excitation Fig. 3 and that of the LC mode. To reproduce the light-matter interaction energy, we consider the capacitors C and Cₚ as nested. (b) Charges induced on the two nested capacitors. (c) Electric field distribution in the nested capacitors.</figcaption></figure>

In the above equation, we recognize the contribution from each individual circuit, whereas the last term on the right-hand side describes the electrostatic coupling between the two circuits. We now return to the quantum description and express this term with the help of the quantum variables $\hat{Q} = \sqrt{\hbar\omega_{LC}C_{eq}/2}(a^{\dagger} + a)$ and $\hat{q} = \sqrt{\hbar\omega_{p}^{2}C_{p}/\tilde{\omega}_{21}}(p^{\dagger} + p)$. As a result, we recognize the light-matter coupling Hamiltonian:

$$ \hat{H}_{\text{int}} = \frac{\hat{q}\hat{Q}}{C_p} = \frac{\hbar\omega_p}{2} \sqrt{\frac{\omega_{LC} C_{eq}}{\tilde{\omega}_{21} C_p}} (a^{\dagger} + a)(p^{\dagger} + p). \quad (24) $$

To complete the identification with Eq. (2), we define the filling factor to be $F = C_{eq}/C_p$, which can be rewritten as $F = (C/C_p)(C_{eq}/C)$. Using Eq. (21) and $C = \varepsilon\varepsilon_0S/d$, we have $C/C_p = N_{QW}L_{eff}/d$; we thus recover the expression of the filling factor from Sec. II B. The equivalent circuit described in Fig. 4 thus leads to the same quantum Hamiltonian as the one discussed in Sec. II A.

## E. Peierls substitution for the USC circuit

We now establish the $P(E)$ function for the polaritonic system based on the equivalent circuit. To this end, we replace

<figure><img src="image_8.png" /><figcaption>FIG. 5. (a) Differential conductance of the system normalized on the conductance quantum $G_K = e^2/h$ ($h = 2\pi\hbar$) as a function of the applied bias U. The blue curve corresponds to the uncoupled system, $\omega_p = 0$, while the red curve corresponds to the polaritonic system with a coupling parameter, $\omega_p = 1.2$ THz. (b) Plot of the amplitudes of the $P(E)$ function for the case of the uncoupled (blue) and the coupled system (red). The bars are equal to the prefactors of the Dirac $\delta$'s in Eq. (39).</figcaption></figure>

the flux $\hat{\Phi}$ in the tunneling coupling operator $\hat{T}$, Eq. (13), with a more general variable $\hat{\mathcal{F}}$, which accounts for the complexity of the circuit. Our approach is strongly inspired from Ref. [18], where the flux autocorrelation function on an inductance is related to the real part of the impedance seen by the inductance. In this vision, the junction probes the quantum fluctuations of the quantity $\int V_{ab}(t)dt$, where $V_{ab}$ is the voltage across the junction [55].

Let us consider a very general case where the electromagnetic environment of the tunnel junction inserted between the nodes a and b is an arbitrary linear circuit composed of capacitors and inductances. Any such circuit can be mapped into a system of coupled harmonic oscillators, with a finite number N of normal modes at frequencies $\omega_j$ that can be obtained by a Hopfield-Bogoliubov transformation. Alternatively, the circuit frequencies can be obtained as the poles of the impedance $Z_{ab}(\omega)$ seen by the junction [65]. The circuit Hamiltonian $\hat{H}$ can then be decomposed into the form $\sum_j \hbar\omega_j(a_j^\dagger a_j + 1/2)$, where all the circuit variables will appear as linear combinations of the operators $a_j$ and $a_j^\dagger$. Then it is always possible to define a quantity $\hat{\mathcal{F}}$ homogeneous to a flux, such that the voltage observable $\hat{V}_{ab}$ between the nodes a and b is expressed as a time derivative by virtue of the Heisenberg relation:

$$ \hat{V}_{ab} = \frac{d\hat{\mathcal{F}}}{dt} = \frac{1}{i\hbar}[\hat{\mathcal{F}}, \hat{H}]. \quad (25) $$
Using Eq. (25), we can define a generalized flux $\hat{\mathcal{F}}_j$ for each mode, which will be expressed as a linear combination of $a_j$ and $a_j^\dagger$ only. All the $\hat{\mathcal{F}}_j$ then commute with each other, and we can define a generalized multimode effective flux through the sum $\hat{\mathcal{F}} = \sum_j \hat{\mathcal{F}}_j$. The Peierls phase factor in the tunneling Hamiltonian is then defined as

$$ \hat{T} = \sum_{l,r} \tau_{l,r} c_l^\dagger c_r e^{i e \hat{\mathcal{F}} / \hbar}. \qquad (26) $$

Since $\hat{\mathcal{F}}_j$ commute with each other, we can treat each mode individually and split the Peierls phase factor into products $e^{i e \hat{\mathcal{F}}_j / \hbar}$. For each mode we can then follow exactly the same method as for the single-mode LC circuit [34], and obtain a corresponding DCB parameter as

$$ \rho_j = \frac{e^2}{\hbar^2} \langle \hat{\mathcal{F}}_j^2 \rangle. \qquad (27) $$

Here the mean value is taken in the ground state of the circuit, which is a tensor product of the ground states of the individual oscillators. The multimode zero-temperature $P(E)$ function will finally be expressed as

$$ P(E) = e^{-\prod_{j=1}^{N} \rho_j} \sum_{n_1, n_2, \dots, n_N}^{\infty} \frac{\rho_1^{n_1} \rho_2^{n_2} \cdots \rho_N^{n_N}}{n_1! n_2! \cdots n_N!} \delta \left( E - \sum_{j=1}^{N} n_j \hbar \omega_j \right). \qquad (28) $$

The function now takes into account the possibility to emit any number of circuit excitations during the tunneling process.

As an example, let us consider the LC resonator alone from Fig. 2, where we momentarily discard the presence of quantum wells in the capacitor C. The circuit seen by the junction with a capacitance $C_T$ is then the inductance L in series with the capacitor C. The voltage operator [Eq. (25)] becomes

$$ \hat{V}_{ab} = \frac{1}{i\hbar}[\hat{\Phi}, \hat{H}] + \frac{\hat{Q}}{C}. \qquad (29) $$

Here $\hat{\Phi}$ is the flux variable of the circuit [see Eq. (1)]. Using the commutator identity,

$$ a + a^{\dagger} = - \frac{1}{i\hbar\omega_{LC}} [i(a - a^{\dagger}), \hat{H}], \qquad (30) $$

as well as the expression for the charge variable from Eq. (1), we obtain the effective flux variable:

$$ \hat{\mathcal{F}} = \hat{\Phi}(1 - C_{\text{eq}}/C) = \hat{\Phi}/(1 + \eta), \qquad (31) $$

with $\eta = C_T/C$. The DCB parameter is then

$$ \rho_0 = \frac{e^2}{2C_T(1 + \eta)\hbar\omega_{LC}}. \qquad (32) $$

We can easily verify that in the case in which the capacitor C is replaced with a wire, $C \to \infty$ and $\eta \to 0$, we recover the expression Eq. (19) for the circuit with the tunnel junction capacitance alone. For any finite value for the capacitance C, the DCB parameter is always smaller than Eq. (19), which is quite intuitive, as one would expect that additional capacitive elements will tend to reduce the charge fluctuations seen by the junction. Such a reduction has been discussed in Ref. [66] in the context of multijunction circuits. Equation (32) indicates

that the case with a large capacitance $C \gg C_T$ is favorable for the DCB regime. Eventually the only request is that the tunnel junction capacitance $C_T$ is made very small. However, considering the case in which C is filled with quantum wells, the geometrical overlap factor between the electromagnetic mode and the electronic polarization is $F = C_{eq}/C_P$ and it is proportional to $C_T/(C_T + C)$. To achieve USC, we typically need values of F as close to unity as possible, therefore the capacitor C cannot be made very large with respect to $C_T$. In the following, for the sake of example, we shall consider the situation in which $C = C_T$, where both DCB and USC regimes can be realized within the same meta-atom.

Considering the full polaritonic circuit from Fig. 4, the voltage observable is

$$ \hat{V}_{ab} = \frac{1}{i\hbar}[\hat{\Phi}, \hat{H}] + \frac{\hat{Q}}{C} + \frac{\hat{q}}{C_P}. \qquad (33) $$

It is convenient to work directly with the polariton operators, Eqs. (3) and (5). Noting $j = +, -, $ the charge and flux variables are expressed as

$$ \hat{\Phi}_j = \sqrt{\frac{L\hbar\omega_{LC}}{2}} i(\Pi_j - \Pi_j^\dagger)(x_j + y_j), \qquad (34) $$

$$ \hat{Q}_j = \sqrt{\frac{C_{\text{eq}}\hbar\omega_{\text{LC}}}{2}} (\Pi_j + \Pi_j^\dagger)(x_j - y_j), \qquad (35) $$

$$ \hat{q}_j = \sqrt{\frac{C_P \hbar \omega_P^2}{2 \tilde{\omega}_{21}}} (\Pi_j + \Pi_j^\dagger)(m_j - t_j). \qquad (36) $$

Using bosonic identities such as the one from Eq. (30) and the algebra of the Hopfield coefficients, one obtains the following result for the effective flux operator $\hat{\mathcal{F}}_j$ for each polariton mode:

$$ \hat{\mathcal{F}}_j = \sqrt{\frac{L\hbar\omega_{LC}}{2}} i(\Pi_j - \Pi_j^\dagger) \frac{(x_j - y_j) \omega_{LC}}{1 + \eta \omega_j}. \qquad (37) $$

Details for this calculation are provided in Appendix A. Eventually, the DCB parameters for the two polariton states are obtained as

$$ \rho_{\pm} = \rho_0 \gamma_{\pm} \left( \frac{\omega_{LC}}{\omega_{\pm}} \right)^3. \qquad (38) $$

Here $\rho_0$ is the DCB parameter of the uncoupled circuit as provided by Eq. (32). The $P(E)$ function of the polaritonic system is then

$$ P(E) = e^{-\rho_{+} - \rho_{-}} \sum_{n,m}^{\infty} \frac{\rho_{+}^{n} \rho_{-}^{m}}{n! m!} \delta(E - n\hbar\omega_{+} - m\hbar\omega_{-}). \qquad (39) $$

The normalized conductance $G(U)/G_K$ of the junction is obtained from Eq. (39) by replacing the Dirac $\delta$'s with Heaviside functions where the energy E is replaced with the applied bias, and the excitation energies are divided by the elementary charge.

In Fig. 5(a) we illustrate the conductance of the system [Eq. (17)] for the case of an uncoupled system ($\omega_P = 0$, blue curve) and a polaritonic system (red curve). The $P(E)$ function is illustrated in Fig. 5(b), where the height of the bars is proportional to the prefactors of the Dirac $\delta$'s in
Eq. (39). The parameters for these plots are $\omega_{21} = 2$ THz, $\omega_P = 1.2$ THz, $\omega_{LC} = \tilde{\omega}_{21}$, $F = 0.7$, and $\rho_0 = 0.4$ (we have $C_T = 10$ aF and $\eta = 1$). For the uncoupled system, the low voltage conductance is suppressed and appears as a series of steps at the frequencies that are multiples of $\omega_{LC}$ and with amplitudes $e^{-\rho_0}\rho_0^n/n!$. Similar behavior is observed for the coupled system, where the steps appear at multiples of the polariton frequencies $n\omega_+ + m\omega_-$. The sizes of the steps are expressed from the parameters $\gamma_\pm$, which also quantify the quantum fluctuations of the charge and flux through Eqs. (8) and (9).

Experimentally, one would measure the device current voltage I-U characteristics $I = I(U)$ and obtain the plot from Fig. 5(a) by computing the differential conductance, $dI/dU$. The $P(E)$ function can be evaluated from the second derivative $d^2I/dU^2$. Both plots will feature discontinuities for bias values $U = \hbar(n\omega_+ + m\omega_-)/e$ from which the frequencies of the coupled states can be deduced. Therefore, the device proposed here, on a very basic level, can be used to perform spectroscopy of the light-matter coupled states in all electrical measurements. As mentioned previously, such measurements can be performed on a single device, which is more onerous in optical experiments [67], rather than using large arrays that are typically required in order to obtain a good optical signal in THz systems [59,60].

Furthermore, measurements of the peak amplitudes observed in the second derivative $d^2I/dU^2$ can be used to evaluate the coefficients $\gamma_\pm$, further providing the amount of the squeezing and antisqueezing effects through Eqs. (8) and (9). The DCB thus offers an experimental platform to probe quantum correlations specific to the USC regime. This idea is further developed in Sec. IV. The features of the conductance displayed in Fig. 5 are property only on the ground state of the coupled system. Therefore, as long as we neglect the backaction of the tunnel junction on the system, the DCB probe can be seen as a nondemolition measurement of the quantum-optical properties of the ground state of the USC system.

The experimental realization of this program is clearly demanding, and among other things the most fundamental correction to the purely quantum picture developed so far is the presence of dissipation. Namely, various loss mechanisms will lead to an enlargement of the peaks that we expect in the $P(E)$ function from Fig. 5(b). It is important to understand how the peak amplitudes, which we use for the estimation of the quantum correlations, are affected by linewidth broadening. In principle, dissipation can still be treated within a purely quantum approach by using the Caldeira-Leggett theory [18,68] to describe the resistive elements of the circuits. However, the circuit $P(E)$ theory, which links the $P(E)$ function to the real part of the circuit impedance [34], provides a more practical approach to the problem, as explained in the next section.

### III. MODEL INCLUDING DISSIPATION

#### A. The $P(E)$ function with losses

We will consider three different dissipation paths in the system, as illustrated in Fig. 6. The LC resonator is subject to a

<figure><img src="image_9.png" /><figcaption>FIG. 6. Schematics indicating the various dissipation mechanisms for the system together with the definitions of the corresponding coupling rates.</figcaption></figure>

radiation loss, with a rate $\Gamma_r$ that quantifies its coupling to the electromagnetic modes of free space. This type of loss can be used for probing the optical excitations of the system [69,70]. Actually, for the meta-atom architectures that we consider, the loss rate $\Gamma_r$ can be engineered through the geometry of the capacitive parts [57]. In particular, one can favor high-symmetry multipole designs that can cancel the radiation loss to a great extent [57]. The second type of loss of the LC circuit is the nonradiation loss due to the Ohmic dissipation that occurs in the inductance, quantified by a loss rate $\Gamma_{nr}$. Finally, we also consider the loss of the material excitation, with a rate $\gamma_{12}$, which describes the relaxation of the electronic polarization of the electron gas [69].

To account for these dissipation rates in our model, we supplement the equivalent circuit from the previous sections with resistive elements. The effective resistance in the matter equivalent circuit, $R_{kin}$ (Fig. 3), is expressed from the kinetic inductance as $R_{kin} = \gamma_{12}L_{kin}$. This relation is identical to the one that links the kinetic inductance to the resistance of a metal conductor described by the Drude model [71]. However, we use here the polarization relaxation time, $1/\gamma_{12}$, instead of the collision time, as would be the case for the Drude model. Similarly, the resistance of the LC resonator is expressed from the loss rates as $R = L(\Gamma_r + \Gamma_{nr})$. This expression is consistent with the quantum Langevin equation approach for electrical circuits [72].

The inclusion of a radiation loss channel allows us to probe the system by optical experiments, which is very useful in spectroscopy. Assuming a finite radiation loss for an identical device probed in electrical measurements allows us to compare the spectroscopic experiments with the DCB probe. Therefore, in Sec. III B, we also provide the reflectivity and absorption coefficients of the system. The latter have been obtained without relying on the rotating-wave approximation typically employed to take into account interactions with dissipative baths. Our method, described in Sec. III B, is thus consistent with the USC regime, where now the rotating-wave
approximation is abandoned not only for the treatment of the light-matter interaction, but also for its coupling with the environment.

In the P(E) theory approach it is sufficient to consider the classical dynamics of the circuit. From Fig. 4 we can easily infer the voltage drops $V_{12} = q/C_p + Q/C$ and $V_{34} = (Q + q)/C_p$ between the terminals defined in Fig. 4(b). The circuit diagram from Fig. 4(a) provides the following equations of motion:

$$V_{12} = -L \frac{d^2 Q}{dt^2} - R_{\text{tot}} \frac{dQ}{dt} - \frac{Q}{C_T} = \frac{q}{C_p} + \frac{Q}{C}, \quad (40)$$

$$V_{34} = -L_{\text{kin}} \frac{d^2 q}{dt^2} - R_{\text{kin}} \frac{dq}{dt} - \frac{q}{C_{12}} = \frac{q}{C_p} + \frac{Q}{C_p}. \quad (41)$$

As we are interested in the high-frequency dynamics, we have neglected the constant voltage source in Fig. 4(a) used to probe the tunnel junction conductance. In the harmonic regime $q, Q \propto \text{Re}(e^{-i\omega t})$ the following characteristic equation of the system (40) and (41) is obtained:

$$(\omega^2 - \omega_{\text{LC}}^2 - i\Gamma_{\text{tot}}\omega)(\omega^2 - \tilde{\omega}_{21}^2 - i\gamma_{12}\omega) = F\omega_p^2\omega_{\text{LC}}^2. \quad (42)$$

In the limit of vanishing loss in the system, this equation is identical to the polariton dispersion relation (6). Its complex solutions $\omega_{\pm} + i\Gamma_{\pm}$ provide both the polariton frequencies $\omega_{\pm}$ as well as their linewidths $\Gamma_{\pm}$ as a function of the corresponding quantities of the uncoupled systems.

From Fig. 4(a) and Eqs. (40) and (41) it is straightforward to determine the impedance seen by the tunnel junction $Z_{ab} = V_{ab}/i\omega Q$:

$$Z_{ab}(\omega) = \frac{1}{i\omega C_T} \left[ 1 - \frac{\omega_{\text{LC}}^2 \Delta_{21}(\omega)}{(1 + \eta)\Pi(\omega)} \right], \quad (43)$$

$$\Pi(\omega) = \Delta_{12}(\omega)\Delta_{LC}(\omega) - F\omega_p^2\omega_{\text{LC}}^2, \quad (44)$$

$$\Delta_{12}(\omega) = \tilde{\omega}_{21}^2 - \omega^2 + i\gamma_{12}\omega, \quad (45)$$

$$\Delta_{\text{LC}}(\omega) = \omega_{\text{LC}}^2 - \omega^2 + i\Gamma_{\text{tot}}\omega. \quad (46)$$

The functions $\Delta_{12}(\omega)$ and $\Delta_{LC}(\omega)$ defined above are the characteristic functions of the plasmon and electromagnetic oscillators, respectively, as their zeros provide the eigenfrequencies and the dissipation rates of the uncoupled oscillators. The function $\Pi(\omega)$ defined in Eq. (44) is the characteristic function of the coupled system, and Eq. (42) can be rewritten as $\Pi(\omega) = 0$. The complex poles of the impedance $Z_{ab}(\omega) \propto 1/\Pi(\omega)$ are thus exactly at the complex frequencies of the polariton states, as expected.

Following Ref. [34], the function P(E) at zero temperature is linked to the impedance $Z_{ab}(\omega)$ through the following integral equation:

$$EP(E) = 2 \int_{0}^{E} \frac{\text{Re}[Z_{ab}((E - E')/\hbar)]}{R_{K}} P(E') dE'. \quad (47)$$

Of particular interest is then the real part of the impedance Eq. (43):

$$\text{Re}Z_{ab}(\omega) = -\text{Im}\left[\frac{\Delta_{21}(\omega)}{\Pi(\omega)}\right] \frac{\omega_{\text{LC}}^2}{(1 + \eta)\omega C_T}. \quad (48)$$

We can use this expression to provide a numerical solution to Eq. (47). Further results will be provided in Sec. IV.

In the limit of vanishing losses, this approach should provide the same results as the quantum calculation from the previous sections. In that limit, we have $\Pi(\omega) = (\omega^2 - \omega_+^2)(\omega^2 - \omega_-^2)$ and the following identity holds:

$$\frac{1}{\Pi(\omega)} = \frac{1}{\omega_{\pm}^2 - \omega_{\mp}^2} \left[ \frac{1}{\omega^2 - \omega_{\pm}^2} - \frac{1}{\omega^2 - \omega_{\mp}^2} \right]. \quad (49)$$

Using the formula $\text{Im}(1/(\omega^2 - \omega_{\pm}^2 + 2i\omega_{\pm}\varepsilon)) = (\pi/2\omega_{\pm})\delta(\omega - \omega_{\pm})$ for $\varepsilon \to 0$ and the definitions in Eq. (7), the real part of the impedance is expressed as a sum of Dirac's:

$$\frac{2 \text{Re}(Z_{ab}(\omega))}{R_K} = \frac{\pi}{C_T(1 + \eta)R_K} \frac{\omega_{\text{LC}}^2}{\omega_{\pm}^2} \gamma_+ \delta(\omega - \omega_+) \\ + \frac{\pi}{C_T(1 + \eta)R_K} \frac{\omega_{\text{LC}}^2}{\omega_{\mp}^2} \gamma_- \delta(\omega - \omega_-). \quad (50)$$

This $\delta$-function decomposition makes the resolution of Eq. (47) very straightforward. For that, we take a trial function $P(E)$ as the sum $\sum_{n,m} A_{n,m} \delta(E - n\hbar\omega_+ - m\hbar\omega_-)$. The convolution in Eq. (47) then allows us to determine the unknown coefficients $A_{n,m}$ iteratively, and eventually we arrive exactly at the P(E) function described by the expressions (39) and (32). In the case of vanishing loss, the P(E) function obtained from the integral equation (47) is thus compatible with the quantum approach based on the generalized Peierls substitution in the tunneling Hamiltonian from Sec. II. Note that a general theory that links the poles of the impedance of the environment and the peak amplitudes of the P(E) function has been developed in Ref. [65].

Let us consider at present the case for weak but non-vanishing loss, where an explicit formula for P(E) can be provided as well. To this end, we replace the $\delta$ functions in Eq. (50) by normalized Lorentzians, $L_{\Gamma_{\pm}}(\omega - \omega_{\pm}) = \Gamma_{\pm}/[\pi((\omega - \omega_{\pm})^2 + \Gamma_{\pm}^2)]$. We can use a similar trial function for P(E) to that used above, but now the Dirac's are replaced by Lorentzians with linewidths that are unknown a priori. Equation (47) can then be solved iteratively by exploiting the convolution rule $L_{\Gamma_1}(\omega) * L_{\Gamma_2}(\omega) = L_{\Gamma_1+\Gamma_2}(\omega)$. Assuming sufficiently small linewidths, such that only nearly resonant Lorentzians are considered to contribute to the integral of (47), we obtain the following expression of P(E):

$$P(E) \approx e^{-\rho_{+} - \rho_{-}} \sum_{n,m}^{\infty} \frac{\rho_{+}^{n} \rho_{-}^{m}}{n! m!} L_{n\Gamma_{+} + m\Gamma_{-}}(E - n\hbar\omega_{+} - m\hbar\omega_{-}). \quad (51)$$

This expression predicts that the successive peaks that we expect from the measurement of the P(E) function have progressively increasing broadening. This is an important conclusion as it shows that particular care must be taken when extracting the coefficients $\gamma_{\pm}$ from the peak amplitudes. These results will be further examined in Sec. IV.

## B. Reflectivity and absorption coefficients

In this section, we consider the system as two coupled harmonic oscillators driven by an external field. Our aim is to
obtain the reflectivity and absorption coefficients that would correspond to a spectroscopic probe of the coupled system. However, since we are considering the USC regime, we would like to avoid any secular approximations and solve for the exact input-output dynamics of the system. To accomplish this task, we will work with the dynamic evolution of the oscillator quadratures rather than with that of the creation and annihilation operators. The latter have the merit to have a clear physical interpretation. We can furthermore identify quadratures that dissipate energy and those that do not. As an example, consider a classical mechanical pendulum evolving in a medium with fluid friction. This system is described by two quadratures: the position x and momentum p. Clearly, loss of energy occurs through the viscous friction, which is proportional to the velocity (momentum). Another example is an LC circuit, where energy is dissipated through the current (flux) and not the charge.

To establish the dynamics, the Hamiltonian is supplemented with a driving term, as described in Ref. [73]. We then use the Heisenberg equation of motion $i\hbar d\hat{O}/dt = [\hat{O}, \hat{H}]$, where $\hat{O}$ is any quadrature of interest. Since we consider permanent harmonic excitation, the two oscillators evolve into a tensor product of coherent states $|\alpha\rangle|\pi\rangle$, with $a|\alpha\rangle = \alpha|\alpha\rangle$ and $p|\pi\rangle = \pi|\pi\rangle$. The system evolution is then expressed as coupled first-order differential equations for the corresponding mean values of the observables $\langle\hat{O}\rangle$ in that state. After adding the loss terms, we arrive at the following equation set:

$$ \frac{dA}{dt} = \omega_{LC}E, \qquad (52) $$

$$ \frac{dE}{dt} = -\omega_{LC}A - 2\Omega_R S - \Gamma_{tot}E + 2\sqrt{\Gamma_{rad}}I_{in}, \qquad (53) $$

$$ \frac{dS}{dt} = \tilde{\omega}_{21}D, \qquad (54) $$

$$ \frac{dD}{dt} = -\tilde{\omega}_{21}S - 2\Omega_R A - \gamma_{21}D, \qquad (55) $$

$$ I_{out} = -I_{in} + \sqrt{\Gamma_{rad}}E. \qquad (56) $$

Here we introduced the notations $\Omega_R = (\omega_p/2)\sqrt{F\omega_{LC}/\tilde{\omega}_{21}}$, $A = \langle(a + a^\dagger)\rangle$, $E = i\langle(a^\dagger - a)\rangle$, $S = \langle(p + p^\dagger)\rangle$, and $D = i\langle(p^\dagger - p)\rangle$, and $I_{out}$ and $I_{in}$ are, respectively, the output and input field applied on the coupled system. The quadratures A and E are directly proportional to the charge and flux of the LC circuit, while the quadratures S and D are related to the charge and flux of the quantum-well “circuit” (Fig. 3). This justifies the choice of the dissipation terms in Eqs. (53) and (55): in both cases, energy is dissipated through the currents of the equivalent circuits. The input/output relations in Eq. (56) can also be obtained by modeling the free-space modes as a transmission line [54]. These equations thus generalize the input-output relations for polaritonic systems [74], as well as the equations of the coupled mode theory [75]. Similar results can also be obtained in the framework of the quantum Langevin equations [72].

We can further show that the above equations obey the following conservation law in a permanent harmonic regime where $d(E^2 + A^2)/dt = 0$:

$$ I_{in}^2 = I_{out}^2 + \Gamma_{nr}E^2 + 2\Omega_RSE. \qquad (57) $$

Energy is thus dissipated through the nonradiative loss of the LC resonator $\Gamma_{nr}E^2$ and absorption in the matter part $2\Omega_RSE$. We can thus define reflectivity and absorption coefficients:

$$ R = 1 - (\Gamma_{nr}E^2 + 2\Omega_RSE)/I_{in}^2, \qquad (58) $$

$$ \eta = 2\Omega_RSE/I_{in}^2. \qquad (59) $$

Since the above equations are linear, we can work with complex amplitudes $O(t) = \text{Re}[\tilde{O}(\omega)e^{i\omega t}]$. Then from Eqs. (58) and (59) we obtain the frequency-dependent absorption and reflectivity coefficients:

$$ \eta(\omega) = F \frac{4\Gamma_{rad}\gamma_{21}\omega^2\omega_{LC}^2\omega_p^2}{|\Pi(\omega)|^2}, \qquad (60) $$

$$ R(\omega) = 1 - \frac{4\omega^2\Gamma_{rad}}{|\Pi(\omega)|^2} [\Gamma_{nr}|\Delta_{12}(\omega)|^2 + \gamma_{21}F\omega_p^2\omega_{LC}^2]. \qquad (61) $$

The absorption efficiency from Eq. (60) is well proportional to the intersubband loss $\gamma_{21}$ and the radiation loss $\Gamma_{rad}$: it vanishes if either of the two loss mechanism is absent. Further, it is proportional to the square of the plasma frequency, which is also proportional to the total number of dipoles in the LC resonator. From Eq. (61) the reflectivity goes to unity if either of the loss mechanism vanishes. Both the reflectivity contrast $1 - R(\omega)$ and the absorption efficiency are inversely proportional to the polaritonic denominator $|\Pi(\omega)|^2$: they are maximized at the polariton frequencies $\omega_\pm$. Using the approach of the previous paragraph, we can relate the peak weights for the reflectivity or absorption to the Hopfield coefficients of the coupled system $\gamma_\pm$. Note, however, that, especially in the case of the reflectivity coefficient, the peak weights also depend on the various loss mechanisms present in the system. Finally, the expressions Eqs. (60) and (61) are free from any secular approximations, and therefore they are valid for arbitrary values of the loss parameters as well as the plasma frequency. In particular, the absorption efficiency and the contrast both vanish at very low frequencies as $\omega^2$.

## IV. RESULTS

We now apply our results from the previous section in order to compare systematically the optical response of a USC system with transport experiments in the DCB regime. For the latter we compute the $P(E)$ function by numerically integrating Eq. (47) with system parameters that will be indicated further. Bearing in mind that $P(E)$ should be obtained from conductivity measurements, we will actually consider the function $P(eU)$, which is plotted as a function of the applied bias $U$ on the junction. Nevertheless, without ambiguity we will keep using the notation “$P(E)$” for the vertical axes of the plots.

In Figs. 7 and 8 we compare the spectroscopic and electrical measurements of a polaritonic system. In Fig. 7 we consider the case in which the LC resonance is matched with the energy of the bare electronic transitions of the wells, $\hbar\omega_{LC} = \hbar\omega_{21} = 8.3$ meV, while the plasma energy $\hbar\omega_p$ is increased from 0 to 6 meV. We use the following values of the damping parameters: $\hbar\gamma_{12} = 0.6$ meV, $\hbar\Gamma_{nr} = 0.4$ meV, and $\hbar\Gamma_{rad} = 0.8$ meV. As seen in Fig. 7, for $\hbar\omega_p = 0$ there is no coupling between the LC resonator and the quantum wells;
<figure><img src="image_10.png" /><figcaption>FIG. 7. Comparison between optical and electrical response of a polaritonic system. (a), (b), (c) Reflectivity contrast $1 - R(\omega)$ [Eq. (61) for progressively increasing plasma energy $\hbar\omega_p$], from the weak (c) to the USC regime (a). The LC resonator energy is matched with the bare transition energy in the wells, $\hbar\omega_{LC} = \hbar\omega_{21} = 8.3$ meV. Parts (a) and (c) are spectra for $\hbar\omega_p = 6$ and 0 meV, respectively. (e), (d), (f) The tunneling probability density function $P(E)$ of the polaritonic system, from the weak (f) to the USC regime (d), as obtained numerically from Eq. (47).</figcaption></figure>

in that case, we observe only a single peak that corresponds to the absorption of the LC resonance [Fig. 7(c)]. When the coupling constant $\hbar\omega_p$ is progressively increased [Fig. 7(b)], the LC resonance is split into two polariton modes, $\hbar\omega_-$ and

<figure><img src="image_11.png" /><figcaption>FIG. 8. Polariton dispersion probed in optical and electrical experiments. (a) Reflectivity spectra for variable frequency $\omega_{LC}$ of the uncoupled LC resonator. (b) Plot of the corresponding $P(E)$ function. (c) Polariton dispersion probed in reflectivity for a polaritonic system (see the main text for detailed parameters). (d) $P(E)$ function for the polaritonic system.</figcaption></figure>

$\hbar\omega_+$. The maximum value $\hbar\omega_p = 6$ meV has been chosen so that the Fermi level of the electron gas in the quantum well lies approximately just below the subband 2, $E_F \approx \hbar\omega_{12}$.

In the panels (d), (e), and (f), the same system is probed in transport measurements, and the quantity $\log(P(E))$ is plotted as function of the applied bias $U$ on the junction. In this case we consider that the radiation loss of the LC resonator is zero, $\Gamma_r = 0$, in order to avoid unwanted vacuum fluctuation from free space. This is always possible in our structures by an appropriate engineering of the capacitor parts [57]. For the LC resonator alone [$\hbar\omega_p = 0$, Fig. 7(f)], we observe the main peak at $\hbar\omega_{LC}$ as well as the replicas $n\hbar\omega_{LC}$, which are a signature of the DCB regime. The capacitance of the circuit is $C = C_T = 10$ aF, and $\rho_0 \approx 0.4$. As the coupling parameter $\hbar\omega_p$ is progressively increased [Fig. 7(e)], we start observing peaks that correspond to the two polariton states and their replicas, $n\hbar\omega_- + m\hbar\omega_+$, which appear fully in the USC regime [$\hbar\omega_p = 6$ meV, Fig. 7(d)]. The upper bound for the energy range for this study has been limited to 25 meV, as we want to avoid the optical phonon resonance at 35 meV present in the semiconductor material that is typically used to fabricate the quantum wells [59].

In Fig. 8 we display the polariton dispersion that is probed in optical and electrical measurements. Now the LC resonance is swept between 2 and 20 meV [Fig. 8(a)], while the bare intersubband energy is kept constant, $\hbar\omega_{21} = 8$ meV. Note that the $\rho_0$ parameter from Eq. (32) is no longer a constant in that case; we assume that the capacitive parts of the circuits do not change, while the inductance is continuously varied [57]. In the case of the uncoupled resonator ($\hbar\omega_p = 0$ meV), we observe the typical DCB behavior with the replicas of the main peak [Fig. 8(b)]. The peak amplitudes of the $P(E)$ function are stronger for low-frequency resonators, as the $\rho_0$ parameter is higher in that case. At low bias the presence of resistors in the equivalent circuit also affects the $P(E)$ function [34], as commented further.

In the reflectivity maps of the polaritonic system with $\hbar\omega_p = 6$ meV, we observe the typical anticrossing behavior with a gap signaling the USC regime [Fig. 8(c)]. The anticrossing features become more complex in the map of the $P(E)$ function [Fig. 8(d)]. The polariton dispersion is replicated at higher energies by the DCB effect. Quite interestingly, we observe a “twisting” pattern between the upper polariton and the first replica of the lower polariton in Fig. 8(d). These results show that, in the THz range, the USC regime can be probed by transport experiments; it is very appealing to perform such studies on completely dark circuit modes ($\Gamma_r = 0$) which cannot be accessed by optical means, but that are particularly interesting in the DCB regime as the corresponding circuit quality factors are higher.

Already in Figs. 7(b) and 7(e) we observe an asymmetry of the polariton branches as the plasma frequency $\hbar\omega_p$ increases; this effect is due to the renormalization of the matter excitation contained in the formula $\tilde{\omega}_{21} = \sqrt{\omega_{21}^2 + \omega_p^2}$. It is interesting to examine the case of the parabolic potential, which can formally be described as the case $\tilde{\omega}_{21} = \text{const}$. This case corresponds to Kohn's theorem [76], which has been numerically proven for intersubband transitions [77] and observed experimentally [78,79]. In Fig. 9 we provide such a comparison: as
<figure><img src="image_12.png" /><figcaption>FIG. 9. DCB effect in a polariton system obtained with square quantum wells (a), (b) and parabolic quantum wells (c), (d). (a), (c) Reflectivity spectra for progressively increasing plasma energy $\hbar\omega_p$. (b), (d) The $P(E)$ function for both systems. (e), (f) Extrapolation to a coupled system that can feature quantum phase transition, QPT.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 10. Comparison of different models for the computation of $P(E)$ for the cases of (a) low-quality and (b) high-quality factors (see the text for the values). Full blue line: full model based on numerical integration of Eq. (47). Dashed red line: analytical expression from Eq. (51).</figcaption></figure>

seen from Figs. 9(c) and 9(d), the polariton asymmetry is reversed in that case with respect to the case of a square potential [Figs. 9(a) and 9(b)]. In particular, the lower polariton tends faster to lower frequencies with respect to the case of a square potential. This is interesting from the point of view of DCB, as the effective parameter $\rho = e^2/2C\hbar\omega_-$ increases faster for the lower polariton. Of particular interest is the asymptotic limit $\omega_- \to 0$, which must be accompanied by singularities of the vacuum fluctuations of the system [80]. To examine this point, we have extrapolated our model to the case of a parabolic potential $\tilde{\omega}_{21} = \text{const} = \Omega_0$ with a high filling factor, $F = 0.95$, and a plasma frequency reaching the maximum value $\omega_p \to \Omega_0$. The corresponding plots for the reflectivity contrast and the $P(E)$ function are shown in Figs. 9(e) and 9(f), where we also have considered smaller values for the dissipation parameters: $\gamma_{12} = 0.06$ meV, $\hbar\Gamma_{nr} = 0.04$ meV. As seen in Fig. 9(f), we observe strong suppression of the tunneling probability density as the critical point $\omega_- \to 0$ is approached. The critical point cannot be traversed in that case because of a no-go theorem for intersubband systems [49]; however, this limitation is not necessarily true for other systems [6,81]. From Fig. 9(f) we can argue that conductance measurements in the DCB regime could indicate the onset of a critical point, and more generally DCB could be used as an experimental probe for quantum phase transitions, a question that has not yet been addressed theoretically to our knowledge.

We now compare systems with different quality factors. The comparison is displayed in Fig. 10, where the system parameters are $C = C_T = 10$ aF, $\hbar\omega_p = 6$ meV, $\hbar\omega_{LC} = \hbar\tilde{\omega}_{21} = 10.2$ meV, and $\rho_0 = 0.4$. In Fig. 10(a) we consider the case of relatively strong dissipation, with $\hbar\Gamma_{nr} = 0.41$ meV (quality factor $Q \approx 20$ for the LC circuit, assuming only nonradiative loss) and $\hbar\gamma_{12} = 0.62$ meV. In Fig. 10(b), we consider weak dissipation, with loss parameters $\Gamma_{nr}$, $\gamma_{12}$ that are set to be ten times lower. In Fig. 10 the full blue lines correspond to the full circuit model based on the numerical integration of Eq. (47). The red dashed lines correspond to the $P(E)$ function obtained from the analytical expression Eq. (51). One of the effects of the dissipation, that is visible in Fig. 10, is the strong increase of the $P(E)$ function at low frequencies: this is due to the resistive element of the equivalent circuit [34,82]. This effect is not present in the simplified model from Eq. (51), and, as seen from Fig. 10(a), it leads to a substantial decrease of the spectral weight of the lower polariton peak at $\hbar\omega_-/e$. At higher frequencies, the analytical expression given by Eq. (51) (red dashed curve) closely matches the outcome of the full model. In particular, the gradual increase in the linewidths according to the formula $n\Gamma_+ + m\Gamma_-$ is clearly observed. Not surprisingly, for low dissipation the analytical expression from Eq. (51) is sufficient to model the expected features of the $P(E)$ function everywhere except in the low-voltage/energy limit.
<figure><img src="image_14.png" /><figcaption>FIG. 11. (a) Full line: exact parameters $\rho_{\pm}$ from Eq. (38) (case without dissipation) for gradually increasing plasma frequency $\omega_{p}$ (in THz). Dots: $\rho_{\pm}$ parameters extracted from fits of the $P(E)$ function for different quality factors of the system. (b) Charge quantum fluctuations enhancement factor [Eq. (8)], as well as the sum $\gamma_{+} + \gamma_{-}$ obtained from the values $\rho_{\pm}$ from (a) (open symbols). (c) Estimation of the charge antisqueezing from the second line of Eq. (63).</figcaption></figure>

We now analyze the information contained in the absolute values of the peaks of the $P(E)$ function. We consider the case in which $\omega_{LC} = \tilde{\omega}_{21} = 2$ THz and the only parameter that we vary is the plasma frequency $\omega_{p}$ while keeping $\tilde{\omega}_{21}$ constant [the case of the parabolic quantum well, Figs. 8(c) and 8(d)]. Using our model, we generate plots of the $P(E)$ function, such as those shown in Fig. 10, both in the regimes of weak and strong dissipation. The resulting curves are fitted with a sum of Lorentzian functions in order to extract both the amplitudes $A_{i}$ and the linewidths $\Gamma_{i}$ of the peaks of the $P(E)$ function. In this approach, there is a natural error margin that arises from fitting the peaks with small amplitudes or from overlapping peaks.

The parameter $\rho_{0}$ can be inferred from the case of an uncoupled resonator [see Fig. 7(f)] by comparing peak amplitudes which vary as $\rho_{0}^{n}/n!$: for instance, $\rho_{0} = 2A_{2}/A_{1}$, where $A_{1}$ and $A_{2}$ are, respectively, the amplitude of the main peak and the $2\hbar\omega_{LC}/e$ peak in Fig. 7(f). Once the parameter $\rho_{0}$ is known or measured, we proceed in the same way to evaluate the $\rho_{+}$ and $\rho_{-}$ parameters. In addition to the two polariton peaks at $\hbar\omega_{\pm}/e$, with amplitudes $A_{-}$ and $A_{+}$, it is sufficient to observe at least one replica: for instance, $2\hbar\omega_{-}/e$ with an amplitude $A_{2-}$. Then we can deduce $\rho_{-} = 2A_{2-}/A_{-}$ and subsequently $\rho_{+} = \rho_{-}A_{+}/A_{-}$. If more replicas are observed, the uncertainty in determining $\rho_{+}$ and $\rho_{-}$ is reduced. For instance, the peak at the voltage $\hbar(\omega_{+} + \omega_{-})/e$ can be useful: it fits better the simplified model for high dissipation [Fig. 10(a)] and it is less affected by the background due to the resistive loss around zero bias. In Fig. 11(a), we compare the values for $\rho_{+}$ and $\rho_{-}$ from such a numerical experiment, in comparison with the analytical result [Eq. (38)]. The agreement is very good and, as expected, even excellent for high-quality factors. The values $\rho_{+}$ and $\rho_{-}$ tend to $\rho_{0}/2$ when $\omega_{p} \to 0$, as the coefficients $\gamma_{+}$ and $\gamma_{-}$ [Eq. (7)] tend to $1/2$ for $\omega_{p} \to 0$ for an LC resonance matched with the matter excitation.

Our method thus relies on comparing LC resonators that are coupled and uncoupled to the quantum wells. Since the light-matter coupling in such systems depends only on the

electronic density in the quantum wells, we can compare a doped and undoped sample; note that the AlGaAs samples that we consider here are typically obtained from epitaxial growth, which produces identical heterostructure layers with high fidelity. Another possibility is applying a gate voltage in order to vary the charge density in the quantum-well region, which is a more direct method as it is applicable to the same device but increases the complexity of the device architecture.

Using Eq. (38) we can further extract the coefficients $\gamma_{\pm} = (\rho_{\pm}/\rho_{0})(\omega_{\pm}/\omega_{LC})^{3}$ as well as the charge quantum fluctuations enhancement factor [Eq. (8)]. We have the following identities:

$$ \gamma_{+} + \gamma_{-} = \frac{\rho_{+}\omega_{+}^{3} + \rho_{-}\omega_{-}^{3}}{\rho_{0}\omega_{LC}^{3}} = 1, \quad (62) $$

$$ \frac{\langle G | \hat{Q}^2 | G \rangle}{\langle 0 | \hat{Q}^2 | 0 \rangle} = \frac{\rho_{+} \omega_{+}^2 + \rho_{-} \omega_{-}^2}{\rho_{0} \omega_{LC}^2} = \omega_{LC} \frac{\rho_{+} \omega_{+}^2 + \rho_{-} \omega_{-}^2}{\rho_{+} \omega_{+}^3 + \rho_{-} \omega_{-}^3}. \quad (63) $$

Equation (62) can be considered as a normalization condition for the extracted values $\rho_{\pm}$, whereas Eq. (63) gives an estimate for the antisqueezing of the charge vacuum fluctuation in the new ground state. The second line in Eq. (63) is a direct consequence of Eq. (62). As seen in Fig. 11(b), the values estimated from the fits of the full model closely match the ideal case without dissipation. The antisqueezing effect becomes very significant in the case of USC and can be reliably demonstrated even for strong dissipation. The second line from Eq. (63) shows that the quantities $\rho_{\pm}$ and $\omega_{\pm}$ that are obtained from a transport experiment, together with the knowledge of the uncoupled frequencies $\omega_{LC}$, are sufficient to infer the absolute values of the ground-state charge fluctuations of the system. As shown in Fig. 11(c), this method has much higher precision, as the ratio that appears in the second line of Eq. (63) reduces strongly the effect of the uncertainties in the values $\rho_{\pm}$. Remarkably, with this expression both the low- and high-quality factor cases yield excellent estimates for the quantum vacuum fluctuations.
V. CONCLUSION

In the present work, we discussed the possibility to
probe quantum vacuum fluctuations, predicted by Ciuti et al.
in ultrastrongly coupled systems [3], in electronic trans-
port experiments. Such experiments can be enabled by THz
metamaterial architectures with nanoscopic capacitive parts
[46,57]. As shown in Sec. II A, in an LC circuit the antisqueez-
ing of charge vacuum fluctuations is another manifestation of
USC, along with the increasing ground-state photon occupa-
tion ⟨G|a†a|G⟩, which was noted in Ref. [3].

In our theoretical analysis, we employed first a purely
quantum model, which neglects all dissipation in the system.
Although nonrealistic, such a model has a conceptual strength
as it shows that the polaritonic states of the system and their
quantum correlations can be observed without any optical
absorption present. This is in stark contrast with spectroscopic
experiments, which rely on energy dissipation. The proposed
system thus allows for a quantum nondemolition measure-
ment of the polariton states, in the limit where the backaction
of the tunneling junction on the environment can be ignored.

A more realistic approach, including dissipation has been
considered in the framework of the P(E) theory [34]. The
latter can also be generalized at finite temperatures. To im-
plement this approach, we developed an equivalent circuit
mapping for the electronic excitation, which can be seen as
an LC circuit comprising a purely kinetic inductance. Such an
approach can clearly be useful for probing other systems, such
as localized surface plasmons in nanoparticles in tunneling
spectroscopy experiments [44]. The circuit approach is in
principle equivalent to the case in which the electromagnetic
environment of the tunnel junction can be described as an en-
semble of bosonic excitation. However, in our meta-material
architecture the capacitive parts can be reduced, in principle,
to arbitrarily small dimensions. In this limit, the fermionic
nature of the electronic excitation must be taken into account
[50]. To our knowledge, the DCB in such a mixed fermionic-
bosonic environment has never been envisioned before, and it
presents an exciting topic for future studies.

Whereas the present work is a theoretical study, all tech-
nological bricks exist separately, and it is possible to combine
them into an experimental device. In addition to the aforemen-
tioned metamaterial structures combined with semiconductor
quantum wells, tunnel junctions can also be fabricated from
semiconductor heterostructures [58]. The proposed system
offers a new paradigm for semiconductor optoelectronic de-
vices in the THz frequency range where the electromagnetic
quantum fluctuations regulate the electronic transport in the
device. An immediate perspective is a novel class of quantum
detectors that are sensitive to the quantum state of the incident
THz photons [21,35].

ACKNOWLEDGMENT

We acknowledge funding from the ERC-COG-863487
UNIQUE.

APPENDIX A: ALGEBRA OF THE
HOPFIELD COEFFICIENTS

Here, we recall briefly the analytical expressions of the
Hopfield coefficients that are obtained from the Bogoliubov

transform of the Hamiltonian Eq. (2). These results are very
general and can be used for any system of quantum harmonic
oscillators described by Eq. (2). It is convenient to set Ω<sub>R</sub> =
(ω<sub>P</sub>/2)√Fω<sub>LC</sub>/tilde{ω}<sub>21</sub>. The normal modes of the coupled system
are obtained using the definition of Eq. (5) and requiring
that ħωΠ = [Π, Ĥ] (the indexes “±” have been dropped for
convenience). Using the linear independence of the bosonic
operators, the following relations are obtained:

$$
H = \Omega_R (m - t), \quad (A1)
$$

$$
(\omega - \omega_{LC})x = (\omega + \omega_{LC})y = H, \quad (A2)
$$

$$
(\omega - \tilde{\omega}_{21})m = (\omega + \tilde{\omega}_{21})t = \Omega_R(x - y), \quad (A3)
$$

$$
x^2 - y^2 + m^2 - t^2 = 1. \tag{A4}
$$

Here the last equation arises from the condition
[Π, Π†] = 1. The above equations lead to the following re-
lationship:

$$
\frac{4H^2\omega\omega_{LC}}{(\omega^2 - \omega_{LC}^2)^2} \frac{2\omega^2 - \omega_{LC}^2 - \tilde{\omega}_{21}^2}{\omega^2 - \tilde{\omega}_{21}^2} = 1. \quad (A5)
$$

These relationships allow us to express the constant *H* and
then all the Hopfield coefficients explicitly as a function of the
frequencies *ω*<sub>LC</sub>, *tilde*ω<sub>21</sub>, *Ω*<sub>R</sub>. The frequencies *ω*<sub>±</sub> are provided
analytically as the roots of the bi-quadratic equation (4). These
roots satisfy the useful relation

$$
\omega_{+}^{2} + \omega_{-}^{2} = \omega_{LC}^{2} + \tilde{\omega}_{21}^{2}. \quad (A6)
$$

Next, the bosonic operators Π+, Π†+ must commute with
Π−, Π†−. As a consequence, we have the following matrix
relation between the Hopfield coefficients (see also Ref. [3]):

$$
\begin{bmatrix} x_+ & y_+ & m_+ & t_+ \\ y_+ & x_+ & t_+ & m_+ \\ x_- & y_- & m_- & t_- \\ y_- & x_- & t_- & m_- \end{bmatrix} \begin{bmatrix} x_+ & -y_+ & x_- & -y_- \\ -y_+ & x_+ & -y_- & x_- \\ m_+ & -t_+ & m_- & -t_- \\ -t_+ & m_+ & -t_- & m_- \end{bmatrix} = \mathbb{I}. \qquad (A7)
$$

The diagonal products express the commutation relations
[Π±, Π†±] = 1. The left matrix of Eq. (A7) is the one that links
the vector (Π+, Π†+, Π−, Π†−) to the vector (a, a†, p, p†). The
matrix on the right of Eq. (A7) then defines the inverse Bogoli-
ubov transform:

$$
a = x_{+} \Pi_{+} - y_{+} \Pi_{+}^{\dagger} + x_{-} \Pi_{-} - y_{-} \Pi_{-}^{\dagger}, \quad (A8)
$$

$$
a^{\dagger} = -y_{+} \Pi_{+} + x_{+} \Pi_{+}^{\dagger} - y_{-} \Pi_{-} + x_{-} \Pi_{-}^{\dagger}. \quad (A9)
$$

Similar expressions can be written for (p, p†). We are
thus in a position to express the ground-state average for
any observable in the system. For that, we use the fact that
the ground state is a tensor product |G⟩ = |G+⟩ ⊗ |G−⟩ with
Π±|G±⟩ = 0 and ⟨G±|Π±Π†±|G±⟩ = 1. For instance, in the
case ⟨G|a†a|G⟩ we obtain easily (see also Ref. [3])

$$
\langle G | a^{\dagger} a | G \rangle = y_{+}^{2} + y_{-}^{2}. \qquad (A10)
$$

Then from the above equations we can show that y² =
H²/(ω + ω_LC). Using Eqs. (A5) and (A6), we obtain the
expressions from Eqs. (6) and (7).
Let us now discuss the analytical derivation of the effective flux $\hat{\mathcal{F}}$ from Eq. (26) that describes the polaritonic environment of the tunnel junction. Combining Eqs. (34), (35), and (36) into Eq. (33), and extracting the flux $\hat{\mathcal{F}}$ from the definition in Eq. (25) with the help of the identity Eq. (30) applied for the polaritonic operators, we arrive at the following result:

$$ \hat{\mathcal{F}}_+ = i(\Pi_+ - \Pi_+^\dagger)\sqrt{\frac{L\hbar\omega_{LC}}{2}} \left\{ (x_++y_+)-(x_+-y_+)\frac{C_{eq}\omega_{LC}}{C\omega_+} - (m_+-t_+)\frac{1}{\omega_+}\sqrt{\frac{\omega_p^2}{C_pL\omega_{LC}\tilde{\omega}_{21}}} \right\}. \quad (A11) $$

Here we have divided the flux into an upper and a lower polariton contribution, $\hat{\mathcal{F}} = \hat{\mathcal{F}}_+ + \hat{\mathcal{F}}_-$, that have identical expressions except for the substitution $+ \to -$. Using the relations $\omega_{LC}^2 = 1/LC_{eq}$ and $F = C_{eq}/C_p$, the expression under the last square root is transformed as follows:

$$ \sqrt{\frac{\omega_p^2}{C_p L \omega_{LC} \tilde{\omega}_{21}}} = \sqrt{F \omega_p^2 \frac{\omega_{LC}}{\tilde{\omega}_{21}}} = 2\Omega_R. \quad (A12) $$

From Eq. (A3) we express $m_+$ and $t_+$ as a function of $x_+ - y_+$ to obtain

$$ \begin{aligned} (m_+ - t_+) \frac{2\Omega_R}{\omega_+} &= (x_+ - y_+) \frac{4\Omega_R^2 \tilde{\omega}_{21}}{\omega_+ (\omega_+^2 - \tilde{\omega}_{21}^2)} \\ &= \frac{F \omega_p^2 \omega_{LC}^2}{\omega_+ \omega_{LC} (\omega_+^2 - \tilde{\omega}_{21}^2)} \\ &= \frac{\omega_+^2 - \omega_{LC}^2}{\omega_+ \omega_{LC}}. \end{aligned} \quad (A13) $$

For the last line we used the fact that $\omega_+$ is a solution of the bi-quadratic equation (4). Next, we use Eq. (A2), which leads to the relation $(x_+ + y_+) = (\omega_+/\omega_{LC})(x_+ - y_+)$. Combining all the pieces in Eq. (A11), we arrive at Eq. (37) in the main text. Finally, multiplying the autocorrelation function $\langle \hat{\mathcal{F}}_+^2 \rangle$ by $e^2/\hbar^2$, the prefactor is transformed as

$$ \begin{aligned} \frac{e^2}{\hbar} \frac{L \omega_{LC}}{2(1+\eta)^2} \gamma_+ \left( \frac{\omega_{LC}}{\omega_+} \right)^3 &= \frac{e^2}{\hbar} \frac{1}{2\omega_{LC}C_{eq}(1+\eta)^2} \gamma_+ \left( \frac{\omega_{LC}}{\omega_+} \right)^3 \\ &= \frac{e^2}{2\hbar\omega_{LC}C_T(1+\eta)} \gamma_+ \left( \frac{\omega_{LC}}{\omega_+} \right)^3. \end{aligned} \quad (A14) $$

This is exactly the expression of $\rho_+$, Eq. (38), from the main text.

## APPENDIX B: FEASIBILITY OF DCB IN THz METAMATERIALS

Here we illustrate a terahertz metamaterial resonator architecture where the DCB can be observed. This architecture is based on Refs. [46,57]. The exact geometry for a geometrical realization might be conditioned by various factors, such as an adequate choice of a fabrication procedure, or the necessity to isolate electrically the quantum-well stack from the tunnel junction. Beyond the two-capacitor structure from Fig. 2, one

<figure><img src="image_15.png" /><figcaption>FIG. 12. (a) A typical geometry for a metamaterial resonator used for our simulations. (b), (c) Illustration of the field components $E_z$ (b) and $H_z$ (c) for an $xy$ plane placed in the middle of the capacitors. (d) Side view of the geometry and a definition of a coordinate system. (e) Closeup view of one of the capacitors showing the $E_z$ field. The arrows indicate the direction of the field.</figcaption></figure>

can envision resonators with multiple capacitor pairs, etc. Therefore, rather than performing a systematic assessment of the DCB as a function of the metamaterial geometry, our aim here is to show that typical values of $\rho_0 \approx 0.3-0.5$ can be achieved with the meta-atoms concept that has already been reported in the literature [57,64].

A typical geometry, similar to the one discussed in the main text, is shown in Fig. 12, together with the corresponding simulations of the magnetic and electric fields. Commercial finite-difference domain software is used in order to obtain the electromagnetic eigenmodes of the structure. The dielectric filling of the capacitors is considered here to be gallium arsenide (GaAs) with a dielectric constant $\varepsilon = 12.25$. With the geometric parameters indicated in Fig. 12(a), the fundamental resonance is obtained to be at 1.8 THz. In Figs. 12(b) and 12(c) we illustrate, respectively, the $z$-component of the electric and magnetic fields of this resonance for an $xy$-plane that passes through the middle of the capacitors. The coordinate system is indicated in Fig. 12(d) together with a side view of the system. From Fig. 12(b) it is clear that the electric field is strongly localized in the capacitors, as anticipated, however there is clearly a leakage of the electric energy density along the inductive part of the resonator. This is due to propagation effects, as discussed in Refs. [57,64]. Therefore, the effective capacitance of the resonator cannot be obtained only from the geometrical capacitance of the double-metal region, but the whole structure should be considered.

We can see from Figs. 12(b) and 12(c) that the electric field distribution is asymmetric with respect to an $xz$-plane passing through the middle of the resonator. To estimate the capacitance of the structure, we compute numerically the total electric energy:

$$ W_e = \frac{1}{2} \iiint \varepsilon \varepsilon_0(\mathbf{r}) \mathbf{E}^2 d^3\mathbf{r} \quad (B1) $$
<figure><img src="image_16.png" /><figcaption>FIG. 13. (a) Indication of the charges induced in the top $Q_2$ and bottom $Q_1$ parts of the resonator. (b) Equivalent circuit which accounts for $Q_2 > Q_1$ by introducing an additional capacitor $C_{inter}$.</figcaption></figure>

as well as the charges induced onto the metallic parts. For the latter we use the expression $\sigma = \varepsilon\varepsilon_0 E_n$, which links the surface charge density $\sigma$ to the normal component of the electric field $E_n$, which has been illustrated in Fig. 12(e). The total charges on the top $Q_2$ and bottom $Q_1$ surfaces are then

$$Q_{1,2} = \iint_{1,2} \varepsilon\varepsilon_0(\mathbf{r})\mathbf{E} \cdot d\mathbf{S}. \qquad (B2)$$

The integral runs either along the half-part of the bottom metal part (“1”) or the top metal part (“2”); see Fig. 13(a). We obtain the numerical values $W_e = 3.3 \times 10^{-13}$ J, $Q_1 = 11 \times 10^{-16}$ C, and $Q_2 = 50 \times 10^{-16}$ C. The fact that $Q_2$ is larger than $Q_1$ indicates that there is an additional electric energy on the parallel sections of the inductance, as described in Ref. [64]. Indeed, these sections can be modeled as a short portion of a transmission line [64]. To take into account this effect, we introduce a capacitance $C_{inter}$ as described in the equivalent circuit shown in Fig. 13(b). The circuit analysis provides the following relations:

$$W_e = \frac{Q_1^2}{C_T} + \frac{(Q_2 - Q_1)^2}{2C_{inter}}, \qquad (B3)$$

$$\frac{Q_2 - Q_1}{C_{inter}} = \frac{2Q_1}{C_T}. \qquad (B4)$$

From these equations we obtain numerical values $C_T = 17$ aF and $C_{inter} = 2.5$ aF. As an example, let us consider the situation in which there is a tunnel junction inserted into one of the capacitors $C_T$. By considering the equivalent circuit seen by the junction, the DCB parameter can be estimated to be

$$\rho = \frac{e^2}{2C_T \hbar \omega_{LC}} \left( \frac{1}{2} + \frac{C_{inter}}{C_T} \right). \qquad (B5)$$

The THz resonance corresponds to $\hbar\omega_{LC} = 7.4$ meV, while the charging energy of a single electron is found to be 3 meV. A value $\rho \approx 0.4$ is obtained, while other simulations provide values in the range 0.3–0.5 as a function of the device geometry. Such meta-atom architectures thus have strong potential for the observation of DCB high ratios $\rho$.

Let us further discuss the feasibility of electrical measurements with such structures. The devices must operate typically in the voltage range between 0 and 20 meV. We have to respect several constraints: (i) keeping a very small

<figure><img src="image_17.png" /><figcaption>FIG. 14. (a) Schematics of the conduction-band profile for a biased semiconductor tunnel junction based on a GaAlAs/GaAs heterostructure. (b) Current-voltage characteristics through a junction with an area $100 \times 100$ nm$^2$ and various thicknesses $W_B$. The potential height is $\phi_B = 80$ meV. (c) Differential conductance normalized on the conduction quantum for the three junctions shown in (b).</figcaption></figure>

capacitor area, on the order of $10^4$ nm$^2$, (ii) the overall current through the system should be measurable, and (iii) the differential conductance should be kept below the value of the conductance quantum $G_K$. The tunnel junction must therefore be relatively transparent: such junctions can be realized from semiconductor GaAs/AlGaAs heterostructures, as illustrated in Fig. 14(a). The current voltage characteristics of the junction is provided by the Simmons formula [83]:

$$I(U) = \frac{W^2 G_K}{W_B^2 2\pi e} \left[ (\phi_B - eU/2)e^{-2\sqrt{\frac{\phi_B - eU/2}{E_B}}} - (\phi_B + eU/2)e^{-2\sqrt{\frac{\phi_B + eU/2}{E_B}}} \right], \qquad (B6)$$

$$E_B = \frac{\hbar^2}{2m^* W_B^2}. \qquad (B7)$$

Here we consider the case $\phi_B = 80$ meV, which corresponds to a Ga$_{0.85}$Al$_{0.15}$As barrier. The barrier height thus exceeds several times the applied bias ($\approx 20$ meV) as required in our studies. In Eq. (B6), $W^2$ is the electrical area
of the junction and $m^*$ is the effective electron mass. Results of Eq. (B6) are plotted in Fig. 14(b) for various values of the barrier thickness $W_B$, and the corresponding curves for the differential conductance $dI/dU$ normalized on conductance quantum $G_K$ are plotted in Fig. 14(c). We recall that the condition $(1/G_K)dI/dU < 1$ is required in order to observe the DCB effect. As seen from Fig. 14(c), this condition is

well satisfied. At the same time, by an appropriate choice of the thickness $W_B$ we can find the values of the tunneling current on the order of $nA$ that are measurable despite the very small area of the junction $W^2$. Such heterostructures thus have a sufficient number of degrees of freedom [58] in order to allow experimental observation of the DCB even with a single metamaterial resonator.

[1] E. Purcell, Spontaneous emission probabilities at radio frequencies, Phys. Rev. **69**, 674 (1946).

[2] M. Fox, *Quantum Optics: An Introduction*, Oxford Master Series in Atomic, Optical, and Laser Physics (Oxford University Press, Oxford, 2006).

[3] C. Ciuti, G. Bastard, and I. Carusotto, Quantum vacuum properties of the intersubband cavity polariton field, Phys. Rev. B **72**, 115303 (2005).

[4] P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Ultrastrong coupling regimes of light-matter interaction, Rev. Mod. Phys. **91**, 025005 (2019).

[5] A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, Ultrastrong coupling between light and matter, Nat. Rev. Phys. **1**, 19 (2019).

[6] P. Nataf, T. Champel, G. Blatter, and D. M. Basko, Rashba cavity QED: A route towards the superradiant quantum phase transition, Phys. Rev. Lett. **123**, 207402 (2019).

[7] D. Guerci, P. Simon, and C. Mora, Superradiant phase transition in electronic systems and emergent topological phases, Phys. Rev. Lett. **125**, 257604 (2020).

[8] F. Schlawin, D. M. Kennes, and M. A. Sentef, Cavity quantum materials, Appl. Phys. Rev. **9**, 011312 (2022).

[9] M. A. Sentef, M. Ruggenthaler, and A. Rubio, Cavity quantum-electrodynamical polaritonically enhanced electron-phonon coupling and its influence on superconductivity, Sci. Adv. **4**, eaau6969 (2018).

[10] F. Schlawin, A. Cavalleri, and D. Jaksch, Cavity-mediated electron-photon superconductivity, Phys. Rev. Lett. **122**, 133602 (2019).

[11] M. Artoni and J. L. Birman, Quantum-optical properties of polariton waves, Phys. Rev. B **44**, 3736 (1991).

[12] M. Artoni and J. L. Birman, Detection of optical squeezing and photon statistics in polaritons, Opt. Commun. **89**, 324 (1992).

[13] S. De Liberato, C. Ciuti, and I. Carusotto, Quantum vacuum radiation spectra from a semiconductor microcavity with a time-modulated vacuum Rabi frequency, Phys. Rev. Lett. **98**, 103602 (2007).

[14] C. M. Wilson, G. Johansson, A. Pourkabirian, M. Simoen, J. R. Johansson, T. Duty, F. Nori, and P. Delsing, Observation of the dynamical Casimir effect in a superconducting circuit, Nature (London) **479**, 376 (2011).

[15] V. V. Dodonov, Current status of the dynamical casimir effect, Phys. Scr. **82**, 038105 (2010).

[16] S. W. Hawking, Particle creation by black holes, Commun. Math. Phys. **43**, 199 (1975).

[17] R. Loudon, *The Quantum Theory of Light* (Oxford University Press, Oxford, 2000).

[18] M. H. Devoret, Quantum fluctuations in electrical circuits, in *Fluctuations Quantiques/Quantum Fluctuations: Les Houches Session LXIII*, 27 June-28 July 1995, edited by S. Reynaud, E. Giacobino, and J. Zinn-Justin (Elsevier, Amsterdam, The Netherlands, 1997), p. 351.

[19] A. Cottet, T. Kontos, and B. Douçot, Electron-photon coupling in mesoscopic quantum electrodynamics, Phys. Rev. B **91**, 205417 (2015).

[20] G. L. Paravicini-Bagliani, F. Appugliese, E. Richter, F. Valmorra, J. Keller, M. Beck, N. Bartolo, C. Rössler, T. Ihn, K. Ensslin, C. Ciuti, G. Scalari, and J. Faist, Magneto-transport controlled by Landau polariton states, Nat. Phys. **15**, 186 (2019).

[21] F. Valmorra, K. Yoshida, L. Contamin, S. Messelot, S. Massabeau, M. Delbecq, M. Dartailh, M. Desjardins, T. Cubaynes, Z. Leghtas, K. Hirakawa, J. Tignon, S. Dhillon, S. Balibar, J. Mangeney, A. Cottet, and T. Kontos, Vacuum-field-induced THz transport gap in a carbon nanotube quantum dot, Nat. Commun. **12**, 5490 (2021).

[22] F. Appugliese, J. Enkner, G. L. Paravicini-Bagliani, M. Beck, C. Reichl, W. Wegscheider, G. Scalari, C. Ciuti, and J. Faist, Breakdown of topological protection by cavity vacuum fields in the integer quantum Hall effect, Science **375**, 1030 (2022).

[23] K. Kuroyama, J. Kwoen, Y. Arakawa, and K. Hirakawa, Coherent interaction of a few-electron quantum dot with a terahertz optical resonator, Phys. Rev. Lett. **132**, 066901 (2024).

[24] J. Feist and F. J. Garcia-Vidal, Extraordinary exciton conductance induced by strong coupling, Phys. Rev. Lett. **114**, 196402 (2015).

[25] D. Hagenmüller, J. Schachenmayer, S. Schütz, C. Genes, and G. Pupillo, Cavity-enhanced transport of charge, Phys. Rev. Lett. **119**, 223601 (2017).

[26] D. Hagenmüller, S. Schütz, J. Schachenmayer, C. Genes, and G. Pupillo, Cavity-assisted mesoscopic transport of fermions: Coherent and dissipative dynamics, Phys. Rev. B **97**, 205303 (2018).

[27] N. Bartolo and C. Ciuti, Vacuum-dressed cavity magnetotransport of a two-dimensional electron gas, Phys. Rev. B **98**, 205301 (2018).

[28] E. Orgiu, J. George, J. A. Hutchison, E. Devaux, J. F. Dayen, B. Doudin, F. Stellacci, C. Genet, J. Schachenmayer, C. Genes, G. Pupillo, P. Samorì, and T. W. Ebbesen, Conductivity in organic semiconductors hybridized with the vacuum field, Nat. Mater. **14**, 1123 (2015).

[29] T. W. Ebbesen, Hybrid light-matter states in a molecular and material science perspective, Acc. Chem. Res. **49**, 2403 (2016).

[30] M. Lagrée, M. Jeannin, G. Quinchard, O. Ouznali, A. Evirgen, V. Trinité, R. Colombelli, and A. Delga, Direct polariton-to-
electron tunneling in quantum cascade detectors operating in the strong light-matter coupling regime, Phys. Rev. Appl. **17**, 044021 (2022).

[31] F. Pisani, D. Gacemi, A. Vasanelli, L. Li, A. Davies, E. Linfield, C. Sirtori, and Y. Todorov, Electronic transport driven by collective light-matter coupled states in a quantum device, Nat. Commun. **14**, 3914 (2023).

[32] M. H. Devoret, D. Esteve, H. Grabert, G.-L. Ingold, H. Pothier, and C. Urbina, Effect of the electromagnetic environment on the coulomb blockade in ultrasmall tunnel junctions, Phys. Rev. Lett. **64**, 1824 (1990).

[33] S. M. Girvin, L. I. Glazman, M. Jonson, D. R. Penn, and M. D. Stiles, Quantum fluctuations and the single-junction Coulomb blockade, Phys. Rev. Lett. **64**, 3183 (1990).

[34] G.-L. Ingold and Y. V. Nazarov, in *Single Charge Tunneling*, edited by H. Grabert and M. H. Devoret, of NATO ASI Series B: Physics Vol. 294 (Plenum, New York, 1992), pp. 21–107.

[35] J. R. Souquet, M. J. Woolley, J. Gabelli, P. Simon, and A. A. Clerk, Photon-assisted tunnelling with nonclassical light, Nat. Commun. **5**, 5562 (2014).

[36] P. Delsing, K. K. Likharev, L. S. Kuzmin, and T. Claeson, Effect of high-frequency electrodynamic environment on the single-electron tunneling in ultrasmall junctions, Phys. Rev. Lett. **63**, 1180 (1989).

[37] A. N. Cleland, J. M. Schmidt, and J. Clarke, Influence of the environment on the Coulomb blockade in submicrometer normal-metal tunnel junctions, Phys. Rev. B **45**, 2950 (1992).

[38] T. Holst, D. Esteve, C. Urbina, and M. H. Devoret, Effect of a transmission line resonator on a small capacitance tunnel junction, Phys. Rev. Lett. **73**, 3455 (1994).

[39] M. Hofheinz, F. Portier, Q. Baudouin, P. Joyez, D. Vion, P. Bertet, P. Roche, and D. Esteve, Bright side of the Coulomb blockade, Phys. Rev. Lett. **106**, 217005 (2011).

[40] C. Altimiras, O. Parlavecchio, P. Joyez, D. Vion, P. Roche, D. Esteve, and F. Portier, Dynamical Coulomb blockade of shot noise, Phys. Rev. Lett. **112**, 236803 (2014).

[41] F. D. Parmentier, A. Anthore, S. Jezouin, H. le Sueur, U. Gennser, A. Cavanna, D. Mailly, and F. Pierre, Strong back-action of a linear circuit on a single electronic quantum channel, Nat. Phys. **7**, 935 (2011).

[42] A. Chepelianskii, P. Delplace, A. Shailos, A. Kasumov, R. Deblock, M. Monteverde, C. Ojeda-Aristizabal, M. Ferrier, S. Guéron, and H. Bouchiat, Phonon-assisted dynamical Coulomb blockade in a thin suspended graphite sheet, Phys. Rev. B **79**, 235418 (2009).

[43] C. Brun, K. H. Müller, I.-P. Hong, F. Patthey, C. Flindt, and W.-D. Schneider, Dynamical Coulomb blockade observed in nanosized electrical contacts, Phys. Rev. Lett. **108**, 126802 (2012).

[44] J. Senkpiel, J. C. Klöckner, M. Etzkorn, S. Dambach, B. Kubala, W. Belzig, A. L. Yeyati, J. C. Cuevas, F. Pauly, J. Ankerhold, C. R. Ast, and K. Kern, Dynamical Coulomb blockade as a local probe for quantum transport, Phys. Rev. Lett. **124**, 156803 (2020).

[45] Y. Todorov, A. M. Andrews, R. Colombelli, S. De Liberato, C. Ciuti, P. Klang, G. Strasser, and C. Sirtori, Ultrastrong light-matter coupling regime with polariton dots, Phys. Rev. Lett. **105**, 196402 (2010).

[46] M. Jeannin, G. Mariotti Nesurini, S. Suffit, D. Gacemi, A. Vasanelli, L. Li, A. G. Davies, E. Linfield, C. Sirtori,

and Y. Todorov, Ultrastrong light-matter coupling in deeply subwavelength THz LC resonators, ACS Photon. **6**, 1207 (2019).

[47] M. Geiser, C. Walther, G. Scalari, M. Beck, M. Fischer, L. Nevou, and J. Faist, Strong light-matter coupling at terahertz frequencies at room temperature in electronic LC resonators, Appl. Phys. Lett. **97**, 191107 (2010).

[48] E. Strupiechonski, G. Xu, M. Brekenfeld, Y. Todorov, N. Isac, A. M. Andrews, P. Klang, C. Sirtori, G. Strasser, A. Degiron, and R. Colombelli, Sub-diffraction-limit semiconductor resonators operating on the fundamental magnetic resonance, Appl. Phys. Lett. **100**, 131113 (2012).

[49] Y. Todorov and C. Sirtori, Intersubband polaritons in the electrical dipole gauge, Phys. Rev. B **85**, 045304 (2012).

[50] Y. Todorov and C. Sirtori, Few-electron ultrastrong light-matter coupling in a quantum LC circuit, Phys. Rev. X **4**, 041031 (2014).

[51] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, *Photons et Atomes* (EDP Sciences/CNRS Editions, Paris, 2001).

[52] J. J. Hopfield, Theory of the contribution of excitons to the complex dielectric constant of crystals, Phys. Rev. **112**, 1555 (1958).

[53] M. Helm, *Intersubband Transitions in Quantum Wells*, edited by H. C. Liu and F. Capasso (Academic, San Diego, 2000).

[54] C. Mora, C. Altimiras, P. Joyez, and F. Portier, Quantum properties of the radiation emitted by a conductor in the Coulomb blockade regime, Phys. Rev. B **95**, 125311 (2017).

[55] Y. V. Nazarov, Anomalous current-voltage characteristics of tunnel junctions, Sov. Phys. JETP **68**, 561 (1989).

[56] M. Jeannin, T. Bonazzi, D. Gacemi, A. Vasanelli, S. Suffit, L. Li, A. G. Davies, E. Linfield, C. Sirtori, and Y. Todorov, High temperature metamaterial terahertz quantum detector, Appl. Phys. Lett. **117**, 251102 (2020).

[57] A. Mottaghizadeh, Y. Todorov, M. Cameau, D. Gacemi, A. Vasanelli, and C. Sirtori, Nanoscale electromagnetic confinement in THz circuit resonators, Opt. Express **25**, 28718 (2017).

[58] S. L. Feng, J. Krynicki, M. Zazoui, J. C. Bourgoin, P. Bois, and E. Rosencher, Electron transport through GaAlAs barriers in GaAs, J. Appl. Phys. **74**, 341 (1993).

[59] Y. Todorov, L. Tosetto, J. Teissier, A. M. Andrews, P. Klang, R. Colombelli, I. Sagnes, G. Strasser, and C. Sirtori, Optical properties of metal-dielectric-metal microcavities in the THz frequency range, Opt. Express **18**, 13886 (2010).

[60] R. Singh, C. Rockstuhl, and W. Zhang, Strong influence of packing density in terahertz metamaterials, Appl. Phys. Lett. **97**, 241108 (2010).

[61] S. Datta, *Electronic Transport in Mesoscopic Systems* (Cambridge University Press, Cambridge, 1995).

[62] L. Wendler and E. Kändler, Intra- and intersubband plasmon-polaritons in semiconductor quantum wells, Phys. Status Solidi B **177**, 9 (1993).

[63] S. Luryi, Quantum capacitance devices, Appl. Phys. Lett. **52**, 501 (1988).

[64] M. Jeannin, D. Gacemi, A. Vasanelli, L. Li, A. G. Davies, E. Linfield, G. Biasol, C. Sirtori, and Y. Todorov, Quasi-static and propagating modes in three-dimensional THz circuits, Opt. Express **28**, 16982 (2020).

[65] G.-L. Ingold, H. Grabert, and U. Eberhardt, Cooper-pair current through ultrasmall Josephson junctions, Phys. Rev. B **50**, 395 (1994).
[66] H. Grabert, G.-L. Ingold, M. H. Devoret, D. Estève, H. Pothier, and C. Urbina, Single electron tunneling rates in multijunction circuits, Z. Phys. B **84**, 143 (1991).

[67] S. Rajabali, S. Markmann, E. Jö, M. Beck, C. A. Lehner, W. Wegscheider, J. Faist, and G. Scalari, An ultrastrongly coupled single terahertz meta-atom, Nat. Commun. **13**, 2528 (2022).

[68] A. Caldeira and A. Leggett, Quantum tunnelling in a dissipative system, Ann. Phys. **149**, 374 (1983).

[69] M. Jeannin, T. Bonazzi, D. Gacemi, A. Vasanelli, L. Li, A. G. Davies, E. Linfield, C. Sirtori, and Y. Todorov, Absorption engineering in an ultrasubwavelength quantum system, Nano Lett. **20**, 4430 (2020).

[70] E. Rodriguez, T. Bonazzi, H. Dely, M. Mastrangelo, K. Pantzas, G. Beaudoin, I. Sagnes, A. Vasanelli, Y. Todorov, and C. Sirtori, Metamaterial engineering for optimized photon absorption in unipolar quantum devices, Opt. Express **30**, 20515 (2022).

[71] J. B. Khurgin and A. Boltasseva, Reflecting upon the losses in plasmonics and metamaterials, MRS Bull. **37**, 768 (2012).

[72] G. W. Ford, J. T. Lewis, and R. F. O'Connell, Quantum langevin equation, Phys. Rev. A **37**, 4419 (1988).

[73] W. Louisell, *Quantum Statistical Properties of Radiation*, A Wiley-Interscience Publication (Wiley, New York, 1973).

[74] C. Ciuti and I. Carusotto, Input-output theory of cavities in the ultrastrong coupling regime: The case of time-independent cavity parameters, Phys. Rev. A **74**, 033811 (2006).

[75] H. A. Haus, *Waves and Fields in Optoelectronics* (Prentice-Hall, Suffolk, 1983).

[76] W. Kohn, Cyclotron resonance and de Haas-van Alphen oscillations of an interacting electron gas, Phys. Rev. **123**, 1242 (1961).

[77] A. Delteil, A. Vasanelli, Y. Todorov, C. Feuillet Palma, M. Renaudat St-Jean, G. Beaudoin, I. Sagnes, and C. Sirtori, Charge-induced coherence between intersubband plasmons in a quantum structure, Phys. Rev. Lett. **109**, 246808 (2012).

[78] M. Geiser, F. Castellano, G. Scalari, M. Beck, L. Nevou, and J. Faist, Ultrastrong coupling regime and plasmon polaritons in parabolic semiconductor quantum wells, Phys. Rev. Lett. **108**, 106402 (2012).

[79] P. Goulain, C. Deimert, M. Jeannin, S. Piotta, W. J. Pasek, Z. Wasilewski, R. Colombelli, and J.-M. Manceau, THz ultra-strong light-matter coupling up to 200 K with continuously-graded parabolic quantum wells, Adv. Opt. Mater. **11**, 2202724 (2023).

[80] C. Emary and T. Brandes, Chaos and the quantum phase transition in the Dicke model, Phys. Rev. E **67**, 066203 (2003).

[81] P. Nataf and C. Ciuti, No-go theorem for superradiant quantum phase transitions in cavity QED and counter-example in circuit QED, Nat. Commun. **1**, 72 (2010).

[82] C. Altimiras, U. Gennser, A. Cavanna, D. Mailly, and F. Pierre, Experimental test of the dynamical Coulomb blockade theory for short coherent conductors, Phys. Rev. Lett. **99**, 256805 (2007).

[83] J. G. Simmons, Generalized formula for the electric tunnel effect between similar electrodes separated by a thin insulating film, J. Appl. Phys. **34**, 1793 (1963).