Probing the diamagnetic term in light–matter interaction

The ultra-strong coupling (USC) regime of light–matter interaction has recently attracted much interest, thanks to impressive advances in theory and experiments [1, 2]. A number of interesting phenomena, of fundamental physical appeal and potential impact on future quantum technologies, arise in this regime, including the formation of a nontrivial ground state [3], dynamical Casimir-like effects [4], and the possibility of quantum phase transitions [5, 6].

Loosely speaking, the USC regime is entered when the light–matter coupling constant, say λ, is a non-negligible fraction of the bare frequencies of light and matter. In this situation the rotating-wave approximation breaks down, and one cannot neglect the so-called 'counter-rotating terms' in the interaction Hamiltonian [3, 7]. While this brings about exciting new physics, it results also in formidable computational challenges. Effective models featuring a small number of degrees-of-freedom (modes), such as the Dicke Hamiltonian [5], are extremely useful in this context: they can reveal crucial features of the new regime while keeping computations tractable.

However, there is disagreement in the literature regarding the specific form that these few-degrees-of-freedom models should take, in particular concerning the presence (or otherwise) of the so-called diamagnetic or A² term. We recall that such a term originates from the general form of the Coulomb gauge minimal coupling Hamiltonian

$$\begin{eqnarray}&&{{ \mathcal H }}_{\min }=\displaystyle \sum _{j}\displaystyle \frac{{({\hat{{\bf{p}}}}_{j}-{q}_{j}\hat{{\bf{A}}})}^{2}}{2{m}_{j}}+V({\hat{{\bf{x}}}}_{1},{\hat{{\bf{x}}}}_{2},\ldots )+{{ \mathcal H }}_{\mathrm{EM}},\end{eqnarray} \tag{ 1 }$$

which describes the interaction of the quantized electromagnetic field with nonrelativistic charged particles. In the above equation m_j, q_j, ${\hat{{\bf{x}}}}_{j}$ and ${\hat{{\bf{p}}}}_{j}$ indicate respectively the mass, charge, position and canonical momentum of the jth particle, $\hat{{\bf{A}}}$ is the vector potential operator, V is the instantaneous Coulomb potential energy, while ${{ \mathcal H }}_{\mathrm{EM}}$ is the free Hamiltonian of the transverse radiation field [7]. ${{ \mathcal H }}_{\min }$ clearly contains a contribution proportional to the squared vector potential ${\hat{{\bf{A}}}}^{2}$, i.e. the diamagnetic term.

In cavity QED (and other similar set-ups) many authors derive their effective models by a direct few-mode truncation of Hamiltonian (1), thus explicitly retaining an A²–like term affecting the radiation modes of interest. The importance of such term is well recognized in preventing the Dicke phase transition [8–11] and in limiting the effectively achievable light–matter coupling [12]. Within such framework, the coupling strength associated with the A² term obeys precise constraints (see below), whose origin can be traced back to the Thomas–Reiche–Kuhn (TRK) sum rule [5].

In contrast, other authors hold that the few-mode approximation is better carried out in the electric dipole gauge, that is, they apply a canonical transformation to Hamiltonian (1) before the derivation of an approximate model. This yields an effective Hamiltonian without diamagnetic term, and seemingly capable of a phase transition [13–15]. Yet, it has been argued that such derivations neglect the contribution of the so-called transverse P² term, which in the electric dipole gauge plays a role mathematically analogue to A² [10].

In the context of circuit QED, [5] maintains that the diamagnetic term can be avoided by appropriately tuning the circuit parameters, a claim which has generated further debate [16–19]. In fact, in circuit QED, there is even disagreement on the most appropriate microscopic description underlying the Dicke model. It is indeed not clear whether the most reliable starting point for the derivation of an effective model should be a minimal coupling Hamiltonian [18] or a direct quantization of macroscopic circuit variables [19].

In this letter we take a complementary approach to the A²-term debate: instead of trying to address the question theoretically, by deriving a Dicke model from first principles, we propose that it may be settled experimentally via an appropriate measurement. We note that steps in this direction have been taken in recent theoretical proposals [20, 21]. Moreover, present-day experimental results seem to support the thesis of Nataf and Ciuti [5], i.e. that the A² term may be neglected in some circuit QED systems [22]. In this work we contribute to these efforts by exploiting the tools of quantum metrology. We address the estimation of the A² coupling constant via measurements on the coupled ground state, thus deriving a theoretical benchmark relevant to any experimental implementation of the model (circuit QED and cavity QED being two relevant special cases). Our work is a first step towards setting precise metrological standards for any experiment probing the diamagnetic parameter.

We consider a general Dicke model of two coupled oscillators featuring counter-rotating terms and an A²-like term with unknown coupling constant D. As the establishment of a nontrivial ground state is perhaps the simplest signature of the USC, we find it natural to ask how much information about D is contained in this state (i.e., how sensitive the state is to variations in D), and how efficiently such information can be extracted by typical measurements (see figure 1 for an illustration). Exploiting the tools of quantum estimation theory [23–28], we are able to give quantitative answers to these questions.

Zoom In Zoom Out Reset image size

Besides deriving the ultimate bounds on the estimation precision of D, we find that typical measurements on the transverse cavity field, such as homodyne detection or photon counting, allow for the precise estimation of the diamagnetic term in a wide range of model parameters. Namely, homodyne detection is optimal when the coupling constant λ is vanishing, still providing a close-to-optimal performance for $\lambda \lesssim 20 \%$ of the bare mode frequencies. This regime of relatively small coupling (yet within the USC) is relevant in both cavity QED [29] and circuit QED experiments [30, 31], and is well described within the two-mode approximation [32]. Moreover, if the Dicke model of interest describes a collection of N two-level emitters, this is also a regime in which low-lying excitations are well described by a coupled oscillator model, even for relatively small values of N: for example, this is the case for a spin-2 Dicke model (N = 5) [32].

