Rashba cavity QED: a route towards the superradiant quantum phase transition

TL;DR

This paper presents a theoretical framework for Rashba cavity QED, showing that a superradiant quantum phase transition can occur in a 2D electron gas with spin-orbit coupling under realistic conditions.

Contribution

It introduces a novel theory of cavity QED with Rashba spin-orbit coupling, predicting a superradiant phase transition driven by soft spin-flip transitions.

Findings

Lowest polaritonic frequency can vanish, indicating a superradiant phase transition.

The transition is linked to magnetostatic instability from non-zero wave vector spin-flip modes.

The theory applies to realistic experimental parameters.

Abstract

We develop a theory of cavity quantum electrodynamics for a 2D electron gas in the presence of Rashba spin-orbit coupling and perpendicular static magnetic field, coupled to spatially nonuniform multimode quantum cavity photon fields. We demonstrate that the lowest polaritonic frequency of the full Hamiltonian can vanish for realistic parameters, achieving the Dicke superradiant quantum phase transition. This singular behaviour originates from soft spin-flip transitions possessing a non-vanishing dipole moment at non-zero wave vectors and can be viewed as a magnetostatic instability.

Equations

$${\mathcal{H}}_{\text{2DEG}}=$$

$${ϵ}_{\eta ,s}={\omega }_c\left(\eta +\frac{s}{2}\sqrt{1+8\eta {\gamma }^2}\right),$$

$${\alpha }^2=\frac{{\eta }_1+{\eta }_2-\sqrt{4{\eta }_1{\eta }_2+1}}{2(m^{\ast }{l}_B{)}^2}.$$

$${\mathbf{A}}^{\mathrm{cav}}(\mathbf{r})=$$

$$\times \left(a_{q_x,n_z}{e}^{i{q}_{x}x}+a_{q_x,n_z}^{†}{e}^{-i{q}_{x}x}\right),$$

$$\frac{c}{2\pi }\sqrt{c^2{q}_x^2-\epsilon {\omega }^2}\coth \frac{L_z\sqrt{q_x^2-\epsilon {\omega }^2/c^2}}{2}=-Q_{yy}(q_x,\omega ),$$

$$Q_{yy}(q_x,\omega )=\frac{n_{\mathrm{qw}}{n}_e{e}^2}{m^{\ast }}\left[1-\frac{{\omega }_c}{\nu }\sum _{\ell \le \nu <{\ell }^{\prime }}\frac{{\omega }_{{\ell }^{\prime }\ell }^{\mathrm{LL}}[g_{q_x}^{{\ell }^{\prime }\ell }{]}^2}{({\omega }_{{\ell }^{\prime }\ell }^{\mathrm{LL}}{)}^2-{\omega }^2}\right].$$

$$g_{q_x}^{{\ell }^{\prime }\ell }=$$

$$+\sin {\theta }_{{\ell }^{\prime }}\sin {\theta }_{\ell }\left(\sqrt{\eta -1}{\Theta }_{{\eta }^{\prime }-1}^{\eta -2}-\sqrt{\eta }{\Theta }_{{\eta }^{\prime }-1}^{\eta }\right)$$

$$+\cos {\theta }_{{\ell }^{\prime }}\cos {\theta }_{\ell }\left(\sqrt{\eta }{\Theta }_{{\eta }^{\prime }}^{\eta -1}-\sqrt{\eta +1}{\Theta }_{{\eta }^{\prime }}^{\eta +1}\right),$$

$${\Theta }_{n_2}^{n_1}=\sqrt{\frac{m!}{M!}}{e}^{-\frac{\xi }{2}}{\xi }^{\frac{M-m}{2}}{L}_m^{M-m}(\xi )S^{n_2-n_1},$$

$$\Delta \alpha \sim \frac{1}{{\nu }^2}\frac{n_{\mathrm{qw}}{n}_e{e}^2}{(m^{\ast }c{)}^2},$$

$$F[\delta \mathbf{B}]=$$

$$+\int \left[-{\mathbf{M}}^{\mathrm{eq}}\cdot \delta \mathbf{B}-\frac{1}{2}{\chi }_{kl}\delta B_k\delta B_l+\dots \right]d^3\mathbf{r},$$

$${\chi }_{zz}(\mathbf{r},{\mathbf{r}}^{\prime })=$$

$$\times \delta (y-y^{\prime })\delta (z-L_z/2)\delta (z^{\prime }-L_z/2).$$

$${\mathbf{B}}^{\mathrm{ext}}+\delta \mathbf{B}(\mathbf{r})=4\pi \mathbf{M}(\mathbf{r})+\nabla \int \frac{d^3{\mathbf{r}}^{\prime }}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}{\nabla }^{\prime }\cdot \mathbf{M}({\mathbf{r}}^{\prime }),$$

$$M_k(\mathbf{r})=M_k^{\mathrm{ext}}(\mathbf{r})+M_k^{\mathrm{eq}}(\mathbf{r})+\int {\chi }_{kl}(\mathbf{r},{\mathbf{r}}^{\prime })\delta B_l({\mathbf{r}}^{\prime })d^3{\mathbf{r}}^{\prime }.$$

$$\delta \mathbf{M}(\mathbf{r})=M_{q_x}{\mathbf{u}}_z\delta (z-L_z/2)e^{i{q}_{x}x},$$

$$\delta B_x(\mathbf{r})=-2i\pi q_x{M}_{q_x}{e}^{-q_x|z-L_z/2|-i{q}_{x}x}\mathrm{sign}(z-L_z/2),$$

$$\delta B_z(\mathbf{r})=2\pi q_x{M}_{q_x}{e}^{-q_x|z-L_z/2|-i{q}_{x}x},$$

$$-2\pi q_x\frac{Q_{yy}(q_x,\omega =0)}{c^2{q}_x^2}=1.$$

$${\psi }_{\eta ,s,p_x}(x,y)=\frac{e^{i{p}_{x}x}}{\sqrt{L_x}}\left(\begin{array}{c}\sin {\theta }_{\eta s}{\varphi }_{\eta -1}(y-p_x{l}_B^2)\\ i\cos {\theta }_{\eta s}{\varphi }_{\eta }(y-p_x{l}_B^2)\end{array}\right),$$