While we tackle the problem in its full generality, we also obtain compact analytical results for three important limiting cases: (i) $D\simeq 0;$ (ii) $D\simeq {D}_{\mathrm{TRK}}$, where ${D}_{\mathrm{TRK}}$ is the value derived from the TRK sum rule [5, 32]; (iii) $D\simeq {D}_{\mathrm{crit}}$, where ${D}_{\mathrm{crit}}$ is the threshold for the Dicke phase transition. In case the model admits a phase transition, we further find that the considered measurements become approximately optimal near the critical point. However we point out that finite-size effects, not included in our model, may become increasingly relevant as criticality is approached [33–35].

We complement our analysis by discussing the problem of binary discrimination [36, 37] between the choices D = 0 and $D={D}_{\mathrm{TRK}}$, which may be relevant in circuit QED [5, 16]. We then discuss the experimental feasibility of the proposed measurements, and conclude the letter by providing an outlook of further research questions stemming from our work.

Before proceeding we note that in this letter we adopt a language appropriate to the Coulomb gauge, such that field operators describe excitations of the transverse (and gauge invariant) radiation field, while what we call 'matter mode' in fact describes a combination of matter and field properties [7]. In this form, our results can be used to determine the most appropriate coupling constant of the diamagnetic term. In appendix C we show that our results can be readily adapted to describe the estimation of the P² term in dipole gauge [10].

Let us start by introducing the Hamiltonian of the system. We denote with a the photonic cavity mode, of bare frequency ${\omega }_{a}$, and with b the matter mode, with frequency ${\omega }_{b}$. After setting ${\hslash }=1$, the Hamiltonian reads

$$\begin{eqnarray}&&{ \mathcal H }={\omega }_{a}{a}^{\dagger }a+{\omega }_{b}{b}^{\dagger }b+\lambda (a+{a}^{\dagger })(b+{b}^{\dagger })+D{(a+{a}^{\dagger })}^{2},\end{eqnarray} \tag{ 2 }$$

where λ is the light–matter coupling strength and the last is the diamagnetic term, with the constant D being the object of our discussion.

Being quadratic in the field operators, ${ \mathcal H }$ is easily diagonalized and can be written, up to an irrelevant constant, as ${ \mathcal H }={\omega }_{U}{p}_{U}^{\dagger }{p}_{U}+{\omega }_{L}{p}_{L}^{\dagger }{p}_{L},$ where ${p}_{U,L}$ are bosonic operators for the upper and lower polaritonic modes. The eigenfrequencies ${\omega }_{U,L}$ are reported in [32] and in appendix A.

If the diamagnetic term is absent, i.e. we set D = 0 in equation (2), we recover the standard Dicke Hamiltonian, with a quantum phase transition at $\lambda ={\lambda }_{\mathrm{crit}}\equiv \sqrt{{\omega }_{a}{\omega }_{b}}/2$. On the other hand, if the diamagnetic term is included, it must satisfy the inequality

$$\begin{eqnarray}&&D\geqslant {D}_{\mathrm{crit}}\equiv {\lambda }^{2}/{\omega }_{b}-{\omega }_{a}/4\end{eqnarray} \tag{ 3 }$$

in order for the Hamiltonian in equation (2) to be bounded from below. If D assumes the 'TRK value' [32],

$$\begin{eqnarray}&&D={D}_{\mathrm{TRK}}\equiv {\lambda }^{2}/{\omega }_{b},\end{eqnarray} \tag{ 4 }$$

then equation (3) is always satisfied and the phase transition is suppressed.

The ground state of the Hamiltonian is the vacuum state of the polaritonic modes, i.e. a Gaussian state with zero mean and covariance matrix ${\sigma }_{0}={\mathbb{I}}/2$. Expressing this state in terms of the initial modes a and b, we obtain a nontrivial, two-mode squeezed Gaussian state with covariance matrix $\sigma =S{\sigma }_{0}{S}^{T}$, where the suitable symplectic transformation S is reported in appendix A.

Let us now address the problem of estimating D, and in doing so we shall review the basics of quantum estimation theory. In order to determine the value of the unknown parameter, one collects measurement outcomes from experiments and builds an estimator $\hat{D}$, i.e. a function of the data that will return our best guess for the value of D. Assuming that the estimator is unbiased, its precision, in the limit of a large number n of measurements, is bounded from below by the Cramér–Rao bound [38] $\mathrm{Var}\hat{D}\geqslant {[{nF}(D)]}^{-1}$, where F(D) is the Fisher information (FI) of the probability distribution $p(x| D)$, defined as $F(D)=\int {\rm{d}}{xp}(x| D){[{\partial }_{D}\mathrm{ln}p(x| D)]}^{2}$. Here $p(x| D)$ is the conditional probability that the outcome of the measurement is x, if the parameter value is D. If we define the quantum Fisher information (QFI) H(D) as the supremum of the FI over all possible quantum measurements described by positive operator-valued measures, we obtain the quantum Cramér–Rao bound, i.e. the ultimate lower bound on the achievable precision in the estimation of D [28].

The QFI with respect to D can be evaluated analytically for Gaussian states [34, 39, 40] and for the case at hand is given by

$$\begin{eqnarray}&&H(D)=\mathrm{Tr}[{{\rm{\Omega }}}^{T}\partial \sigma {\rm{\Omega }}{\rm{\Phi }}],\end{eqnarray} \tag{ 5 }$$

where ∂ denotes the derivative with respect to D, ${\rm{\Omega }}={\oplus }_{i=1}^{2}\left(\genfrac{}{}{0em}{}{\ 0\ \ 1}{-1\ 0}\right)$ is the symplectic matrix, and the matrix Φ satisfies the linear equation $\partial \sigma =2\sigma {\rm{\Omega }}{\rm{\Phi }}{{\rm{\Omega }}}^{T}\sigma -\tfrac{1}{2}{\rm{\Phi }}.$ The expression for the QFI is quite cumbersome and we report it in appendix B. In the following we analyze its behavior in the cases of interest for our discussion. We also note that the QFI is a homogeneous function of the parameters of the Hamiltonian: $H(\alpha D,\alpha {\omega }_{a},\alpha {\omega }_{b},\alpha \lambda )$ $=\,{\alpha }^{-2}H(D,{\omega }_{a},{\omega }_{b},\lambda )$.

The QFI for the two-mode ground state is a relevant benchmark for the precision of any measurement on the system. In realistic experiments, however, it will only be possible to measure a limited number of system observables. For definiteness, we shall focus here on typical measurements that can be performed on the transverse radiation field a. It is thus interesting to compute also the QFI H_a(D) for the reduced state of mode a, which is a squeezed thermal state, see appendix B.

Homodyne detection corresponds to the measurement of the quadrature operator $\hat{x}(\phi )$ $=\,(a{{\rm{e}}}^{-{\rm{i}}\phi }+{a}^{\dagger }{{\rm{e}}}^{{\rm{i}}\phi })/$ $\sqrt{2}$, where ϕ is an arbitrary angle. For the reduced state ${\rho }_{a}$, the probability distribution for the outcome of $x(\phi )$ is a normal distribution with zero mean and covariance ${\sigma }_{11}{\cos }^{2}\phi +{\sigma }_{22}{\sin }^{2}\phi ;$ the corresponding FI, ${F}^{(\mathrm{hd})}$ is reported in appendix B. It is further argued that the FI has a maximum at $\phi =0$: this has the important consequence that the best measurement angle does not depend on the (yet unknown) D parameter. On the other hand, the FI for photon counting cannot be determined analytically as it involves an infinite sum of the terms $p(n)=\langle n| \rho | n\rangle$, i.e. the probabilities of finding a photon in the Fock state $| n\rangle$ [34, 41]; we evaluate it numerically in a truncated Fock space.

The main results of this paper, namely the analytical expressions for the QFI and FI for the homodyne detection (reported in appendix B) are valid for any values of the model parameters, provided the Hamiltonian (2) admits a ground state. In what follows, we focus on the parameter regimes that are most relevant from a theoretical and from an experimental point of view.

Let us start by focusing on the regime of relatively small coupling within the USC. The QFI has a non vanishing limiting value

$$\begin{eqnarray}&&H(D)=2{(4D+{\omega }_{a})}^{-2}+O({\lambda }^{2}),\end{eqnarray} \tag{ 6 }$$

showing how, in this regime, smaller values of D can be estimated more efficiently, yet a large number of measurements is needed to neutralize the contribution of ${\omega }_{a}$ and obtain a precise estimation. Remarkably, the FI for the homodyne measurement saturates H(D) up to second order in λ, indicating that, no matter what the true value of D is, we can estimate it with optimal efficiency by detecting the quadrature of mode a when λ is sufficiently small. In detail, if we consider the plausible scenario $D\lesssim {D}_{\mathrm{TRK}}$, the ratio between FI and H(D) is

$$\begin{eqnarray}&&\displaystyle \frac{{F}^{(\mathrm{hd})}(D\lesssim {D}_{\mathrm{TRK}})}{H(D\lesssim {D}_{\mathrm{TRK}})}=1-\displaystyle \frac{8{\omega }_{a}^{2}}{{({\omega }_{a}+{\omega }_{b})}^{4}}{\lambda }^{2}+O({\lambda }^{4}).\end{eqnarray} \tag{ 7 }$$

We can now investigate the estimation properties of the diamagnetic parameter D without constraining the value of the coupling constant λ and thus, when possible, also at the critical point. For the case D = 0 the phase transition appears when the coupling reaches the critical value ${\lambda }_{\mathrm{crit}}=\sqrt{{\omega }_{a}{\omega }_{b}}/2$. In this case the QFI for D diverges as

$$\begin{eqnarray}&&H(D=0)\sim {\omega }_{b}/[8{\omega }_{a}{(\lambda -{\lambda }_{\mathrm{crit}})}^{2}],\end{eqnarray} \tag{ 8 }$$

and the homodyne and the photon counting measurements both saturate the QFI in this limit, as shown in figure 2 (left panel). If D approaches the TRK value $D\simeq {D}_{\mathrm{TRK}}$ (see equation (4)), the phase transition is suppressed and there is no criticality in the system. The QFI for the two-mode state (for any λ) reads

$$\begin{eqnarray}&&H({D}_{\mathrm{TRK}})=\displaystyle \frac{2}{{\omega }_{a}^{2}}-\displaystyle \frac{16{\lambda }^{2}{\omega }_{a}{\omega }_{b}}{{[4{\lambda }^{2}{\omega }_{a}+{\omega }_{b}({\omega }_{a}+{\omega }_{b}){}^{2}]}^{2}}.\end{eqnarray} \tag{ 9 }$$