$${\varphi }_{\eta }(y)=\frac{e^{-y^2/(2l_B^2)}}{\sqrt{2^{\eta }\eta !l_B\sqrt{\pi }}}{H}_{\eta }(y/l_B),$$

$$|F⟩=\prod _{j=1}^{n_{\mathrm{qw}}}\prod _{\ell \le \nu }\prod _{p_x}{c}_{j,\ell ,p_x}^{†}|0⟩,$$

$$b_{q_x,{\ell }^{\prime }\ell }^{†}=\sqrt{\frac{2\pi l_B^2}{n_{\mathrm{qw}}{L}_x{L}_y}}\sum _{j,p_x}{c}_{j,{\ell }^{\prime },p_x+q_x}^{†}{c}_{j,\ell ,p_x},$$

$$H=H_{\mathrm{cav}}+H_{\mathrm{dia}}+H_{2DEG}+H_{\mathrm{int}}.$$

$$H_{\mathrm{cav}}=\sum _{q_x,n_z}{\omega }_{q_x,n_z}^{\mathrm{cav}}{a}_{q_x,n_z}^{†}{a}_{q_x,n_z}.$$

$$H_{\mathrm{dia}}=\sum _{q_x,n_z,n_z^{\prime }}{D}_{q_x}^{n_z,n_z^{\prime }}{a}_{q_x,n_z}\left(a_{-q_x,n_z^{\prime }}+a_{q_x,n_z^{\prime }}^{†}\right)+\text{mbox}h.c.,$$

$$D_{q_x}^{n_z,n_z^{\prime }}=\frac{2\pi n_{\mathrm{qw}}{n}_e{e}^2}{\epsilon L_z{m}^{\ast }\sqrt{{\omega }_{q_x,n_z}^{\mathrm{cav}}{\omega }_{q_x,n_z^{\prime }}^{\mathrm{cav}}}}\sin \frac{n_z\pi }{2}\sin \frac{n_z^{\prime }\pi }{2}.$$

$$H_{2DEG}=\sum _{q_x}\sum _{\ell \le \nu <\ell }{\omega }_{{\ell }^{\prime }\ell }^{\mathrm{LL}}{b}_{q_x,{\ell }^{\prime }\ell }^{†}{b}_{q_x,{\ell }^{\prime }\ell }.$$

$$H_{\mathrm{int}}=\sum _{q_x,n_z}\sum _{\ell \le \nu <{\ell }^{\prime }}i{\Omega }_{q_x,n_z}^{{\ell }^{\prime }\ell }(a_{q_x,n_z}+a_{-q_x,n_z}^{†})b_{q_x,{\ell }^{\prime }\ell }^{†}+\text{mbox}h.c.,$$

$${\Omega }_{q_x,n_z}^{{\ell }^{\prime }\ell }$$

$$p_{q_x}=W_{q_x}{a}_{q_x}+Y_{q_x}{a}_{-q_x}^{†}$$

Full text

I Introduction

The Dicke model, describing an ensemble of identical two-level systems (matter excitations) coupled to a single bosonic (cavity photon) mode, is a prototypical model of cavity quantum electrodynamics (QED) Dicke (1954). In the so-called ultrastrong coupling regime Ciuti et al. (2005); Forn-Díaz et al. (2019), when the coupling strength (Rabi frequency) becomes comparable to the energy splitting of the two-level system and that of the photon, the Dicke model model was shown to exhibit the so-called superradiant quantum phase transition (SQPT) towards a ground state characterised by a finite static average of the photon field Hepp and Lieb (1973); Emary and Brandes (2003). To the best of our knowledge, this phase transition has never been observed at equilibrium, although ultra-strong coupling regime has been reached in a two-dimensional electron gas (2DEG) in a semiconductor nanostructure placed in a cavity and subject to a perpendicular static magnetic field, so that the matter excitations were represented by the cyclotron resonance Scalari et al. (2012). Moreover, a softening of the lowest polaritonic excitation has been recently observed in this system Keller et al. (2018).

The Dicke model has an intrinsic flaw: it must be obtained by a reduction of a full microscopic model of some matter system coupled to the electromagnetic field, and typically, the assumptions used to justify this reduction, break down when the model is pushed to the ultrastrong coupling. As a consequence, the SQPT is usually prevented by the so-called “no-go” theorems Rzażewski et al. (1975); Nataf and Ciuti (2010a); Todorov and Sirtori (2012); Hayn et al. (2012); Chirolli et al. (2012); Bamba and Ogawa (2014); Rousseau and Felbacq (2017); Andolina et al. (2019). They express the simple fact that a physical system cannot respond to a static uniform vector potential which can be simply removed by a gauge tranformation. On top of the driven-dissipative scenario Dimer et al. (2007); Damanet et al. (2019), which has been successfully realized with cold atoms Baumann et al. (2010), different suggestions have been proposed to circumvent “no-go” theorems at equilibrium. These include systems with magnetic-dipole interactions due to the presence of cavity magnetic fields Knight et al. (1978) or its circuit QED analog with an inductive coupling Nataf and Ciuti (2010b, 2011) that can be of much larger magnitude Devoret et al. (2007). Notably, in the past two decades it has been shown in several works for different physical systems that upon a proper microscopic treatment, the mysterious SQPT assumes a more familiar shape of a ferroelectric Keeling (2007); Vukics et al. (2015) or an excitonic insulator Pellegrino et al. (2016); Mazza and Georges (2019) instability. In these studies, the crucial role of the Coulomb interaction has been pointed out. In addition, the instability occurred at length scales much shorter than the cavity size, thereby questioning the very role of the cavity.

In this work, we present a model without Coulomb interaction and still exhibiting a SQPT. Namely, we consider a 2DEG with Rashba spin-orbit coupling, placed inside an optical cavity, and subject to a perpendicular magnetic field $B$. In the decoupled 2DEG, the Landau levels can cross at certain values of $B$ corresponding to dipole-allowed excitations with zero energy. The presence of such intrinsic soft excitations greatly enhances the effect of the coupling to the transverse electromagnetic field. We develop a theory of Rashba cavity QED for integer filling factors and show that this coupling leads to further softening of the system and appearance of some “superradiant” phases. Crucially, the instability occurs at a finite wave vector of the cavity field; to describe it, all high-energy cavity modes must be included without any truncation in energy. This instability is of magnetostatic nature; the resulting “superradiant” phase is a remote relative of Condon domains of spontaneous magnetization known since long ago for bulk metals in a magnetic field Shoenberg (1984). Moreover, it turns out that this instability can also occur without the cavity: the coupling to the free vacuum field appears sufficient.