$H({D}_{\mathrm{TRK}})$ is clearly higher for small frequencies of the radiation mode, whereas it decreases in the case of high frequencies. With respect to λ, $H({D}_{\mathrm{TRK}})$ has a minimum at ${\lambda }^{2}={\omega }_{b}{({\omega }_{a}+{\omega }_{b})}^{2}/4{\omega }_{a}$ and the asymptotic value $2/{\omega }_{a}^{2}$ for large and small λ. The ratio between the FI for the two measurements and the QFI is shown in figure 3 as a function of λ and ${\omega }_{a}$. One can also verify that the expansion in equation (7) is in fact a good guideline for remarkably large values of λ, implying that homodyne would be an excellent measurement strategy in most experimentally relevant situations. For example, setting $D={D}_{\mathrm{TRK}}$ and ${\omega }_{a}={\omega }_{b}=2\lambda$, we still get a ratio of $81/92\simeq 0.88$, while the approximation of equation (7) would yield $7/8\simeq 0.87$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In general, in the vicinity of the critical point ${D}_{\mathrm{crit}}$, the asymptotic behavior of the two-mode QFI is

$$\begin{eqnarray}&&H(D)\mathop{\sim }\limits_{D\to {D}_{\mathrm{crit}}}\displaystyle \frac{1}{8{(D-{D}_{\mathrm{crit}})}^{2}},\end{eqnarray} \tag{ 10 }$$

that is, the QFI diverges, meaning that the amount of information that can be extracted per measurement increases enormously. This is another confirmation of the role of quantum criticality in the enhancement of estimation [34, 42, 43]. We point out, however, that we neglect finite-size effects which may become relevant near the critical point. Nonetheless, it is reasonable to believe that these effects translate into a smoothing of the QFI, for either the atomic ensemble or the radiation mode (see [33] for the finite-size case without the diamagnetic term).

The QFI for the reduced state of one mode saturates the two-mode QFI H(D) when the system is close to the critical point. This is remarkable: optimal estimation of D around criticality can be achieved by measuring only a part of the system. Moreover, it can effectively be achieved with feasible experiments such as homodyne detection or photon counting (see figure 2). The FI for the x-quadrature measurement indeed saturates the QFI as D gets close to the critical value ${D}_{\mathrm{crit}}$,

$$\begin{eqnarray}&&\displaystyle \frac{{F}^{(\mathrm{hd})}(D)}{H(D)}\sim 1-\displaystyle \frac{16{\lambda }^{2}{\omega }_{a}^{3/2}}{4{\lambda }^{2}{\omega }_{a}{\omega }_{b}+{\omega }_{b}^{4}}{(D-{D}_{\mathrm{crit}})}^{1/2}.\end{eqnarray} \tag{ 11 }$$

We further find numerically that the photon counting measurements can also saturate the QFI near the critical point.

We finally tackle the problem of discriminating between the two values of D that are most argued in the recent circuit-QED literature: D = 0 and ${D}_{\mathrm{TRK}}$. By making a measurement on a subsystem, typically the cavity field, we want to decide which of the two hypotheses is correct. Quantum mechanics poses a limit to the discriminating ability, quantified by the Helstrom bound [36, 37]: if we label ${\rho }^{(0)}$ and ${\rho }^{(\mathrm{TRK})}$, respectively, the ground state in the two hypotheses, the lowest probability of error in the discrimination between these two states is ${P}_{e}=\tfrac{1}{2}[1-{\mathbb{D}}({\rho }^{(0)},{\rho }^{(\mathrm{TRK})})]$, where ${\mathbb{D}}$ denotes the trace distance. The Helstrom bound is based on the assumption that the system has the same probability to be in either of the two states, reflecting the lack of any a priori information on which hypothesis is the correct one. The Helstrom bound calculated on the two-mode state can only be used as an ideal benchmark, as it is saturated by a (typically unfeasible) collective measurement on the two modes. A more practical benchmark is given by the Helstrom bound calculated on the reduced state of mode a, i.e. in the discrimination between ${\rho }_{a}^{(0)}$ and ${\rho }_{a}^{(\mathrm{TRK})}$. Being these states mixed, the trace distance cannot be computed analytically, and a diagonalization of the density operators expanded on a truncated Fock basis is required. The ensuing bound on the probability of error, ${P}_{e}^{(a)}$, is less stringent than the two-mode bound.

We show the dependence of P_e and ${P}_{e}^{(a)}$ on $\lambda /{\lambda }_{c}$ in figure 4: the probability of error is close to 1/2 for small λ and vanishes when $\lambda \to {\lambda }_{c}=\sqrt{{\omega }_{a}{\omega }_{b}}/2$ in the case D = 0. The dependence on λ in the neighborhood of ${\lambda }_{c}$ is ${P}_{e}\sim {({\lambda }_{c}-\lambda )}^{1/4}$. For ${P}_{e}^{(a)}$ a numerical fit shows an exponent around 1/5. Thus we find that the discrimination between the two hypotheses is hard in the regime of small coupling, while criticality allows for a great improvement, should the optimal strategy be found.

Zoom In Zoom Out Reset image size

To this end, we checked the performance of two feasible discrimination strategies using either the $\hat{x}(0)$ quadrature measurement or photon counting on the mode a. The former exploits the fact that the probability distribution of the outcome of the quadrature measurement is a Gaussian with variance ${\sigma }_{11}$. In the case D = 0, when $\lambda \to {\lambda }_{\mathrm{crit}}$, ${\sigma }_{11}\sim {({\lambda }_{\mathrm{crit}}-\lambda )}^{-1/2}$. On the other hand, for $D={\lambda }^{2}/{\omega }_{b}$, ${\sigma }_{11}$ does not depart much from its limit for $\lambda \to 0$, which is 1/2. Thus, close to the critical value, the Gaussian distribution for $D={\lambda }^{2}/{\omega }_{b}$ is very narrow compared to the distribution for D = 0.

We can thus set up a readily feasible discrimination experiment as follows: if the outcome of the experiment is⁶ $| x| \lt 2{\sigma }_{11}^{(\mathrm{TRK})}$ then the state of the system is ${\rho }^{(\mathrm{TRK})}$, otherwise it is ${\rho }^{(0)}$. The corresponding probability of error is

$$\begin{eqnarray}&&{P}_{e}^{(\mathrm{hd})}={\int }_{0}^{2{\sigma }_{11}^{(\mathrm{TRK})}}{ \mathcal N }(0,{\sigma }_{11}^{(0)}){\rm{d}}x+{\int }_{2{\sigma }_{11}^{(\mathrm{TRK})}}^{\infty }{ \mathcal N }(0,{\sigma }_{11}^{(\mathrm{TRK})}){\rm{d}}x,\end{eqnarray} \tag{ 12 }$$

where the first (second) term is the probability of detecting ${\rho }^{(0)}$ (${\rho }^{(\mathrm{TRK})}$) when the actual state was ${\rho }^{(\mathrm{TRK})}$ (${\rho }^{(0)}$), and ${ \mathcal N }(\mu ,\sigma )$ is the normal distribution with mean μ and variance σ.

In a photon counting experiment, we can exploit the fact that, when D = 0, the average photon number diverges near the phase transition [34], while it is close to zero if the phase transition is suppressed by the presence of the A² term. Thus, we can discriminate between the two hypotheses by setting a threshold photon number n_T and assigning any outcome below that threshold to the hypothesis ${\rho }^{(\mathrm{TRK})}$ and any outcome above to ${\rho }^{(0)}$. The probability of error in this case would be

$$\begin{eqnarray}&&{P}_{e}^{(\mathrm{pc})}=\displaystyle \frac{1}{2}\left(\displaystyle \sum _{n=0}^{{n}_{{\rm{T}}}}{p}_{0}(n)+1-\displaystyle \sum _{n=0}^{{n}_{{\rm{T}}}}{p}_{\mathrm{TRK}}(n)\right).\end{eqnarray} \tag{ 13 }$$

We find that the optimal threshold value is n_T = 0, i.e. the discrimination is between no photons and any number of photons. This is a harder question to settle in a realistic experiment.

The resulting ${P}_{e}^{(\mathrm{hd})}$ and ${P}_{e}^{(\mathrm{pc})}$ are compared to the Helstrom bounds in figure 4. We see that the homodyne scheme is slightly better than the photon counting one. Interestingly, these two feasible discrimination schemes have the same behavior as the Helstrom bound ${P}_{e}^{(a)}$ on the reduced state when approaching the critical point, although neither of them is optimal.

Before closing, we would like to comment on the feasibility of the proposed measurements, which require access to the cavity mode a (see figure 1). One possible way to access the intra-cavity field is to suddenly switch off the coupling λ, a drastic yet experimentally feasible procedure [1]. Allowing one of the cavity mirrors to have a small but finite transmissivity, one may subsequently collect the cavity output field (i.e. the radiation that gradually leaks out into the external world). In absence of light–matter coupling and other losses, the cavity output field can be used to extract the full quantum statistics of mode a just before the switch-off [44]. A crucial open question, however, is how coherent the switching process can be, i.e., how well it can preserve the quantum state of light. Importantly, our results can be generalized in future work to take into account experimental imperfections in the switching process: due to finite transients, losses and other decoherence mechanisms, one would end up measuring a deteriorated version of the original cavity field a. It is reasonable to assume that the whole process could be modeled as a quantum channel acting on the reduced state of the cavity mode, just before the latter is measured. Insofar as these imperfections can be described by a Gaussian channel, the problem could be attacked by a straightforward generalization of the tools employed here. Importantly, the relevant noise parameters must be known in advance for this approach to work within the paradigm of single-parameter estimation. As an example, we illustrate the special case of a Gaussian pure-loss channel [45], described by a loss parameter $\eta \in [0,1]$ quantifying the probability of single-photon loss (with $\eta =1$ meaning complete loss, i.e. the output state is always the vacuum). The effect on the estimability of D is shown in figure 5. The QFI and FI inevitably decrease with η, but the ratio does not, signifying that homodyne detection may be a robust measurement to estimate the diamagnetic term in presence of decoherence.

Zoom In Zoom Out Reset image size

We investigated the detection of the diamagnetic term in a Dicke model of light–matter interaction, formulating the problem in terms of quantum parameter estimation or quantum state discrimination. We obtained the ultimate quantum limits to the rate at which information can be extracted from the ground state of the coupled system, and we discussed the performance of two typical measurements on the cavity field, homodyne and photon counting, that allow for efficient estimation of the parameter of interest. This efficiency becomes optimal in experimentally relevant regimes, such as the small-coupling regime in which the coupling is up to 20% of the mode frequencies.

To conclude, we would like to indicate some possible directions for future work. By employing multiparameter quantum estimation theory [28], our study can be expanded to cover more realistic situations in which other parameters of the model, e.g. loss rates or the coupling λ, are not known exactly. The role of finite-size effects could also be explored, in particular near criticality where a coupled-oscillator model might be inaccurate [33] (this, however, will require techniques beyond the Gaussian formalism employed here). The estimability of the A²-term from a thermal state of the coupled system could also be addressed: in this scenario radiation will be continuously emitted by the cavity, without the need to have a fast modulation of the coupling constant.