II The model

It is well known Weisbuch et al. (1992) that the effective strength of the light-matter coupling is enhanced if multiple copies of the material system are present. We therefore consider $n_{\mathrm{qw}}$ identical quantum wells, each hosting a 2DEG with the single-electron Hamiltonian containing a Rashba coupling term Bychkov and Rashba (1984),

[TABLE]

Here $\mathbf{p}=-i(\partial_{x},\partial_{y})$ is the 2D in-plane electron momentum (we set $\hbar=1$), $m^{*}$ is the effective mass, $\boldsymbol{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is the vector of Pauli matrices, and $\alpha$ is the Rashba spin-orbit coupling constant. Typically, for some existing InSb samples Becker et al. (2010), $m^{*}\simeq 0.02m_{0}$ ($m_{0}$ being the free electron mass), $\alpha\simeq 0.7\>\mbox{eV}\cdot\mbox{\AA}$. We assume $e>0$, the electron charge being $-e$. Finally, the vector potential $\mathbf{A}=\mathbf{A}^{\mathrm{ext}}+\mathbf{A}^{\mathrm{cav}}$ consists of two parts that we discuss separately.

First, $\mathbf{A}^{\mathrm{ext}}(\mathbf{r})=(-By,0,0)$ corresponds to the external magnetic field $B$, applied perpendicularly to the 2DEG plane (in the $z$ direction). The resulting single-particle spectrum consists of Landau levels (LLs) with energies

[TABLE]

where $\eta=0,1,2,\ldots$ is the LL index, $s=\pm 1$ for $\eta\geq 1$ and $s=1$ for $\eta=0$ is what remains of the spin index, $\omega_{c}=eB/(m^{*}c)$ is the cyclotron frequency, and $\gamma\equiv\alpha m^{*}l_{B}$ with $l_{B}\equiv(eB/c)^{-1/2}$ being the magnetic length. In Fig. 1, we show the LL energies $\epsilon_{\eta,s}$ for parameters consistent with InSb Becker et al. (2010); the spectrum exhibits crossings that occur between LLs $(\eta_{1},s)$ and $(\eta_{2},-s)$ satisfying the conditions Hernangómez-Pérez et al. (2013) $|\eta_{1}-\eta_{2}|>1$ and

[TABLE]

Note that levels with the same $s$ never cross. Each Landau level has a degeneracy $L_{x}L_{y}/(2\pi{l}_{B}^{2})$ where $L_{x}L_{y}$ is the sample area. We assume to be at zero temperature, at a fixed electron density $n_{e}$, and at a magnetic field corresponding to an integer filling factor $\nu\equiv 2\pi{l}_{B}^{2}n_{e}$. Indeed, the SQPT is associated with a change in the character of the non-degenerate ground state of a gapped system. Lifting the ground state degeneracy, which occurs at fractional fillings, represents a totally different problem.

The vector potential $\mathbf{A}^{\mathrm{cav}}(\mathbf{r})$ of the cavity field is defined by the mode expansion, determined by the cavity shape. For simplicity, here we consider a perfect metallic cavity whose dimensions satisfy $L_{x}\gg L_{z}\gg L_{y}$, filled by a material with a dielectric constant $\varepsilon$. Then, one can consider only resonator modes with wave vectors $\mathbf{q}=(q_{x},0,q_{z})$, where $q_{x}$ is varying continuously and $q_{z}=\pi n_{z}/L_{z}$, $n_{z}=1,2,3,\ldots$ The corresponding mode frequencies are $\omega^{\mathrm{cav}}_{q_{x},n_{z}}=(c/\sqrt{\varepsilon})\sqrt{q_{x}^{2}+q_{z}^{2}}$. The cavity vector potential then reads Kakazu and Kim (1994); Hagenmüller et al. (2010)

[TABLE]

where $a_{q_{x},n_{z}}^{\dagger}$ and $a_{q_{x},n_{z}}$ are the photon creation and annihilation operators and $\mathbf{u}_{y}$ is the unit vector in the $y$ direction. We assume the whole 2DEG sample with $n_{\mathrm{qw}}$ quantum wells to be much thinner than $L_{z}$ and placed in the middle of the cavity. Then, what enters Eq. (1), is $\mathbf{A}_{\mathrm{cav}}(z=L_{z}/2)$ and the modes with even $n_{z}$ are decoupled.

III Polariton modes and instability

The “superradiant” instability is signalled by the vanishing of the lowest polariton frequency. To find the polariton modes – the excitations of the coupled 2DEG-cavity system – one can proceed in several ways. For example, similarly to that adopted in Ref. Hagenmüller et al. (2010) for the same problem without spin-orbit coupling, one writes the 2DEG many-body Hamiltonian in terms of creation and annihilation operators for inter-LL excitations, which are approximately bosonic; then the full Hamiltonian of the 2DEG and the cavity becomes bilinear in the bosonic operators and thus can be diagonalized by the appropriate Bogoliubov transformation. Alternatively, one can write the (zero-temperature or Matsubara) action for coupled electron and photon fields, integrate out the electrons, and expand the resulting bosonic action to the second order in the photon vector potential. Both (rather standard) calculations are given in the respective appendices A and B, and their equivalence is checked explicitly.

As a result, the polariton frequencies are given by the solutions of the equation

[TABLE]

where $Q_{yy}(q_{x},\omega)$ is the susceptibility determining the linear response $j_{y}=Q_{yy}(\delta{A}_{y}/c)e^{iq_{x}x-i\omega{t}}$ of the 2D electron current density $j_{y}$ to a perturbing vector potential $\delta\mathbf{A}$ on top of $\mathbf{A}^{\mathrm{ext}}$ included in the unperturbed system, $\mathbf{A}=\mathbf{A}^{\mathrm{ext}}+\delta\mathbf{A}e^{iq_{x}x-i\omega{t}}$. The susceptibility consists of two contributions, the diamagnetic one and the sum over all inter-LL transitions in all quantum wells,

[TABLE]

Here the LL indices $(\eta,s)\equiv\ell$ are combined into a single label, ordered according to the LL energies $\epsilon_{\ell}$, Eq. (2), so that LLs with $\ell\leq\nu$ are filled, and those with $\ell>\nu$ are empty. The transition energy $\omega^{\mathrm{LL}}_{\ell^{\prime}\ell}\equiv\epsilon_{\ell^{\prime}}-\epsilon_{\ell}$, and the reduced coupling constants $g^{\ell^{\prime}\ell}_{q_{x}}$ (dipole matrix elements) are defined as

[TABLE]

where $\tan{\theta}_{l}=[-1+s\sqrt{1+8\eta\gamma^{2}}]/(\sqrt{8\eta}\gamma)$, and the overlap function $\Theta_{n_{2}}^{n_{1}}$ containing the $q_{x}$ dependence is given by

[TABLE]

with $L^{M-m}_{m}(\xi)$ the generalized Laguerre polynomial of $\xi\equiv l_{B}^{2}q_{x}^{2}/2$, $S=\operatorname{sign}\left[q_{x}(n_{2}-n_{1})\right]$, $m=\min\{n_{1},n_{2}\}$, and $M=\max\{n_{1},n_{2}\}$. Since at $q_{x}=0$ we have $\Theta_{n_{2}}^{n_{1}}=\delta_{n_{1},n_{2}}$, the reduced coupling constants $g^{\ell^{\prime}\ell}_{q_{x}=0}$ are non-zero only between consecutive LLs, $\eta^{\prime}=\eta\pm 1$, with no restriction on $s$. In contrast, at finite $q_{x}$, this selection rule is relaxed.

Equations (5)–(8) represent the main analytical result of this paper. Note that in Eq. (5) all information about the cavity is on the left-hand side, while all information about the 2DEG is on the right. At $\omega=0$ (i.e., at the sought superradiant transition) the left-hand side is proportional to $c^{2}$ which is much larger than any velocity scale occurring in a typical solid. Moreover, when $\omega=0$, the second term of the right-hand side of Eq. (6) is nothing but the total sum of the oscillator strengths $|g^{\ell^{\prime}\ell}_{q_{x}}|^{2}/\omega^{\mathrm{LL}}_{\ell^{\prime}\ell}$, which is the fundamental quantum optics quantity that determines the occurrence of the SQPT in multilevel systems coupled to cavity fields Hayn et al. (2012). It diverges at the level crossing, balancing the large $c^{2}$ factor in the left-hand side and allowing the existence of a solution to Eq. (5).

The key reason is that for the spin-flip transitions at $q_{x}\neq 0$, the dipoles $g^{\ell^{\prime},\ell}_{q_{x}}$ can be non-zero even at a crossing between $\epsilon_{\ell^{\prime}}$ and $\epsilon_{\ell}$, opening the possibility of a diverging oscillator strength. In Fig. 2, we present an example of such crossing. This is in sharp contrast with what happens at $q_{x}=0$, where gauge invariance demands the vanishing of the dipole between eigenstates of equal energy (i.e., at the crossings of energy levels), which is at the heart of many no-go theorems regarding the Dicke SQPT for spatially uniform cavity fields Rzażewski et al. (1975); Nataf and Ciuti (2010a); Todorov and Sirtori (2012); Hayn et al. (2012); Bamba and Ogawa (2014); Andolina et al. (2019).

Indeed, for a generic single-electron Hamiltonian $\mathcal{H}$, after the minimal coupling replacement of the electron momentum $\mathbf{p}\rightarrow\mathbf{p}+(e/c)\mathbf{A}_{\mathrm{cav}}$, where $\mathbf{A}_{\mathrm{cav}}$ is the uniform cavity field, the matrix element of the linear light-matter coupling term between two arbitrary eigenstates $|1\rangle$, $|2\rangle$ of $\mathcal{H}$, is proportional to that of the electron velocity $\mathbf{v}=\partial\mathcal{H}/\partial\mathbf{p}=i[\mathcal{H},\mathbf{r}]$, so that $\mathbf{v}_{12}=\langle 1|[\mathcal{H},\mathbf{r}]|2\rangle=(\epsilon_{1}-\epsilon_{2})\langle 1|\mathbf{r}|2\rangle$. Then, at the crossing $\epsilon_{1}-\epsilon_{2}$ vanishes, and both the velocity matrix element $\mathbf{v}_{12}$ and the oscillator strength $\propto|\mathbf{v}_{12}|^{2}/(\epsilon_{1}-\epsilon_{2})$ vanish. Crucially, this argument does no apply to spatially non-uniform fields.

From a different perspective, gauge invariance imposes a constraint on $Q_{yy}(q_{x},\omega)$: a physical quantity (current) cannot respond to a static spatially homogeneous vector potential, thus $Q_{yy}(0,0)=0$. This prohibits the instability at $q_{x}=0$; models or approximations violating this constraint can give wrong results. $Q_{yy}(0,0)=0$ implies a cancellation between the two terms in Eq. (6). This cancellation is usually ensured by sum rules such as the famous Thomas-Reiche-Kuhn strength of localized atomic systems Rzażewski et al. (1975); Nataf and Ciuti (2010a); Todorov and Sirtori (2012); Chirolli et al. (2012); Hayn et al. (2012); Bamba and Ogawa (2014); Andolina et al. (2019). Here, we have checked it numerically.

In Fig. 3(a), we plot the two sides of Eq. (5) at $\omega=0$ as functions of the cavity field wave-vector $q_{x}$ and show six values of $q_{x}>0$ (shown in magenta dashed lines) when $\omega=0$ is a solution. In Fig. 3(b), we display the lowest polariton frequency $\omega_{q_{x}}^{\mathrm{pol}}$ which indeed vanishes at the indicated values of $q_{x}$. Between two values of $q_{x}$ where Eq. (5) is satisfied at $\omega=0$, the left-hand side of Eq. (5) is smaller than the right-hand side and the system is in the “superradiant” phase. In the inset of Fig. 3(a), we plot the two sides of Eq. (5) at $\omega=0$ as a function of $B$ for a given value of $q_{x}$, for $n_{\mathrm{qw}}=1$ and $n_{\mathrm{qw}}=10^{3}$. Essentially, we exploit a divergence appearing around level crossings to make Eq. (5) have a solution, and then use large $n_{\mathrm{qw}}$ to have an extended regime for the superradiant phase.

Fig. 4 shows the phase diagram in the plane $(\alpha,q_{x})$, for fixed magnetic field $B$ and filling $\nu$. The “superradiant” phase appears for nonzero $q_{x}$ and for values of $\alpha$ very close to those given by Eq. (3) for some integers $\eta_{1}\neq\eta_{2}$ satisfying $\eta_{1}+\eta_{2}=\nu$, and depicted in dashed red lines in Fig. 4. The characteristic width $\Delta\alpha$ of the “superradiant” regions on the phase diagram can be estimated as (see Appendix C)

[TABLE]

and the typical scale of $q_{x}$ is given by the inverse cyclotron radius, $(\sqrt{\nu}\,l_{B})^{-1}$. The “superradiant” regions are very narrow; this happens because the mechanism for the instability can be traced to the magnetostatic interaction, as we discuss below. Typically, one arrives at the Dicke model assuming the light-matter coupling via the cavity electric field. However, this electric field, $\propto\partial\mathbf{A}_{\mathrm{cav}}/\partial{t}$, vanishes at $\omega=0$. The remaining magnetic interaction is intrinsically weak. These simple physical arguments are not obvious from the equations.

From Eq. (9) and Fig. 4 we see that small filling factors are favoring the “superradiant” phase. This is in stark contrast to the condition of $\nu\gg 1$ formulated in Ref. Hagenmüller et al. (2010) to achieve the ultrastrong coupling regime. Again, the reason for this difference is that the SQPT obtained here is determined by the magnetic coupling and not the electric one.

IV SQPT as a magnetostatic instability

Consider the 2DEG in free space (no cavity). When placed in a static homogeneous magnetic field $\mathbf{B}^{\mathrm{ext}}(\mathbf{r})=B\mathbf{u}_{z}$, it develops a static equilibrium magnetization $\mathbf{M}^{\mathrm{eq}}(\mathbf{r})=M^{\mathrm{eq}}\mathbf{u}_{z}\,\delta(z-L_{z}/2)$. Let us see at what conditions the system can spontaneously develop an additional inhomogeneous magnetization $\delta\mathbf{M}(\mathbf{r})$ which would produce an inhomogeneous magnetic field. Equivalently, we study the stability of the described equilibrium configuration with respect to a small perturbation, $\mathbf{B}=\mathbf{B}^{\mathrm{ext}}+\delta\mathbf{B}(\mathbf{r})$.

We look for static magnetic field configurations $\delta\mathbf{B}(\mathbf{r})$ which minimize the classical free energy functional

[TABLE]

under the constraint $\boldsymbol{\nabla}\cdot\delta\mathbf{B}=0$. The first line contains the energy density of the free field as well as the contribution due to $\mathbf{M}^{\mathrm{ext}}$, the magnetization of the currents which produce $\mathbf{B}^{\mathrm{ext}}$. The second line is the free energy of the 2DEG.

Its quadratic part is determined by the static magnetic susceptibility $\chi_{ij}(\mathbf{r},\mathbf{r}^{\prime})$, which gives the response of the magnetization $\delta\mathbf{M}$ to a magnetic field perturbation $\delta\mathbf{B}$ on top of $\mathbf{B}^{\mathrm{ext}}$. Since $\mathbf{B}=\boldsymbol{\nabla}\times\mathbf{A}$, and the static current density $\mathbf{j}$ can be written in terms of the magnetization $\mathbf{M}$ as $\mathbf{j}=c\boldsymbol{\nabla}\times\mathbf{M}$, the susceptibility $\chi$ is related to the current response $Q(\omega=0)$ found earlier on,

[TABLE]

Here $\delta(y-y^{\prime})$ appears because our starting model was restricted to a $y$-independent $\mathbf{A}_{\mathrm{cav}}(\mathbf{r})$, see Eq. (II).

Varying $F+\int\lambda(\mathbf{r})\,\boldsymbol{\nabla}\cdot\mathbf{B}\,d^{3}\mathbf{r}$ with the Lagrange multiplier $\lambda(\mathbf{r})$, we arrive at the equations

[TABLE]

In Eq. (12a), the terms containing $\mathbf{B}^{\mathrm{ext}}$ and $\mathbf{M}^{\mathrm{ext}}$ cancel out by construction; the terms with $\mathbf{M}^{\mathrm{eq}}$ drop out because a magnetization $\propto\mathbf{u}_{z}\,\delta(z-L_{z}/2)$ does not produce any magnetic field except near the sample edges; thus, Eq. (12a) becomes a linear homogeneous equation for $\delta\mathbf{B}(\mathbf{r})$. The existence of a non-trivial solution corresponds to a direction along which the quadratic form in Eq. (IV) is flat.

Let us try the modulation

[TABLE]

with some $q_{x}>0$. The solution of Eq. (12a) for $\delta\mathbf{B}$,

[TABLE]

when substituted in Eq. (12b), gives the magnetostatic instability condition

[TABLE]

This condition is identical to Eq. (5) at $\omega=0$ in the limit $L_{z}\to\infty$, that is, in the absence of the cavity.

The obtained modulational instability of the magnetization is long known for bulk paramagnetic materials in an external magnetic field, where it gives rise to the so-called Condon domains of different magnetization Shoenberg (1984). Indeed, such a domain structure has the characteristics of the “superradiant” phase: the magnetic field produced by the domains corresponds to a non-zero expectation value of the photon field $\langle{a}_{q_{x},n_{z}}\rangle$. On the other hand, while the superradiant phase in the Dicke model spontaneously breaks the $Z_{2}$ symmetry associated with the sign of the spontaneously generated uniform cavity field Emary and Brandes (2003) (variations of the model with higher symmetries have also been considered Hepp and Lieb (1973); Nataf et al. (2012); Baksic and Ciuti (2014)), the domain structure obtained here spontaneously breaks at least the continuous translation symmetry. In particular, for a sinusoidal modulation the $Z_{2}$ sign flip is equivalent to translating the pattern by half a period. To make precise statements about the broken symmetries, one would need a detailed quantitative study of the “superradiant” phase and of the resulting magnetization profile; such a study is beyond the scope of the present paper.

V Conclusions and outlook

We have proved that the SQPT can be reached in a cavity QED system with Rashba spin-orbit coupling and non-uniform cavity resonator fields. Within our model, the appearance of the SQPT is a consequence of the singularity of the spin-flip transitions, for which the spin-flip dipole at the field wave-vector $q_{x}\neq 0$ can be non-zero even if the transition energy vanishes. Consequently, the SQPT must occur close to LLs energy crossings (see Figs. 3 and 4), and requires relatively fine tuning of the magnetic field and electron density, as well as large $n_{\mathrm{qw}}$. We have shown that the SQPT can be viewed as a magnetostatic instability of the 2DEG, the static cavity field corresponding to the magnetic field induced by a spatially modulated 2DEG magnetization. Moreover, the presence of the cavity is not necessary : it can also happen via the coupling to the free vacuum field.

The microscopic model we studied in this paper has the minimal number of ingredients necessary to produce SQPT. To make a connection with state of the art experiments Scalari et al. (2012); Keller et al. (2018), several other ingredients must be introduced. Zeeman coupling of the 2DEG to the magnetic filed is likely to enhance the effect since the spin contribution to the magnetic susceptibility is usually paramagnetic. Coulomb interaction is also likely to further soften the excitations due to the excitonic effect. A very important ingredient is the disorder which lifts the LL degeneracy and broadens the cyclotron resonance. Coherent state based methods have been quite successful in describing the local density of states in a 2DEG with smooth disorder in a strong magnetic field Champel and Florens (2009); Hernangómez-Pérez et al. (2013); the effect of smooth disorder on the inter-LL transitions remains an open problem. Effects of strain, as well as mixing between bands with different spins (through, for instance, the multiband Luttinger $6\times 6$ $\mathbf{k}\cdot\mathbf{p}$ model Luttinger (1956)), which both importantly impact the amplitude and the nature of the Rashba coupling Winkler et al. (1996); Winkler (2000); Moriya et al. (2014), could also be studied. Finally, adapting our calculations to some other cavities, like the Split Ring Resonator (SSR) and Complementary SRR Maissen et al. (2014), where an additionnal geometric factor can enhance the light matter interaction, could also be useful. Finally, in this paper, we focused on the instability, leaving aside the study of the “superradiant” phase itself, a topic that deserves future investigation as well. Moreover, the new physical ingredients listed above, will also be important in determining the properties of the “superradiant” phase.

Acknowledgements.

We acknowledge fruitful discussions with Tobias Wolf, Matteo Biondi, Jérome Faist, Giacomo Scalari, Janine Keller, Maksym Myronov, Rai Moriya, Takaaki Koga, Roland Winkler, Cristiano Ciuti and David Hagenmüller.

Appendix A Bogoliubov transformation approach

The single-electron eigenstates of Hamiltonian (1) with the external vector potential $\mathbf{A}_{\mathrm{ext}}=(-By,0,0)$ are labeled by $\ell\equiv(\eta,s)$ and $p_{x}$, the momentum in the $x$ direction, taken as an integer multiple of $2\pi/L_{x}$, with $L_{x}$ being the sample length the $x$ direction. It is convenient to order $\ell$ according to the energies $\epsilon_{\ell}$, given by Eq. (2). The spinor wave functions of the eigenstates are

[TABLE]

with $\tan{\theta}_{\eta{s}}=-u_{\eta}+s\sqrt{1+u_{\eta}^{2}}$, $u_{\eta}=1/(\sqrt{8\eta}\gamma)$ and the harmonic oscillator wave functions

[TABLE]

where $H_{\eta}$ is the Hermite polynomial of degree $\eta$. The many-body ground state of the 2DEG without coupling to the cavity is

[TABLE]

where $|0\rangle$ is the vacuum, $c^{{\dagger}}_{j,\ell,p_{x}}$ is the fermionic creation operator for an electron on the level $\ell$ with momentum $p_{x}$ in the $j$th quantum well. The product over $p_{x}$ goes over all $L_{x}L_{y}/(2\pi{l}_{B}^{2})$ allowed values. Cavity photons induce collective electronic excitations which can be described by the approximatively bosonic bright modes creation operator, defined for $\ell,\ell^{\prime}$ such that $\epsilon_{\ell}\leq\epsilon_{\nu}<\epsilon_{\ell^{\prime}}$:

[TABLE]

where the sum runs over the $n_{\mathrm{qw}}L_{x}L_{y}/(2\pi{l}_{B}^{2})$ individual excitation between filled and empty LLs $\ell$ and $\ell^{\prime}$, respectively.

The full light-matter Hamiltonian $H$ contains four terms:

[TABLE]

The free cavity photon part is

[TABLE]

The diamagnetic part of the light-matter interaction arises from the $\mathbf{A}_{\mathrm{cav}}^{2}$ term in the single-particle Hamiltonian (1):

[TABLE]

$H_{\mathrm{2DEG}}$ stands for the electronic part of the Hamiltonian, which can be written using the operators in (19), as:

[TABLE]

The corresponding energy differences $\omega^{\mathrm{LL}}_{\ell^{\prime}\ell}=\epsilon_{\ell^{\prime}}-\epsilon_{\ell}$ are associated to each transition across the Fermi level. Finally, the linear in $\mathbf{A}_{\mathrm{cav}}$ term produces

[TABLE]

where the $\ell,\ell^{\prime}$ sum is again over $\epsilon_{\ell}\leq\mu\leq\epsilon_{\ell^{\prime}}$, and the Rabi frequencies for each transition are determined by matrix elements of $e^{iq_{x}x}(-i\partial_{y}/m^{*}-\alpha\sigma_{x})$ between different Landau level states:

[TABLE]

and the reduced coupling constants (dipole matrix elements) $g^{\ell^{\prime}\ell}_{q_{x}}$ are defined in Eq. (III).

In the following, we prove that the eigenfrequencies of the quadratic light-matter Hamiltonian $H$, shown in (20) can vanish for some wave vectors $q_{x}$. Looking for the bosonic magnetopolariton modes in the form

[TABLE]

which satisfies $[p_{q_{x}},H]=\omega_{q_{x}}^{\mathrm{pol}}p_{q_{x}}$, one has to diagonalize the corresponding Bogoliubov matrix $\mathcal{M}_{q_{x}}$. The polaritonic frequencies $\omega_{q_{x}}^{pol}$ are the positive eigenvalues of $\mathcal{M}_{q_{x}}$, given by the solutions of the equation $\det(\mathcal{M}_{q_{x}}-\omega)=0$. The determinant can be written as

[TABLE]

where we defined $\mathcal{I}_{q_{x}}(\omega)$ as

[TABLE]

Using the fact that

[TABLE]

we can evaluate the $n_{z}$ sum in Eq. (27) explicitly:

[TABLE]

Equating to zero the first prefactor in Eq. (27), we obtain the equation for the polariton frequencies, the main analytical result of the present paper:

[TABLE]

which is equivalent to Eqs. (5) and (6) of the main text.

In the numerical calculations, we diagonalized the Bogoliubov matrix in a truncated basis of LLs and cavity modes, and checked for convergence. Typically, this required 50 LLs and $n_{z}$ up to $10^{4}$ for a given $q_{x}$.

Appendix B Effective action approach

Let us describe the system by its Euclidean action $S$ (weight $e^{-S}$) in the imaginary time $\tau$ varying in the interval $0\leq\tau\leq\beta$, where $\beta=1/T$ is the inverse temperature that we will eventually send to infinity ($T\to 0$). The action consists of two parts, $S=S_{\mathrm{cav}}+S_{\mathrm{2DEG}}$. The action of the free cavity field is

[TABLE]

Here $\mathbf{A}^{\mathrm{cav}}$ is the vector potential of the cavity field in the Coulomb gauge, $\boldsymbol{\nabla}\cdot\mathbf{A}^{\mathrm{cav}}=0$, while $\varepsilon$ is the dielectric constant of the material filling the cavity. The electrons of the 2DEG, whose number is assumed to be fixed, are described by the action

[TABLE]

where $\psi(x,y)$, $\psi^{\dagger}(x,y)$ are $2n_{\mathrm{qw}}$-component vectors of Grassmann variables, corresponding to two electron spin projections in each of the $n_{\mathrm{qw}}$ quantum wells. The single-electron Hamiltonian $\mathcal{H}_{\mathrm{2DEG}}$ is given by Eq. (1) with $\mathbf{A}=\mathbf{A}^{\mathrm{ext}}+\mathbf{A}^{\mathrm{cav}}$. $\mathcal{H}_{\mathrm{2DEG}}$ a $2\times 2$ matrix in the spin space, and in Eq. (33) it should be replicated $n_{\mathrm{qw}}$ times which is indicated by the superscript $\otimes{n}_{\mathrm{qw}}$.

Let us integrate out the electrons and obtain the effective action $S_{\mathrm{eff}}[\mathbf{A}^{\mathrm{cav}}]$ for the field only. It is convenient to pass to Matsubara frequencies

[TABLE]

where $\omega_{n}=2\pi{n}T$, and $\mathbf{A}^{\mathrm{cav}}(-\omega_{n})=[\mathbf{A}^{\mathrm{cav}}(\omega_{n})]^{*}$. Then the effective action can be written as a series in powers of $\mathbf{A}^{\mathrm{cav}}$:

[TABLE]

where the kernels $\mathcal{Q}^{(n)}$ are given by the sum of electronic loops where the vertices $\psi^{*}\psi\mathbf{A}^{\mathrm{cav}}$ and $\psi^{*}\psi(\mathbf{A}^{\mathrm{cav}})^{2}$ appear $n_{1}$ and $n_{2}$ times, respectively, and $n_{1}+2n_{2}=n+1$. Generally speaking, the system is allowed to have an equilibrium current density, $\mathbf{j}^{\mathrm{eq}}(\mathbf{r})$, since it is placed in an external magnetic field and time reversal invariance is broken; in our specific case $\mathbf{j}^{\mathrm{eq}}=0$.

Analytical continuation in frequency of the kernels $\mathcal{Q}^{(n)}$ from the positive imaginary semiaxis $\omega=i\omega_{n}$, $\omega_{n}>0$, to the real axis $\omega$, gives the response functions which determine the response of the current density $\delta\mathbf{j}$ to a change in the vector potential, $\mathbf{A}=\mathbf{A}^{\mathrm{ext}}+\delta\mathbf{A}$; in the linear order,

[TABLE]

The polariton modes of the coupled system are determined by the quadratic part of $S_{\mathrm{eff}}$ whose kernel is

[TABLE]

Namely, the polariton frequencies are solutions of the equation $\det\mathcal{D}^{-1}(\omega)=0$, where $\mathcal{D}^{-1}(\omega)$ is the analytical continuation of $\mathcal{D}^{-1}(i\omega_{n})$ from the positive imaginary semiaxis and the determinant should be understood in the operator sense. Equivalently, the polariton eigenmodes are the nonzero solutions of

[TABLE]

with $\boldsymbol{\nabla}\cdot\mathbf{A}=0$, which is just the third Maxwell’s equation with a source current. The SQPT corresponds to the appearance of an unstable direction in the quadratic form with the kernel $\mathcal{D}^{-1}(\omega=0)$.

It should be noted that the zero-frequency limit of the Kubo susceptibility $\mathcal{Q}(\omega)$, which was obtained by the analytical continuation of the Matsubara susceptibility $\mathcal{Q}(i\omega_{n})$ from the positive imaginary semiaxis $\omega_{n}>0$, is, generally speaking, different from the value of the Matsubara susceptibility $\mathcal{Q}(i\omega_{n})$ taken directly at $\omega_{n}=0$. The Kubo limit $\mathcal{Q}(\omega\to{0})$ corresponds to evaluating the current with LL wave functions perturbed by a static $\delta\mathbf{A}$ and keeping the LL populations unperturbed; the Matsubara limit $\mathcal{Q}(i\omega_{n}=0)$ also takes into account the change in the LL populations which occurs because of the shift of the LL energies while the temperature and the chemical potential are kept fixed. In other words, they correspond to different order of limits $\omega\to 0$ and population relaxation time to infinity. In our case, we are working at zero temperature and fixed electron density, rather than at fixed chemical potential. Thus, we take the Kubo limit.

By gauge invariance, a static vector potential $\mathbf{A}(\mathbf{r})$ can affect observable quantities only via the associated magnetic field $\mathbf{B}=\boldsymbol{\nabla}\times\mathbf{A}$. Also, at zero frequency the continuity equation for the current density requires $\boldsymbol{\nabla}\cdot\mathbf{j}=0$, so one can write $\mathbf{j}=c\boldsymbol{\nabla}\times\mathbf{M}$, where $\mathbf{M}$ is called magnetization. Then it is convenient to express the response function $\mathcal{Q}_{ij}(\mathbf{r},\mathbf{r}^{\prime},\omega=0)$ in terms of the static magnetic susceptibility $\chi_{ij}(\mathbf{r},\mathbf{r}^{\prime})$, which determines the response of the magnetization $\delta\mathbf{M}$ to a magnetic field perturbation $\delta\mathbf{B}$ on top of $\mathbf{B}^{\mathrm{ext}}$. Then,

[TABLE]

where $e_{ikl}$ is the Levi-Civita antisymmetric tensor.

Now we proceed to evaluation of the response function $Q_{kl}(\mathbf{r},\mathbf{r}^{\prime},\omega)$ for the 2DEG with the Rashba spin-orbit coupling under an external magnetic field. The single-electron 2D current density operator, corresponding to Hamiltonian (1) for a single quantum well is given by

[TABLE]

where $\mathbf{u}_{x,y,z}$, denotes the unit vectors in the respective directions. As $\mathbf{A}=\mathbf{A}^{\mathrm{ext}}+\delta\mathbf{A}$, the term proportional to $\delta\mathbf{A}$ gives the diamagnetic contribution to the response, while the rest should be plugged in the Kubo formula. We neglect the 2DEG thickness compared to all other length scales, so all quantum wells contribute additively and the response has the form

[TABLE]

We are interested in the $Q_{yy}$ component at $\mathbf{q}=q_{x}\mathbf{u}_{x}$:

[TABLE]

where the first term comes from the Kubo formula and the last is the diamagnetic one.

Evaluation of the current matrix elements between the LL eigenstates (16) gives

[TABLE]

where $g^{\ell^{\prime}\ell}_{q_{x}}$ is given by Eq. (III). As a result,

[TABLE]

where $\mathcal{I}_{q_{x}}(\omega)$ is given by Eq. (28).

Equation (38) becomes

[TABLE]

with the boundary conditions $A_{y}(0)=A_{y}(L_{z})=0$. It has one family of solutions, $A_{y}(z)\propto\sin(\pi{n}_{z}z/L_{z})$ with even $n_{z}$, corresponding to the cavity modes which remain decoupled from the 2DEG since $A_{y}(L_{z}/2)=0$. The other family can be represented as

[TABLE]

It has a jump of the derivative at $z=L_{z}/2$, which should be matched to the right-hand side of Eq. (45), yielding

[TABLE]

which is equivalent to Eq. (5) of the main text and Eq. (31).

Appendix C Estimate for the “superradiant” phase width

Let us choose some LL crossing, $\epsilon_{\eta_{1},+}=\epsilon_{\eta_{2},-}$, and focus on $\alpha$ close to that given by Eq. (3) which is denoted by $\alpha_{\mathrm{c}}$. The LL indices $\eta_{1},\eta_{2}$ satisfy $\eta_{1}+\eta_{2}=\nu$, and we denote $\eta_{1}-\eta_{2}=\delta$. We will assume $\nu\gg{1}$, and the order-of-magnitude estimate is expected to be valid for $\nu\sim{1}$ as well. The dimensionless Rashba coupling strength corresponding to the crossing is then

[TABLE]

When $\alpha$ is detuned away from the crossing, keeping both the magnetic field and the filling constant, the LL energy difference becomes

[TABLE]

To estimate the coupling strength $g_{q_{x}}$ at the crossing, we use the well-known asymptotic expression for the generalized Laguerre polynomials with large index in terms of the Bessel function of the first kind, $J$, which gives

[TABLE]

Next, we need to find the angles $\theta_{\eta_{1},+}$ and $\theta_{\eta_{2},-}$. To the leading order in $1/\sqrt{\nu}$,

[TABLE]

and thus $\theta_{\eta_{2},-}=\theta_{\eta_{1},+}+\pi/2+O(1/\nu)$. As a result, the leading terms in the coupling $g_{q_{x}}$ cancel. This happens because of approximate orthogonality of the spin part of the wave functions for $|\eta_{1}-\eta_{2}|\ll\eta_{1},\eta_{2}$ and $s_{1}=-s_{2}$.

Expansion to the next order will generate many terms; the resulting bulky expressions are not very informative. To obtain an order-of-magnitude estimate, it is sufficient to notice that the coupling strength can be written in the form

[TABLE]

where for $\delta\sim{1}$ the function $\mathcal{G}_{\delta}(x)$ decays faster than $x$ at $x\to{0}$, reaches the first maximum at $x\sim{1}$, and the width of this maximum is also $\sim{1}$. Then, assuming that $Q_{yy}(q_{x},\omega=0)$ is dominated by a single term corresponding to the LL crossing under consideration and setting $\coth(q_{x}L_{z}/2)\to{1}$, we can write Eq. (5) of the main text as

[TABLE]

This equation has no solutions for $q_{x}$ if $\Lambda\gg{1}$; solutions appear when $\Lambda\sim{1}$ [Eq. (9) of the main text]; then the typical scale of $q_{x}$ and the width of the unstable $q_{x}$ interval are determined by $\sqrt{\nu}{q}_{x}l_{B}\sim{1}$.