Finally, we remark that the ideas presented in this letter may be generalized to test the validity of alternative or more complex models of light–matter interaction; for example, alternative Dicke models for circuit QED such as that proposed in [19], or models with additional terms to describe electrostatic contributions (e.g. dipole–dipole interactions between the atoms), or effective contributions from higher harmonics of the cavity field [32].

This work has been supported by EU through the Erasmus+ Placement programme, the Marie Skłodowska-Curie Action H2020 (project ConAQuMe) and the collaborative Project QuProCS (Grant Agreement 641277), the EPSRC UK (grant No. EP/K026267/1), the European Research Council (ERC StG GQCOP, Grant No. 637352), the Royal Society (Grant No. IE150570), and by UniMI through the H2020 Transition Grant 15-6-3008000-625.

As is shown in [32], the USC Hamiltonian in equation (2) can be cast to the diagonal form of equation (2), where the frequencies ${\omega }_{U}$, ${\omega }_{L}$ are

$$\begin{eqnarray}&&{\omega }_{U,L}^{2}=\displaystyle \frac{{\omega }_{a}^{2}+4D{\omega }_{a}+{\omega }_{b}^{2}}{2}\pm \sqrt{{\left(\displaystyle \frac{{\omega }_{a}^{2}+4D{\omega }_{a}-{\omega }_{b}^{2}}{2}\right)}^{2}+4{\lambda }^{2}{\omega }_{a}{\omega }_{b}}.\end{eqnarray} \tag{ A.1 }$$

The symplectic matrix that transforms the vector $(a,b,{a}^{\dagger },{b}^{\dagger })$ of the original-mode creation and annihilation operators into the vector of polaritonic modes $({p}_{U},{p}_{L},{p}_{U}^{\dagger },{p}_{L}^{\dagger })$, is given by

$$\begin{eqnarray}\bar{S}=\left(\begin{array}{cccc}\cos \theta \ {f}_{+}\left(\tfrac{{\omega }_{U}}{{\omega }_{a}}\right) & -\sin \theta \ {f}_{+}\left(\tfrac{{\omega }_{U}}{{\omega }_{b}}\right) & \cos \theta \ {f}_{-}\left(\tfrac{{\omega }_{U}}{{\omega }_{a}}\right) & -\sin \theta \ {f}_{-}\left(\tfrac{{\omega }_{U}}{{\omega }_{b}}\right)\\ \sin \theta \ {f}_{+}\left(\tfrac{{\omega }_{L}}{{\omega }_{a}}\right) & +\cos \theta \ {f}_{+}\left(\tfrac{{\omega }_{L}}{{\omega }_{b}}\right) & \sin \theta \ {f}_{-}\left(\tfrac{{\omega }_{L}}{{\omega }_{a}}\right) & \cos \theta \ {f}_{-}\left(\tfrac{{\omega }_{L}}{{\omega }_{b}}\right)\\ \cos \theta \ {f}_{-}\left(\tfrac{{\omega }_{U}}{{\omega }_{a}}\right) & -\sin \theta \ {f}_{-}\left(\tfrac{{\omega }_{U}}{{\omega }_{b}}\right) & \cos \theta \ {f}_{+}\left(\tfrac{{\omega }_{U}}{{\omega }_{a}}\right) & -\sin \theta \ {f}_{+}\left(\tfrac{{\omega }_{U}}{{\omega }_{b}}\right)\\ \sin \theta \ {f}_{-}\left(\tfrac{{\omega }_{L}}{{\omega }_{a}}\right) & \cos \theta \ {f}_{-}\left(\tfrac{{\omega }_{L}}{{\omega }_{b}}\right) & \sin \theta \ {f}_{+}\left(\tfrac{{\omega }_{L}}{{\omega }_{a}}\right) & +\cos \theta \ {f}_{+}\left(\tfrac{{\omega }_{L}}{{\omega }_{b}}\right)\end{array}\right),\end{eqnarray} \tag{ A.2 }$$

where

$$\begin{eqnarray}&&\cos 2\theta =\displaystyle \frac{{\omega }_{a}^{2}+4D{\omega }_{a}-{\omega }_{b}^{2}}{{\omega }_{U}^{2}-{\omega }_{L}^{2}}\qquad \sin (2\theta )=-\displaystyle \frac{4\lambda \sqrt{{\omega }_{a}{\omega }_{b}}}{{\omega }_{U}^{2}-{\omega }_{L}^{2}}.\end{eqnarray} \tag{ A.3 }$$

And we define

$$\begin{eqnarray}&&{f}_{\pm }(x)=\displaystyle \frac{1}{2}\left(\sqrt{x}\pm \displaystyle \frac{1}{\sqrt{x}}\right).\end{eqnarray} \tag{ A.4 }$$

The ground state of the diagonalized Hamiltonian is the vacuum, i.e. a Gaussian state with no displacement and a covariance matrix ${\sigma }_{0}={\mathbb{I}}/2$ in the basis of the quadratures ${x}_{U,L}=({p}_{U,L}$ $+\,{p}_{U,L}^{\dagger })/\sqrt{2}$, ${y}_{U,L}=-{\rm{i}}$ $({p}_{U,L}-{p}_{U,L}^{\dagger })/\sqrt{2}$.

To obtain the covariance matrix σ in the original modes a and b, transform ${\sigma }_{0}$ using a symplectic matrix S that can easily be obtained from $\bar{S}$. After some simplifications, we have

$$\begin{eqnarray}\sigma =S{\sigma }_{0}{S}^{T}=\left(\begin{array}{cccc}\tfrac{{\omega }_{a}({\omega }_{b}^{2}+{\omega }_{L}{\omega }_{U})}{{\omega }_{L}{\omega }_{U}({\omega }_{L}+{\omega }_{U})} & 0 & -\tfrac{2\lambda {\omega }_{a}{\omega }_{b}}{{\omega }_{L}{\omega }_{U}({\omega }_{L}+{\omega }_{U})} & 0\\ 0 & \tfrac{2({\omega }_{a}^{2}+4D{\omega }_{a}+{\omega }_{L}{\omega }_{U})}{2{\omega }_{a}{\omega }_{L}+2{\omega }_{a}{\omega }_{U}} & 0 & \tfrac{2\lambda }{{\omega }_{L}+{\omega }_{U}}\\ -\tfrac{2\lambda {\omega }_{a}{\omega }_{b}}{{\omega }_{L}{\omega }_{U}({\omega }_{L}+{\omega }_{U})} & 0 & \tfrac{{\omega }_{b}({\omega }_{a}^{2}+4D{\omega }_{a}+{\omega }_{L}{\omega }_{U})}{{\omega }_{L}{\omega }_{U}({\omega }_{L}+{\omega }_{U})} & 0\\ 0 & \tfrac{2\lambda }{{\omega }_{L}+{\omega }_{U}} & 0 & \tfrac{{\omega }_{b}^{2}+{\omega }_{L}{\omega }_{U}}{{\omega }_{b}({\omega }_{L}+{\omega }_{U})}\end{array}\right).\end{eqnarray} \tag{ A.5 }$$

B.1. QFI for the two-mode state

Given σ, we can calculate the QFI using equation (5). If the state is pure, the solution of equation (6) is ${\rm{\Phi }}=-\partial \sigma$. This in turn yields the following expression for the QFI:

$$\begin{eqnarray}&&H(D)=2(\partial {\sigma }_{\mathrm{1,1}}\partial {\sigma }_{\mathrm{2,2}}+2\partial {\sigma }_{\mathrm{2,4}}\partial {\sigma }_{\mathrm{3,1}}+\partial {\sigma }_{\mathrm{3,3}}\partial {\sigma }_{\mathrm{4,4}}).\end{eqnarray} \tag{ B.1 }$$

Here we report the general expression of the QFI (after some manipulations)

$$\begin{eqnarray}H & = & \displaystyle \frac{2{\omega }_{a}^{2}{\omega }_{b}^{2}}{{\omega }_{L}^{4}{\omega }_{U}^{4}{({\omega }_{L}+{\omega }_{U})}^{4}}\{({\omega }_{b}^{2}[16{D}^{2}{\omega }_{a}^{2}+8D{\omega }_{a}({\omega }_{a}^{2}+3{\omega }_{b}^{2}+2{\omega }_{L}{\omega }_{U})\\ & & +{\omega }_{a}^{4}+4{\omega }_{L}{\omega }_{U}({\omega }_{a}^{2}+{\omega }_{b}^{2})+6{\omega }_{a}^{2}{\omega }_{b}^{2}+{\omega }_{b}^{4}]\\ & & -8{\lambda }^{2}{\omega }_{a}{\omega }_{b}(4D{\omega }_{a}+{\omega }_{a}^{2}+2{\omega }_{b}^{2})+32{\lambda }^{4}{\omega }_{a}^{2}\}.\end{eqnarray} \tag{ B.2 }$$

The relevant cases are discussed in the paper.

B.2. QFI for the reduced state of mode a

The covariance matrix ${\sigma }_{a}$ of the reduced state ${\rho }_{a}$ of the photonic mode is simply the upper left diagonal block of σ. The corresponding state is a squeezed thermal state

$$\begin{eqnarray}&&{\rho }_{a}=S(r){\nu }_{\mathrm{th}}{S}^{\dagger }(r)=\displaystyle \frac{1}{1+N}\displaystyle \sum _{m}^{\infty }{\left(\displaystyle \frac{N}{1+N}\right)}^{m}S(r)| m\rangle \langle m| {S}^{\dagger }(r),\end{eqnarray} \tag{ B.3 }$$

where N is the number of thermal photons and r is the squeezing parameter:

$$\begin{eqnarray}&&N=\displaystyle \frac{1}{2}\left(\sqrt{\displaystyle \frac{({\omega }_{b}^{2}+{\omega }_{L}{\omega }_{U})(4D{\omega }_{a}+{\omega }_{a}^{2}+{\omega }_{L}{\omega }_{U})}{{\omega }_{L}{\omega }_{U}({\omega }_{L}+{\omega }_{U}){}^{2}}}-1\right)\end{eqnarray} \tag{ B.4 }$$

$$\begin{eqnarray}&&r=\displaystyle \frac{1}{4}\mathrm{log}\displaystyle \frac{{\omega }_{L}{\omega }_{U}(4D{\omega }_{a}+{\omega }_{a}^{2}+{\omega }_{L}{\omega }_{U})}{{\omega }_{a}^{2}({\omega }_{b}^{2}+{\omega }_{L}{\omega }_{U})}.\end{eqnarray} \tag{ B.5 }$$

The QFI for the reduce state can be obtained easily by solving equation (6) of the paper for Φ. Being ${\sigma }_{a}$ diagonal, we simply find that Φ is diagonal with

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{11}=2\displaystyle \frac{{(\partial {\sigma }_{a})}_{11}+4{({\sigma }_{a})}_{11}^{2}{(\partial {\sigma }_{a})}_{22}}{16{({\sigma }_{a})}_{11}^{2}{({\sigma }_{a})}_{22}^{2}-1},\end{eqnarray} \tag{ B.6 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{22}=2\displaystyle \frac{{(\partial {\sigma }_{a})}_{22}+4{({\sigma }_{a})}_{22}^{2}{(\partial {\sigma }_{a})}_{11}}{16{({\sigma }_{a})}_{11}^{2}{({\sigma }_{a})}_{22}^{2}-1}.\end{eqnarray} \tag{ B.7 }$$

The corresponding QFI is

$$\begin{eqnarray}&&{H}_{a}=\displaystyle \frac{4(2{\sigma }_{11}^{2}{(\partial {\sigma }_{22})}^{2}+2{\sigma }_{22}^{2}{(\partial {\sigma }_{11})}^{2}+\partial {\sigma }_{11}\partial {\sigma }_{22})}{16{\sigma }_{11}^{2}{\sigma }_{22}^{2}-1}.\end{eqnarray} \tag{ B.8 }$$

B.3. FI for the homodyne detection

In a homodyne detection experiment one can measure the field mode quadrature at an arbitrary phase, i.e. the expected value of the operator $\hat{x}(\phi )=(a{{\rm{e}}}^{-{\rm{i}}\phi }+{a}^{\dagger }{{\rm{e}}}^{{\rm{i}}\phi })/\sqrt{2}$.

The probability density for the outcome of a homodyne measurement is easily obtained from the Wigner function of the reduced state of the photon mode ${\rho }_{a}={\mathrm{Tr}}_{b}[\rho ]$, which is a Gaussian distribution with zero mean and covariance matrix $\mathrm{diag}({\sigma }_{11},{\sigma }_{22})$.

The probability ${p}_{\phi }(x)$ of the outcoume $x(\phi )$ for the measurement is the marginal distribution of the Wigner function, obtained after a rotation of an angle ϕ in the phase space:

$$\begin{eqnarray}{p}_{\phi }(x) & = & \displaystyle \int {\rm{d}}{yW}[{\rho }_{a}](x\cos \phi -y\sin \phi ,x\sin \phi +y\cos \phi )\\ & = & { \mathcal N }(0,{\sigma }_{11}{\cos }^{2}\phi +{\sigma }_{22}{\sin }^{2}\phi ),\end{eqnarray} \tag{ B.9 }$$

that is, a normal distribution with variance ${\sigma }_{11}{\cos }^{2}\phi +{\sigma }_{22}{\sin }^{2}\phi$. The FI for this distribution is easily obtained to be

$$\begin{eqnarray}&&{F}_{\phi }(D)=\int {\rm{d}}x\displaystyle \frac{{[\partial {p}_{\phi }(x)]}^{2}}{{p}_{\phi }(x)}=\displaystyle \frac{{(\partial {\sigma }_{11}{\cos }^{2}\phi +\partial {\sigma }_{22}{\sin }^{2}\phi )}^{2}}{2{({\sigma }_{11}{\cos }^{2}\phi +{\sigma }_{22}{\sin }^{2}\phi )}^{2}}.\end{eqnarray} \tag{ B.10 }$$

It is easy to check that the FI ${F}_{\phi }(D)$ has maxima for $\phi =0,\pi /2$, in which

$$\begin{eqnarray}&&{F}_{0}(D)=\displaystyle \frac{{(\partial {\sigma }_{11})}^{2}}{2{\sigma }_{11}},\qquad {F}_{\tfrac{\pi }{2}}(D)=\displaystyle \frac{{(\partial {\sigma }_{22})}^{2}}{2{\sigma }_{22}}.\end{eqnarray} \tag{ B.11 }$$

Numerical analysis indicates that ${F}_{0}(D)\geqslant {F}_{\tfrac{\pi }{2}}(D)$, thus $\phi =0$ is the optimal angle to perform the homodyne measurement.

Our calculations can be readily adapted to the study of Dicke models in the electric dipole gauge. In such a case, our Hamiltonian may be written as

$$\begin{eqnarray}&&{{ \mathcal H }}_{\mathrm{dip}}={\bar{\omega }}_{a}{a}^{\dagger }a+{\bar{\omega }}_{b}{b}^{\dagger }b+\bar{\lambda }(a+{a}^{\dagger })(b+{b}^{\dagger })+\bar{D}{(b+{b}^{\dagger })}^{2},\end{eqnarray} \tag{ C.1 }$$

where ${\bar{\omega }}_{a},{\bar{\omega }}_{b}$, $\bar{\lambda }$, $\bar{D}$ are the bare frequencies and light–matter coupling constant relevant to the new gauge. Note how the coupling parameter $\bar{D}$ is now associated with the P² term [10]. In the dipole gauge, the operators $b,{b}^{\dagger }$ describe physical degrees of freedom of matter, thanks to the equivalence between canonical and kinetic momentum. On the other hand, the field operators $a,{a}^{\dagger }$ no longer describe the transverse radiation field, but are 'contaminated' by matter properties [7]. From equation (C.1) it is evident that the results presented in our manuscript can be easily translated in the dipole gauge, by simply swapping the role of a and b in all calculations. Then, the parameter to be probed becomes the P² coupling constant $\bar{D}$. Specifically, our calculations relative to homodyne detection and photon counting indicate that efficient estimation of the P² term is achievable through measurements on the matter degree of freedom.