††thanks: These three authors contributed equally.††thanks: These three
authors contributed equally.††thanks: These three authors contributed equally.

# Photon condensation, Van Vleck paramagnetism, and chiral cavities

Alberto Mercurio  Dipartimento di Scienze Matematiche e Informatiche, Scienze
Fisiche e Scienze della Terra, Università di Messina, I-98166 Messina, Italy
Gian Marcello Andolina  ICFO-Institut de Ciències Fotòniques, The Barcelona
Institute of Science and Technology, Av. Carl Friedrich Gauss 3, 08860
Castelldefels (Barcelona), Spain  Francesco M. D. Pellegrino  Dipartimento di
Fisica e Astronomia “Ettore Majorana”, Università di Catania, Via S. Sofia 64,
I-95123 Catania, Italy  INFN, Sez. Catania, I-95123 Catania, Italy  CNR-IMM,
Via S. Sofia 64, I-95123 Catania, Italy  Omar Di Stefano  Dipartimento di
Scienze Matematiche e Informatiche, Scienze Fisiche e Scienze della Terra,
Università di Messina, I-98166 Messina, Italy  Pablo Jarillo-Herrero
Department of Physics, Massachusetts Institute of Technology, Cambridge,
Massachusetts, USA  Claudia Felser  Max Planck Institute for Chemical Physics
of Solids, Nöthnitzer Str. 40, Dresden 01187, Germany  Frank H. L. Koppens
ICFO-Institut de Ciències Fotòniques, The Barcelona Institute of Science and
Technology, Av. Carl Friedrich Gauss 3, 08860 Castelldefels (Barcelona), Spain
ICREA-Institució Catalana de Recerca i Estudis Avançats, Passeig de Lluís
Companys 23, 08010 Barcelona, Spain  Salvatore Savasta  Dipartimento di
Scienze Matematiche e Informatiche, Scienze Fisiche e Scienze della Terra,
Università di Messina, I-98166 Messina, Italy  Marco Polini  Dipartimento di
Fisica dell’Università di Pisa, Largo Bruno Pontecorvo 3, I-56127 Pisa, Italy
Istituto Italiano di Tecnologia, Graphene Labs, Via Morego 30, I-16163 Genova,
Italy  ICFO-Institut de Ciències Fotòniques, The Barcelona Institute of
Science and Technology, Av. Carl Friedrich Gauss 3, 08860 Castelldefels
(Barcelona), Spain

###### Abstract

We introduce a gauge-invariant model of planar, square molecules coupled to a
quantized spatially-varying cavity electromagnetic vector potential
$\hat{\bm{A}}({\bm{r}})$. Specifically, we choose a temporally chiral cavity
hosting a uniform magnetic field $\hat{\bm{B}}$, as this is the simplest
instance in which a transverse spatially-varying $\hat{\bm{A}}({\bm{r}})$ is
at play. We show that when the molecules are in the Van Vleck paramagnetic
regime, an equilibrium quantum phase transition to a photon condensate state
occurs.

Introduction.—Basic quantum statistical mechanics dictates that free massless
photons, despite being bosonic particles, cannot condense into a single
macroscopic state [1]. On the other hand, when photons are confined to a
cavity [2, 3, 6, 7, 4, 5] and coupled to matter degrees of freedom (such as
excitons), condensation into a single quantum state can occur [8].

In this Letter, we are interested in the ground state of a photon condensate,
i.e. a state containing a macroscopically large number of coherent photons,
i.e. $\langle\hat{a}\rangle\propto\sqrt{N}$, where $\hat{a}$
($\hat{a}^{\dagger}$) destroys (creates) a cavity photon. Such an equilibrium
state was predicted in 1973 by Hepp and Lieb [9] and Wang and Hioe [10],
independently. These authors pointed out that, in the thermodynamic
$N,V\to\infty$ limit and for sufficiently strong light-matter interactions,
the celebrated Dicke model [11]—which describes $N$ two-level systems coupled
to a spatially uniform vector potential
$\hat{\bm{A}}=A_{0}{\bm{u}}(\hat{a}+\hat{a}^{\dagger})$, with $A_{0}\propto
1/\sqrt{V}$—has a finite-temperature second-order equilibrium phase transition
between a normal and a photon condensate state. To the best of our knowledge,
this phase transition has never been observed [12]. This is in agreement with
a large body of “converging” literature [13, 14, 19, 15, 16, 17, 18, 20, 21,
22, 23, 24, 25, 26] stating that photon condensation in light-matter
interacting systems is forbidden by gauge invariance when the diamagnetic
term, i.e. a term proportional to $\hat{\bm{A}}^{2}$, is included into the
theory. Such a term is missing in the Dicke model and is necessary for
defining any properly gauge-invariant model of light-matter interactions at
equilibrium.

\begin{overpic}[width=433.62pt]{3d_and_phase.pdf}
\put(-2.0,37.0){\footnotesize{(a)}}\put(45.0,37.0){\footnotesize{(b)}}
\end{overpic} Figure 1: (Color online) (a) Pictorial representation of the
setup studied in this Letter. Many planar square molecules are coupled to a
single-mode cavity hosting a quantized magnetic field $\hat{\bm{B}}$ directed
perpendicular to the plane where each molecule lies. Note that, in each
molecule, one of the hopping integrals (yellow) is different from the other
three ones. (b) Phase diagram of the model. Results in this figure have been
obtained by setting with $t=-\omega_{\rm c}$ and $\lambda=0.3$. The photon
condensate order parameter $\tilde{\alpha}$ (color scale) is plotted as a
function of the two microscopic parameters $\tau\in[t,-t]$ and
$\Theta\in[0,\pi[$.

Recent studies, however, have shown that when extended electron systems are
coupled in a gauge-invariant fashion to a transverse spatially-varying
electromagnetic vector potential $\hat{\bm{A}}({\bm{r}})$, equilibrium photon
condensation may occur as a Condon domain instability [27, 28, 29, 30, 31].
Key to this phenomenon is the orbital paramagnetic character of the itinerant
electron system. Orbital paramagnetism, i.e. a positive sign of the orbital
magnetic susceptibility, is a rather rare phenomenon in nature, as
diamagnetism tends to dominate. For example, we know that the orbital response
of a free parabolic-band electron gas has a diamagnetic character (Landau
diamagnetism) [1]. Similarly, the zero-energy Landau level in two-dimensional
(2D) systems of massless Dirac fermions is responsible for the strong
diamagnetic nature of graphite and graphene (McClure diamagnetism) [32].
Nevertheless, itinerant electron systems may display orbital paramagnetism
[33, 34, 35, 36, 37, 38].

Since Refs. [27, 28, 29] are limited to extended electron systems, the key
question we try and answer in this Letter is the following: can photon
condensation occur when a molecular system, which is not translationally
invariant, is coupled to a cavity field? We show that the answer is
affirmative provided that two conditions are satisfied. On the one hand, one
needs to work with a molecular system whose orbital response displays
paramagnetic character. To this end, we exploit a genuinely quantum mechanical
mechanism leading to orbital paramagnetism in a molecular system, which is
often dubbed “Van Vleck paramagnetism” [39]. This mechanism guarantees a
paramagnetic orbital response provided that the molecule has a doublet of
quasi-degenerate levels. On the other hand, one needs to confine these Van
Vleck molecules to a cavity where a transverse spatially-varying
${\bm{A}}({\bm{r}})$ is at play. Since a spatially-varying ${\bm{A}}$ leads to
a finite magnetic field
${\bm{B}}({\bm{r}})=\nabla_{\bm{r}}\times{\bm{A}}({\bm{r}})$, the simplest
cavity one can consider is a cavity hosting a uniform ${\bm{B}}$ field. Since
such field breaks time-reversal symmetry, our cavity falls into the category
of chiral cavities [40] with chirality of temporal (rather than structural)
character.

Model.—We consider a single, spinless electron [41] hopping between the
$n_{\rm s}$ sites of a planar plaquette, such as a square or triangular
plaquette, lying in the $x$-$y$ plane. This system will be considered as a toy
model of a “molecule” since the plaquette contains a small number of sites
($n_{\rm s}=4$ and $n_{\rm s}=3$ below) and the electron roaming in the
plaquette is therefore far away from the lattice regime where electronic Bloch
bands emerge.

We then consider a system of $N$ such molecules described by the following
electronic Hamiltonian, $\hat{\mathcal{H}}_{\rm e}=\sum^{N}_{k=1}\hat{h}_{{\rm
e},k}$ , where

$\hat{h}_{{\rm e},k}=\sum_{j=0}^{n_{\rm
s}-1}\left(t_{j,j+1}\hat{c}_{j,k}^{\dagger}\hat{c}_{j+1,k}+{\rm H.c.}\right)$
(1)

is the single-molecule Hamiltonian. Here, $\hat{c}_{j,k}^{\dagger}$
($\hat{c}_{j,k}$) creates (destroys) an electron on the $j$-th site of the
$k$-th plaquette, $t_{j,j+1}=\delta_{j,0}t^{\prime}+(1-\delta_{j,0})t$ is the
tunnelling amplitude between neighboring sites of a plaquette, and

$t^{\prime}=\tau e^{-i\Theta}~{},$ (2)

where $\tau\in[t,-t]$ (with $t<0$ is assumed to be real, without any loss of
generality) and $\Theta\in[0,\pi[$ are the tunneling “intensity” and phase,
respectively. We emphasize that we have taken $t^{\prime}\neq t$ to simulate
different physical conditions. For example, letting $\tau\to 0$ physically
corresponds to opening up the closed plaquette, breaking the $0$-$1$ link. For
$\Theta\neq 0$, the Hamiltonian is not invariant under time-reversal symmetry,
and the plaquette supports a circulating ground-state current. Viceversa, when
$\Theta=0$ no ground-state currents can flow. Finally, we note that a finite
value of $\Theta$ can be physically obtained by applying an external,
classical magnetic field along the $z$ direction.

We now couple the molecular system to a single-mode cavity with a non-
vanishing magnetic field oriented along the $z$ direction. The dynamics of the
total system is governed by the following Hamiltonian

$\displaystyle\hat{\mathcal{H}}=\omega_{\rm
c}\hat{a}^{\dagger}\hat{a}+\sum_{k,j}\left(t_{j,j+1}e^{-i\theta_{j,j+1}}\hat{c}_{j,k}^{\dagger}\hat{c}_{j+1,k}+{\rm
H.c.}\right)~{},$ (3)

where $\omega_{\rm c}$ is the cavity photon energy ($\hbar=1$), and the
Peierls phase
$\theta_{j,j+1}=(-e/c)\int_{{\bm{r}}_{j}}^{{\bm{r}}_{j+1}}{\bm{A}}({\bm{r}})\cdot
d^{2}{\bm{r}}$ is necessary to minimally couple the matter degrees of freedom
living on the plaquette to the cavity field [41].

In the symmetric gauge, the vector potential of the single photon mode is
${\bm{A}}({\bm{r}})=-By{\bm{e}}_{x}/2+Bx{\bm{e}}_{y}/2$, with
$\nabla_{\bm{r}}\times{\bm{A}}({\bm{r}})={\bm{B}}=B{\bm{e}}_{z}$. Here,
${\bm{e}}_{x}$, ${\bm{e}}_{y}$, and ${\bm{e}}_{z}$ are unit vectors in the
$x$, $y$, and $z$ directions, respectively. Quantization of the cavity field
is carried out in the usual manner by promoting $B$ from a $c$-number to a
bosonic quantum operator $B\rightarrow B_{0}(\hat{a}+\hat{a}^{\dagger})$,
where $\hat{a}^{\dagger}$ ($\hat{a}$) creates (destroys) a cavity photon. The
final Hamiltonian describing light-matter interactions in the cavity is:

$\displaystyle\hat{\mathcal{H}}=\omega_{\rm
c}\hat{a}^{\dagger}\hat{a}+\sum_{k,j}\left[t_{j,j+1}e^{-i\lambda(\hat{a}+\hat{a}^{\dagger})/{\sqrt{N}}}\hat{c}_{j,k}^{\dagger}\hat{c}_{j+1,k}+{\rm
H.c.}\right]~{},$ (4)

where $\lambda=-2\pi[\Phi/(\Phi_{0}n_{s})]\sqrt{N}$ is a dimensionless light-
matter coupling constant, which is proportional to the ratio between the
magnetic flux $\Phi\equiv B_{0}A_{\rm p}$ piercing the plaquette of area
$A_{\rm p}$ and the flux quantum $\Phi_{0}=2\pi c/e$. For example, for a
square plaquette ($n_{\rm s}=4$) of side $d$, $A_{\rm p}=d^{2}$. Note that i)
in the thermodynamic $N\to\infty$ limit, $\lambda$ is independent of $N$,
since $B_{0}$ scales as $\sim 1/\sqrt{N}$ to make sure that the magnetic field
$B=B_{0}(\hat{a}+\hat{a}^{\dagger})$ is an intensive quantity (i.e. it does
not scale with $N$) and ii) the physical flux is
$\hat{\Phi}=\Phi(\hat{a}+\hat{a}^{\dagger})$. The Peierls phase introduced in
Eq. (4) is necessary to satisfy the gauge principle in the presence of
magnetic fields and, more in general, in any theory beyond the dipole
approximation [43, 42, 26].

Mean-field theory and photon condensation criterion.—To the end of studying
the possible emergence of photon condensation, we approximate the ground-state
of $\hat{\cal H}$ as a product state of the form
$\ket{\Psi}=\ket{\psi}\ket{\phi}$, where $\ket{\psi}$ and $\ket{\phi}$ are the
matter and light quantum states [44, 45, 46, 47, 48], respectively. In the
thermodynamic limit, the photonic state can be considered a coherent state
$\hat{a}\ket{\phi_{\alpha}}=\alpha\sqrt{N}\ket{\phi_{\alpha}}$. Photon
condensation occurs when the photonic order parameter
$\alpha\equiv\langle\phi_{\alpha}|\hat{a}|\phi_{\alpha}\rangle/\sqrt{N}$
acquires a non-zero value in the thermodynamic limit. From now on, we will
take $\alpha\in{\mathbb{R}}$, without loss of generality. We hasten to
emphasize that our order parameter $\alpha$ is not plagued by any gauge
ambiguity [49] since it physically corresponds to the magnetic flux—defined as
the expectation value of $\hat{\Phi}$ over the ground state
$\ket{\phi_{\alpha}}$, i.e.
$2\alpha\Phi\sqrt{N}=\langle\phi_{\alpha}|\hat{\Phi}|\phi_{\alpha}\rangle$—which
is measurable e.g. via SQUID or NV-center magnetometry [50, 51].

We now consider the following mean-field matter Hamiltonian, describing only
the electronic degrees of freedom:

$\frac{\hat{\mathcal{H}}_{\rm
MF}(\alpha)}{N}\equiv\frac{\matrixelement{\phi_{\alpha}}{\mathcal{\hat{H}}}{\phi_{\alpha}}}{N}=\omega_{\rm
c}\alpha^{2}+\frac{1}{N}\sum^{N}_{k=1}\hat{h}_{{\rm e},k}(\alpha)~{},$ (5)

where

$\hat{h}_{{\rm e},k}(\alpha)=\sum_{j=0}^{n_{\rm
s}-1}\left(t_{j,j+1}e^{-i2\lambda\alpha}\hat{c}_{j,k}^{\dagger}\hat{c}_{j+1,k}+{\rm
H.c.}\right)~{}.$ (6)

Notice that $\hat{h}_{{\rm e},k}(0)$ coincides with the single-molecule
Hamiltonian $\hat{h}_{{\rm e},k}$ in the absence of cavity, which has been
defined above in Eq. (1). The quantity $\hat{\mathcal{H}}_{\rm MF}(\alpha)/N$
in Eq. (5) can be interpreted as an effective mean-field matter Hamiltonian
per molecule.

We denote by the symbol $\ket{\varphi_{l}(\alpha)}_{k}$ and
$\epsilon_{l}(\alpha)$, with $l=0,1,2,\ldots,n_{\rm s}-1$, the eigenstates and
the corresponding eigenvalues of the Hamiltonian $\hat{h}_{{\rm
e},k}(\alpha)$: $\hat{h}_{{\rm
e},k}(\alpha)\ket{\varphi_{l}(\alpha)}_{k}=\epsilon_{l}(\alpha)\ket{\varphi_{l}(\alpha)}_{k}$.
The spectrum of $\hat{h}_{{\rm e},k}(\alpha)$ does not depend on $k$ since all
molecules are identical.

A generic many-body eigenstate of the effective Hamiltonian
$\hat{\mathcal{H}}_{\rm MF}(\alpha)/N$ can be written as [52]
$\ket{\psi_{n}(\alpha)}=\prod^{N}_{k=1}\ket{\varphi_{l_{k,n}}(\alpha)}_{k}$.
Here, $l_{k,n}$ is a discrete index that specifies which single-particle state
is occupied for the $k$-th molecule and $n=0,1,2,\ldots,n^{N}_{\rm s}-1$ is an
integer. In particular, the many-body ground state (i.e. $l_{k,0}=0~{}\forall
k$) is given by
$\ket{\psi_{0}(\alpha)}=\prod^{N}_{k=1}\ket{\varphi_{0}(\alpha)}_{k}$.
Finally, we introduce the energy per molecule
$\bar{\epsilon}_{n}(\alpha)=\sum^{N}_{k=1}\epsilon_{l_{k,n}}(\alpha)/N$. The
ground-state energy is $\bar{\epsilon}_{0}(\alpha)=\epsilon_{0}(\alpha)$. We
emphasize that the quantities $\bar{\epsilon}_{n}(\alpha)$ and
$\epsilon_{l}(\alpha)$ do not include the electromagnetic energy $\omega_{\rm
c}\alpha^{2}$.

At the onset of the phase transition, $\alpha$ is a small parameter and the
mean-field Hamiltonian can be therefore expanded in a power series of
$\alpha$, retaining only terms up to ${\cal O}(\alpha^{2})$. Using the
magnetization operator derived in Ref. [53],

$\displaystyle\hat{M}_{z}(\alpha)=\frac{\lambda}{B_{0}\sqrt{N}}\sum_{k,j}\left(it_{j,j+1}e^{-i2\lambda\alpha}\hat{c}_{j,k}^{\dagger}\hat{c}_{j+1,k}+{\rm
H.c.}\right)~{},$ (7)

we get

$\displaystyle\frac{\hat{\mathcal{H}}_{\rm MF}(\alpha)}{N}=\omega_{\rm
c}\alpha^{2}+\hat{\mathcal{H}}_{0}+\hat{\mathcal{M}}_{\rm
p}(0)\alpha+\frac{1}{2}\hat{\mathcal{M}}_{\rm d}(0)\alpha^{2}~{},$ (8)

where $\hat{\mathcal{H}}_{0}=\hat{\cal H}_{\rm e}/N$, and we have introduced
the paramagnetic $\hat{\mathcal{M}}_{\rm
p}(\alpha)=\sum^{N}_{k=1}\partial_{\alpha}\hat{h}_{{\rm
e},k}(\alpha)/N=-2B_{0}\hat{M}_{z}(\alpha)/\sqrt{N}$ and diamagnetic
$\hat{\mathcal{M}}_{\rm d}(\alpha)=\partial_{\alpha}\hat{\mathcal{M}}_{\rm
p}(\alpha)$ contributions to the magnetic moment operator. We emphasize that
since $B_{0}\sim 1/\sqrt{N}$ and $\hat{M}_{z}(\alpha)\sim N$, both
$\hat{\mathcal{M}}_{\rm p}(\alpha)$ and $\hat{\mathcal{M}}_{\rm d}(\alpha)$
are intensive quantities. In Ref. [53] we show that $\hat{\mathcal{M}}_{\rm
p}(\alpha)$ is proportional to the paramagnetic current operator.

As detailed in Ref. [53], the conditions for photon condensation to take place
can be divided in two classes:

(i) If matter, decoupled from light, displays a non-zero paramagnetic
magnetization, $\matrixelement{\psi_{0}(0)}{\hat{\mathcal{M}}_{\rm
p}(0)}{\psi_{0}(0)}\neq 0$—where $\ket{\psi_{0}(0)}$ is the ground state of
$\hat{\mathcal{H}}_{0}$ with eigenvalue $\bar{\epsilon}_{0}(0)$—the coupled
light-matter system is always in a photon condensate state for any value of
the light-matter coupling $\lambda$, and the real ground state energy shift is
linear in $\alpha$. The physical reason for this phenomenon is that, since
$\matrixelement{\psi_{0}(0)}{\hat{\mathcal{M}}_{\rm p}(0)}{\psi_{0}(0)}\neq
0$, the matter ground state $\ket{\psi_{0}(0)}$ carries a persistent current,
which, in turn, creates a non-zero magnetic field in the cavity. This self-
generated field corresponding to a finite photonic displacement is the
manifestation of photon condensation. Persistent ground-state currents in our
planar molecules can be obtained, for example, by applying an external
classical magnetic field, which yields a finite $\Theta$. From now on, we will
focus only on $\Theta=0$;

(ii) Conversely, if there are no circulating currents in the uncoupled matter
system, i.e. if $\Theta=0$ and
$\matrixelement{\psi_{0}(0)}{\hat{\mathcal{M}}_{\rm p}(0)}{\psi_{0}(0)}=0$,
photon condensation occurs if and only if [53]

$\chi_{\rm M}\geq\omega_{\rm c}~{},$ (9)

where $\chi_{\rm M}=-\matrixelement{\psi_{0}(0)}{\hat{\mathcal{M}}_{\rm
d}(0)}{\psi_{0}(0)}/2-\chi_{\rm p}(0)/2$ is the magnetic susceptibility, and
$\chi_{\rm p}(0)=-2\sum_{n\neq 0}|\langle{\psi_{n}(0)}|{\hat{\mathcal{M}}_{\rm
p}(0)}|{\psi_{0}(0)}\rangle|^{2}/[\bar{\epsilon}_{n}(0)-\bar{\epsilon}_{0}(0)]$
is the static paramagnetic susceptibility.

\begin{overpic}[width=216.81pt]{criteria.pdf}\put(-1.0,92.0){(a)}\put(-1.0,58.0){(b)}\put(20.0,17.0){\rm\tiny{\color[rgb]{1,0,0}Diamagnetism}}\put(57.0,38.0){\rm\tiny{\color[rgb]{1,0,0}Paramagnetism}}\put(76.95,19.5){$\bullet$~{}{\rm
QCP}}\end{overpic} Figure 2: (Color online) Photon condensation of molecules
in a chiral cavity. Results in this figure have been obtained by setting
$\Theta=0$, $\lambda=0.3$, and $t=-\omega_{\rm c}$. (a) Eigenvalues
$\epsilon_{l}(\alpha)$ of the single-molecule Hamiltonian $\hat{h}_{{\rm
e},k}(\alpha)$ in Eq. (6), evaluated at $\alpha=\tilde{\alpha}$ and plotted as
functions of $\tau/\omega_{\rm c}$ (the eigenvalues do not depend on $k$).
Solid lines denote results obtained for $\lambda=0.3$. Dashed lines denote
results in the absence of cavity (i.e. for $\lambda=0$). Note a degeneracy in
the molecular spectra at $\lambda=0$ (dashed lines) and $\tau/\omega_{\rm
c}=1$ (i.e. $\tau=-t$). The vertical dotted line marks the critical value of
$\tau$ at which a transition to a photon condensate state occurs. (b) The
photon condensate order parameter $\tilde{\alpha}$ (solid blue line) is
plotted as a function of $\tau$. When $\tau$ reaches a critical value,
$\tilde{\alpha}$ becomes finite, signalling a quantum phase transition to a
photon condensate state. The dashed red line denotes the orbital magnetic
susceptibility $\chi_{\rm M}$. The transition occurs when $\chi_{\rm
M}/\omega_{\rm c}>1$—see the criterion (9). Diamagnetism (paramagnetism)
corresponds to $\chi_{\rm M}<0$ ($\chi_{\rm M}>0$). Paramagnetism appears in
the Van Vleck regime, when the ground state and the first excited state of the
molecule are close in energy.

Eq. (9) is the most important result of this Letter. Despite it was derived
for a toy-model molecular Hamiltonian, we believe that its range of validity
is much more ample. For example, the inclusion of electron-electron
interactions is not expected to modify Eq. (9) but to dramatically alter the
dependence of $\chi_{\rm M}$ on the microscopic molecular parameters. Also,
the actual details of the cavity will certainly matter but deviations from our
toy-model temporally-chiral cavity hosting a spatially-uniform fluctuating
${\bm{B}}$ field can be easily encoded in the right-hand side of the
inequality (9), changing $\omega_{\rm c}$ into a more complicated
electromagnetic form factor. Our criterion (9) provides guidance in the
experimental search for photon condensation in polaritonic chemistry [55],
emphasizing that the quest for this exotic state of matter should focus on the
combination between molecular systems with a positive orbital magnetic
susceptibility $\chi_{\rm M}$ and chiral cavities.

Variational theory of the photon condensate state.—The perturbative approach
for $\alpha\ll 1$ described so far served only to reach the criterion (9) for
photon condensation. If one is interested in the actual calculation of the
order parameter $\alpha$ and the spectra ${\bar{\epsilon}}_{l}$ of the coupled
light-molecular system as functions of the microscopic parameters of the
model, a non-perturbative approach is needed. To this end, we fix the
parameter $\alpha$ by imposing that the optimal value $\tilde{\alpha}$ yields
a minimum of the mean-field ground-state energy functional ${\cal
E}(\alpha)\equiv\bra{\psi_{0}(\alpha)}\hat{\cal H}_{\rm
MF}(\alpha)\ket{\psi_{0}(\alpha)}/N$. Hence, we determine $\tilde{\alpha}$ by
imposing that $\partial_{\alpha}{\cal E}(\alpha)|_{\alpha=\tilde{\alpha}}=0$:

$\left.\partial_{\alpha}{\cal
E}(\alpha)\right|_{\alpha=\tilde{\alpha}}=2\omega_{\rm
c}\tilde{\alpha}+\matrixelement{\psi_{0}(\tilde{\alpha})}{\hat{\mathcal{M}}_{\rm
p}(\tilde{\alpha})}{\psi_{0}(\tilde{\alpha})}=0~{}.$ (10)

The solution of Eq. (10) not only determines whether the system is in a normal
($\tilde{\alpha}=0$) or photon condensate ($\tilde{\alpha}\neq 0)$ phase but
also yields $\tilde{\alpha}$ as a function of the microscopic parameters
$\tau$, $\Theta$, and $\lambda$. Results for $\tilde{\alpha}$ as a function of
$\tau/\omega_{\rm c}$ and $\Theta$ for $\lambda=0.3$ have been reported in
Fig. 1.

Figure 2(a) shows the molecular spectra $\epsilon_{l}(\tilde{\alpha})$ as
functions of $\tau/\omega_{\rm c}$. The dashed lines describe the molecular
spectra in the absence of cavity (i.e. for $\lambda=0$), while the solid lines
describe the case of a finite light-matter coupling ($\lambda=0.3$). The
vertical dotted line marks the critical value of $\tau$ beyond which a
transition to a photon condensate state occurs. We clearly see that at this
value of $\tau$ the spectra are largely affected by the cavity. Fig. 2(b)
shows the order parameter $\tilde{\alpha}$ (solid blue line) and the orbital
magnetic susceptibility $\chi_{\rm M}$ (red dashed line) as functions of
$\tau/\omega_{\rm c}$. In agreement with Eq. (9), photon condensation (i.e.
$\tilde{\alpha}\neq 0$) occurs when $\chi_{\rm M}/\omega_{\rm c}>1$.
Diamagnetism (paramagnetism) corresponds to $\chi_{\rm M}<0$ ($\chi_{\rm
M}>0$).

Polaritons.—Measuring directly the molecular spectrum
$\epsilon_{l}(\tilde{\alpha})$ in the presence of the cavity is of course
possible but challenging. Light-matter interactions yield also polaritons,
i.e. hybrid light-matter collective modes, which can be measured in a variety
of ways, including scanning probe methods [56, 57, 58, 59]. In order to find
polaritons in our system, we need to study Gaussian fluctuations around the
mean-field state described by the Hamiltonian (5). To this aim, we write the
photon operators as $\hat{a}\rightarrow\tilde{\alpha}\sqrt{N}+\delta\hat{a}$,
where $\delta\hat{a}$ describes a zero-average fluctuation around the mean-
field solution $\tilde{\alpha}\sqrt{N}$, which has been described above. We
then introduce the collective “bright mode” creation operator [60, 61, 62, 63,
64, 65],

$\hat{b}_{l}^{\dagger}\equiv\frac{1}{\sqrt{N}}\sum^{N}_{k=1}\hat{\varphi}^{\dagger}_{l,k}\hat{\varphi}_{0,k}$
(11)

with $l>0$, where $\hat{\varphi}^{\dagger}_{l,k}$ ($\hat{\varphi}_{0,k}$)
creates (destroys) an electron in the eigenstate
$\ket{\varphi_{l}(\tilde{\alpha})}_{k}$
($\ket{\varphi_{0}(\tilde{\alpha})}_{k}$). The bright mode collective operator
creates an electron-hole transition with a finite electrical dipole moment by
annihilating an electron in the ground state of each molecule and creating an
electron in an excited state $l>0$ with energy $\epsilon_{l}(\tilde{\alpha})$.
The associate transition energy is
$\epsilon_{l}(\tilde{\alpha})-\epsilon_{0}(\tilde{\alpha})$. In the
thermodynamic $N\to\infty$, $\hat{b}^{\dagger}_{l}$ behaves as a quasi-bosonic
operator [60, 61, 62, 63, 64, 65], i.e.
$[\hat{b}_{m},\hat{b}^{\dagger}_{l}]\approx\delta_{m,l}$. Further details are
given in Ref. [53].

Hence, in the thermodynamic limit, and expanding the shifted Hamiltonian up to
the second order in the photonic fluctuations $\delta\hat{a}$ and in the
bright mode bosonic operators $b_{l}^{\dagger}$, we can write an approximate
polaritonic Hamiltonian $\hat{\cal H}_{\rm pol}$, describing the lowest
excited states of the coupled cavity-molecular system:

$\displaystyle\hat{\cal H}_{\rm pol}$ $\displaystyle=$
$\displaystyle\omega_{\rm
c}\delta\hat{a}^{\dagger}\delta\hat{a}+\sum_{l=1}^{n_{\rm
s}-1}[\epsilon_{l}(\tilde{\alpha})-\epsilon_{0}(\tilde{\alpha})]\hat{b}_{l}^{\dagger}\hat{b}_{l}$
(12) $\displaystyle+$
$\displaystyle\frac{1}{2}(\delta\hat{a}+\delta\hat{a}^{\dagger})\sum_{l=1}^{n_{\rm
s}-1}\left[\mathcal{M}_{\rm
p}^{l,0}(\tilde{\alpha})\hat{b}_{l}^{\dagger}+\mathcal{M}_{\rm
p}^{0,l}(\tilde{\alpha})\hat{b}_{l}\right]$ $\displaystyle+$
$\displaystyle\frac{1}{8}\mathcal{M}_{\rm
d}^{0,0}(\tilde{\alpha})(\delta\hat{a}+\delta\hat{a}^{\dagger})^{2}~{},$

where
$\mathcal{M}_{\gamma}^{l,m}(\tilde{\alpha})=\sum^{N}_{k=1}\tensor{[}_{k}]{\langle}{}{\varphi_{l}(\tilde{\alpha})}|\hat{\cal
M}_{\gamma}(\tilde{\alpha})\ket{\varphi_{m}(\tilde{\alpha})}_{k}$ with
$\gamma={\rm p},{\rm d}$. Note that $\mathcal{M}_{\rm
p}^{0,0}(\tilde{\alpha})=\langle\psi_{0}(\tilde{\alpha})|\hat{\cal M}_{\rm
p}(\tilde{\alpha})|\psi_{0}(\tilde{\alpha})\rangle$ and $\mathcal{M}_{\rm
d}^{0,0}(\tilde{\alpha})=\langle\psi_{0}(\tilde{\alpha})|\hat{\cal M}_{\rm
d}(\tilde{\alpha})|\psi_{0}(\tilde{\alpha})\rangle$. Finally, the polariton
frequencies can be derived by diagonalizing the Hopfield matrix $\Xi$, i.e.
the matrix that represents the quadratic polaritonic Hamiltonian in Eq. (12).
Further details are reported in Ref. [53].

\begin{overpic}[width=433.62pt]{polaritons_1.pdf} \put(-1.0,88.0){(a)}
\put(-1.0,58.0){(b)} \end{overpic} Figure 3: (Color online) Polariton
softening at finite $\lambda$. Results in this figure have been obtained by
setting with $\Theta=0$, $\tau=0.5\omega_{\rm c}$, and $t=-\omega_{\rm c}$.
(a) The single-molecule ground-state energy $\epsilon_{0}(\tilde{\alpha})$ is
plotted as a function of $\lambda$. (b) The four polariton energies
$\Omega_{\rm p}$ are plotted as functions of $\lambda$. The most important
feature of this panel is the softening of the lowest polariton mode, which
occurs at the quantum phase transition to a photon condensate state.

Figure 3(a) shows the single-molecule ground-state energy
$\epsilon_{0}(\tilde{\alpha})$ for $\tau=0.5\omega_{\rm c}$, as a function of
the normalized collective coupling $\lambda$. Results in this panel have been
obtained by diagonalizing the single-molecule Hamiltonian $\hat{h}_{{\rm
e},k}(\tilde{\alpha})$ in Eq. (6). When $\lambda$ exceeds a critical value
$\lambda_{\rm c}$, $\tilde{\alpha}$ increases from zero to a finite value.
Physically, this means that at $\lambda>\lambda_{\rm c}$ a magnetic flux
$2\tilde{\alpha}\Phi\sqrt{N}$ appears spontaneously and a paramagnetic
molecule responds to it by reducing its energy. This explains why
$\epsilon_{0}(\tilde{\alpha})$ decreases with increasing $\lambda>\lambda_{\rm
c}$, while it does not change with $\lambda$ for $\lambda<\lambda_{\rm c}$
(when there is no self-generated flux). Figure 3(b) shows the four polariton
frequencies $\Omega_{\rm p}$ as functions of $\lambda$ and for a fixed value
of $\tau$. We clearly see that, at $\lambda=\lambda_{\rm c}$, the lowest
polariton mode softens, signalling that the transition to a photon condensate
state is a second-order quantum phase transition.

\begin{overpic}[width=433.62pt]{polaritons_2.pdf} \put(-1.0,88.0){(a)}
\put(-1.0,58.0){(b)} \end{overpic} Figure 4: (Color online) Polariton
softening at $\lambda\to 0$. Results in this figure have been obtained by
setting with $\Theta=0$, $\tau=\omega_{\rm c}$, and $t=-\omega_{\rm c}$. (a)
The single-molecule ground-state energy $\epsilon_{0}(\tilde{\alpha})$ is
plotted as a function of $\lambda$. (b) The four polariton energies
$\Omega_{\rm p}$ are plotted as functions of $\lambda$. Softening of the
lowest polariton mode occurs for $\lambda\to 0$.

The physics discussed so far does not require large values of the light-matter
coupling $\lambda$. Although reaching the ultra-strong coupling regime [66,
55, 67, 68] is possible in a variety of condensed matter and quantum chemistry
setups, it is highly desirable to have toy models and realistic systems where
photon condensation occurs in the weak-coupling $\lambda\to 0$ limit. In Fig.
4 we show that our model displays such pleasant feature. In this figure we
illustrate the same quantities as in Fig. 3 but for a larger value of $\tau$,
i.e. for $\tau=\omega_{\rm c}$. In this case, polariton softening occurs in
the weak-coupling $\lambda\to 0$ limit. The reason is easy to understand. The
criterion (9) for the occurrence of photon condensation in molecular systems
we derived in this Letter depends on the intrinsic orbital magnetic response
$\chi_{\rm M}$ of the molecular system, i.e. on the dependence of $\chi_{\rm
M}$ on the microscopy of the molecular Hamiltonian $\hat{h}_{{\rm e},k}$ in
the absence of light-matter interactions. For photon condensation to occur is
more important to hunt for molecular systems with a large value of $\chi_{\rm
M}>0$ and design cavities with a suitable electromagnetic vacuum structure so
that $\chi_{\rm M}/\omega_{\rm c}>1$ rather than achieving ultra-strong
coupling.

In summary, we have shown that photon condensation can occur also in molecular
systems (and not only in extended electronic systems [27, 28, 29]). The recipe
for achieving it is encoded in the simple and elegant criterion we derived in
Eq. (9). One needs to find molecules with a large and positive orbital
magnetic susceptibility $\chi_{\rm M}$ and place them in cavities hosting a
significant magnetic component of the electromagnetic field. In order to grasp
the essential physics, we have deliberately analyzed, for the sake of
simplicity, single-electron toy-model molecules placed inside a temporally-
chiral cavity with a uniform magnetic field. We hope that our results will
stimulate future work on real molecules loaded into more complicated chiral
cavities [40], which can be studied with recently developed ab initio
numerical approaches [69]. Results similar to ours, which show the importance
of a cavities with a significant magnetic component, have been also found in a
two-leg ladder model: see Ref. [71].

###### Acknowledgements.

Acknowledgements.—This work was supported by: i) the European Union’s Horizon
2020 research and innovation programme under grant agreement no. 881603 -
GrapheneCore3; ii) the University of Pisa under the “PRA - Progetti di Ricerca
di Ateneo” (Institutional Research Grants) - Project No. PRA_2020-2021_92
“Quantum Computing, Technologies and Applications”; and iii) the Italian
Minister of University and Research (MUR) under the “Research projects of
relevant national interest - PRIN 2020” - Project no. 2020JLZ52N, title
“Light-matter interactions and the collective behavior of quantum 2D materials
(q-LIMA)”. This work has been partially funded by European Union
(NextGeneration EU), through the MUR-PNRR project SAMOTHRACE (ECS00000022).
F.M.D.P. was supported by the PNRR MUR project PE0000023-NQSTI and the “Piano
di Incentivi per la Ricerca di Ateneo 2020/2022” of the University of Catania,
through the projects “QUAPHENE” and “Q-ICT”. P.J.H. acknowledges support by
the Gordon and Betty Moore Foundation’s EPiQS Initiative through grant
GBMF9463, the Fundacion Ramon Areces, and the ICFO Distinguished Visiting
Professor program. F.H.L.K. acknowledges financial support from the ERC
TOPONANOP (726001), the Government of Catalonia trough the SGR grant, the
Spanish Ministry of Economy and Competitiveness through the Severo Ochoa
Programme for Centres of Excellence in R&D (Ref. SEV-2015-0522) and Explora
Ciencia (Ref. FIS2017- 91599-EXP), Fundacio Cellex Barcelona, Generalitat de
Catalunya through the CERCA program, the Mineco grant Plan Nacional (Ref.
FIS2016-81044-P), and the Agency for Management of University and Research
Grants (AGAUR) (Ref. 2017-SGR-1656). S.S. acknowledges support by the Army
Research Office (ARO) through grant No. W911NF1910065. It is a great pleasure
to thank Z. Bacciconi, T. Chanda, G. Chiriacò, M. Dalmonte, G. Mazza and M.
Schiró and for useful discussions.

## References

  * [1] K. H. Huang, Statistical Mechanics (John Wiley, New York, 1987).
  * [2] J. M. Raimond, M. Brune, and S. Haroche, Rev. Mod. Phys. 73, 565 (2001).
  * [3] S. Haroche, Rev. Mod. Phys. 85, 1083 (2013).
  * [4] F. J. Garcia-Vidal, C. Ciuti , and T. W. Ebbesen, Science 373, 6551 (2021).
  * [5] D. N. Basov, A. Asenjo-Garcia, P. J. Schuck, X. Zhu, and A. Rubio, Nanoph. 10, 1 ( 2021).
  * [6] F. Schlawin, D. M. Kennes, and M. A. Sentef, Appl. Phys. Rev. 9, 011312 (2022).
  * [7] J. Bloch, A. Cavalleri, V. Galitski, M. Hafezi, and A. Rubio, Nature 606, 41 (2022).
  * [8] I. Carusotto and C. Ciuti, Rev. Mod. Phys. 85, 299 (2013).
  * [9] K. Hepp and E. H. Lieb, Ann. Phys. 76, 360 (1973) and Phys. Rev. A 8, 2517 (1973).
  * [10] Y. K. Wang and F. T. Hioe, Phys. Rev. A 7, 831 (1973).
  * [11] R. H. Dicke, Phys. Rev. 93, 99 (1954).
  * [12] See, however, experimental work on Dicke superradiance in open systems formed by Bose-Einstein condensates coupled to optical cavities: H. Ritsch, P. Domokos, F. Brennecke, and T. Esslinger, Rev. Mod. Phys. 85, 553 (2013).
  * [13] K. Rzażewski, K. Wódkiewicz, and W. Żakowicz, Phys. Rev. Lett. 35, 432 (1975); I. Bialynicki-Birula and K. Rza̧żewski, Phys. Rev. A 19, 301 (1979); K. Gawȩdzki and K. Rza̧żewski, Phys. Rev. A 23, 2134 (1981).
  * [14] K. Rza̧żewski and K. Wódkiewicz, Phys. Rev. Lett. 96, 089301 (2006).
  * [15] P. Nataf and C. Ciuti, Nature Commun. 1, 72 (2010).
  * [16] O. Viehmann, J. von Delft, and F. Marquardt, Phys. Rev. Lett. 107, 113602 (2011).
  * [17] L. Chirolli, M. Polini, V. Giovannetti, and A. H. MacDonald, Phys. Rev. Lett. 109, 267404 (2012).
  * [18] F. M. D. Pellegrino, L. Chirolli, V. Giovannetti, R. Fazio, and M. Polini, Phys. Rev. B 89, 165406 (2014).
  * [19] T. Tufarelli, K. R. McEnery, S. A. Maier, and M. S. Kim, Phys. Rev. A 91, 063840 (2015).
  * [20] V. Rokaj, D. Welakuh, M. Ruggenthaler, and A. Rubio, J. Phys. B: At. Mol. Opt. Phys. 51, 034005 (2018).
  * [21] D. De Bernardis, P. Pilar, T. Jaako, S. De Liberato, and P. Rabl, Phys. Rev. A 98, 053819 (2018).
  * [22] D. De Bernardis, T. Jaako, and P. Rabl, Phys. Rev. A 97, 043820 (2018).
  * [23] G. Mazza and A. Georges, Phys. Rev. Lett. 122, 017401 (2019).
  * [24] C. Schäfer, M. Ruggenthaler, V. Rokaj, and A. Rubio, ACS Photon. 7, 975 (2020).
  * [25] G. M. Andolina, F. M. D. Pellegrino, V. Giovannetti, A. H. MacDonald, and M. Polini, Phys. Rev. B 100, 121109(R) (2019).
  * [26] G. M. Andolina, F. M. D. Pellegrino, A. Mercurio, O. Di Stefano, M. Polini, and S. Savasta, Eur. Phys. J. Plus 137, 1348 (2022).
  * [27] P. Nataf, T. Champel, G. Blatter, and D. M. Basko, Phys. Rev. Lett. 123, 207402 (2019).
  * [28] G. M. Andolina, F. M. D. Pellegrino, V. Giovannetti, A. H. MacDonald, and M. Polini, Phys. Rev. B 102, 125137 (2020).
  * [29] D. Guerci, P. Simon, and C. Mora, Phys. Rev. Lett. 125, 257604 (2020).
  * [30] G. Manzanares, T. Champel, D. M. Basko, and P. Nataf, Phys. Rev. B 105, 245304 (2022).
  * [31] G. Mazza and M. Polini, arXiv.2302.00475.
  * [32] J. W. McClure, Phys. Rev. 104, 666 (1956).
  * [33] G. Vignale, Phys. Rev. Lett. 67, 358 (1991).
  * [34] C. Bruder and Y. Imry, Phys. Rev. Lett. 80, 5782 (1998).
  * [35] A. Principi, M. Polini, G. Vignale, and M. I. Katsnelson, Phys. Rev. Lett. 104, 225503 (2010).
  * [36] G. Gómez-Santos and T. Stauber, Phys. Rev. Lett. 106, 045504 (2011).
  * [37] A. Raoux, M. Morigi, J.-N. Fuchs, F. Piéchon, and G. Montambaux, Phys. Rev. Lett. 112, 026402 (2014)
  * [38] F. Piéchon, A. Raoux, J.-N. Fuchs, and G. Montambaux, Phys. Rev. B 94, 134423 (2016).
  * [39] G. Grosso and G. P. Parravicini, Solid State Physics (Academic Press, London, 2014).
  * [40] H. Hübener, U. De Giovannini, C. Schäfer, J. Andberger, M. Ruggenthaler, J. Faist, and A. Rubio, Nat. Mater. 20, 438 (2021).
  * [41] In this work we have not included the Zeeman coupling to the fluctuating cavity magnetic field. It is well known—see, for example, Ref. [28] and J. Romá-Roche, F. Luis, and D. Zueco, Phys. Rev. Lett. 127, 167201 (2021)—that such coupling yields a positive contribution to the magnetic susceptibility, therefore helping photon condensation. Working with spinless electrons yields a conservative criterion for photon condensation.
  * [42] S. Savasta, O. Di Stefano, A. Settineri, D. Zueco, S. Hughes, and F. Nori, Phys. Rev. A 103, 053703 (2021).
  * [43] U. J. Wiese, Ann. Phys. 525, 777 (2013).
  * [44] C. Emary and T. Brandes, Phys. Rev. Lett. 90, 044101 (2003) and Phys. Rev. E 67, 066203 (2003); N. Lambert, C. Emary, and T. Brandes, Phys. Rev. Lett. 92, 073602 (2004).
  * [45] G. Chiriacò, M. Dalmonte, and T. Chanda, Phys. Rev. B 106, 155113 (2022).
  * [46] G. Passetti, C. J. Eckhardt, M. A. Sentef, and D. M. Kennes, arXiv:2212.03011.
  * [47] C. J. Eckhardt, G. Passetti, M. Othman, C. Karrasch, F. Cavaliere, M. A. Sentef, and D. M. Kennes, Commun. Phys. 5, 122 (2022).
  * [48] I. Amelio, L. Korosec, I. Carusotto, and G. Mazza, Phys. Rev. B 104, 235120 (2021).
  * [49] A. Stokes and A. Nazir, Phys. Rev. Lett. 125, 143603 (2020).
  * [50] F. Casola, T. van der Sar, and A. Yacoby, Nat. Rev. Mater. 3, 17088 (2018).
  * [51] J. F. Barry, J. M. Schloss, E. Bauch, M. J. Turner, C. A. Hart, L. M. Pham, and R. L. Walsworth, Rev. Mod. Phys. 92, 015004 (2020).
  * [52] One should write down states $\ket{\psi_{n}(\alpha)}$ in the form of Slater determinants to make sure that the Pauli exclusion principle is obeyed. However, this is not necessary in our case since, in the present work, we are neglecting Coulomb electron-electron interactions. In ground-state energy calculations, the Pauli principle manifests through the appearance of an exchange contribution only when such interactions are taken into account by carrying out either first-order perturbation theory or e.g. a fully self-consistent Hartree-Fock calculation—see G. F. Giuliani and G. Vignale, Quantum Theory of the Electron Liquid (Cambridge University Press, Cambridge, 2005). An analysis of the role of Coulomb interactions is left for future work.
  * [53] See the Supplemental Material File, which gives additional information on the derivation of the magnetization operator, the effect of an external classical field, the derivation of the condensation criterion, the polaritonic spectrum, and the microscopic theory of Van Vleck paramagnetism. In the Supplemental Material Ref. [54] is cited.
  * [54] J.D. Jackson, Classical Electrodynamics (Wiley, 1998).
  * [55] J. Fregoni, F. J. Garcia-Vidal, and J. Feist, ACS Photon. 9, 1096 (2022).
  * [56] D. N. Basov, M. M. Fogler, and F. J. García de Abajo, Science 354, aag1992 (2016).
  * [57] T. Low, A. Chaves, J. D. Caldwell, A. Kumar, N. X. Fang, P. Avouris, T. F. Heinz, F. Guinea, L. Martin-Moreno, and F. Koppens, Nat. Mater. 16, 182 (2017).
  * [58] D. N. Basov, A. Asenjo-Garcia, P. James Schuck, X. Zhu, and A. Rubio, Nanophotonics 10, 549 (2021).
  * [59] A. Reserbat-Plantey, I. Epstein, I. Torre, A. T. Costa, P. A. D. Gonçalves, N. A. Mortensen, M. Polini, J. C. W. Song, N. M. R. Peres, and F. H. L. Koppens, ACS Photonics 8, 85 (2021).
  * [60] J. J. Hopfield, Phys. Rev. 112, 1555 (1958).
  * [61] V. Savona, Z. Hradil, A. Quattropani, and P. Schwendimann, Phys. Rev. B 49, 8774 (1994).
  * [62] C. Ciuti, G. Bastard, and I. Carusotto, Phys. Rev. B 72, 115303 (2005).
  * [63] D. Gerace and L. C. Andreani, Phys. Rev. B 75, 235325 (2007).
  * [64] I. Carusotto, M. Antezza, F. Bariani, S. De Liberato, and C. Ciuti, Phys. Rev. A 77, 063621 (2008).
  * [65] D. Hagenmüller, S. De Liberato, and C. Ciuti, Phys. Rev. B 81, 235303 (2010).
  * [66] A. F. Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, Nat. Rev. Phys. 1, 19 (2019).
  * [67] F. J Garcia-Vidal, C. Ciuti, and T. W. Ebbesen, Science 373, 6551 (2021).
  * [68] C. Genet, J. Faist, and T. W. Ebbesen, Phys. Today 74, 42 (2021).
  * [69] T. S. Haugland, E. Ronca, E. F. Kjøonstad, A. Rubio, and H. Koch, Phys. Rev. X 10, 041043 (2020).
  * [70] Z. Kurucz and K. Mølmer, Phys. Rev. A 81, 032314 (2010).
  * [71] Z. Bacciconi, G.M. Andolina, T. Chanda, G. Chiriacò, M. Schiró and M. Dalmonte, arXiv:2302.09901.

Supplemental Material for:
“Photon condensation, Van Vleck paramagnetism, and chiral cavities”

Alberto Mercurio, 1 Gian Marcello Andolina,2 Francesco M. D. Pellegrino,3,4,5
Omar Di Stefano,1 Pablo Jarillo-Herrero,6 Claudia Felser,7 Frank H. L.
Koppens2,8, Salvatore Savasta1, and Marco Polini9,10,2

1Dipartimento di Scienze Matematiche e Informatiche, Scienze Fisiche e
Scienze della Terra, Università di Messina, I-98166 Messina, Italy

2ICFO-Institut de Ciències Fotòniques, The Barcelona Institute of Science and
Technology, Av. Carl Friedrich Gauss 3, 08860 Castelldefels (Barcelona), Spain

3Dipartimento di Fisica e Astronomia “Ettore Majorana”, Università di
Catania, Via S. Sofia 64, I-95123 Catania, Italy

4INFN, Sez. Catania, I-95123 Catania, Italy

5CNR-IMM, Via S. Sofia 64, I-95123 Catania, Italy

6Department of Physics, Massachusetts Institute of Technology, Cambridge,
Massachusetts, USA

7Max Planck Institute for Chemical Physics of Solids, Nöthnitzer Str. 40,
Dresden 01187, Germany

8ICREA-Institució Catalana de Recerca i Estudis Avançats, Passeig de Lluís
Companys 23, 08010 Barcelona, Spain

9Dipartimento di Fisica dell’Università di Pisa and INFN, Largo Pontecorvo 3,
I-56127 Pisa, Italy

10Istituto Italiano di Tecnologia, Graphene Labs, Via Morego 30, I-16163
Genova, Italy

In this Supplemental Material we provide additional information on the
derivation of the magnetization and current operators, the effect of an
external classical field, the derivation of the condensation criterion, the
polaritonic spectrum, and the microscopic theory of Van Vleck paramagnetism.

##  Appendix A Derivation of the magnetization and current operators

In this Section we derive the explicit form of the magnetization operator for
electrons hopping in a single “molecule” in a second-quantized fashion. The
position operator $\hat{\bm{r}}_{k}$, associated with the $k$-th molecule, can
be written in terms of electron creation and annihilation operators,
$\hat{c}^{\dagger}_{j,k}$ and $\hat{c}_{j,k}$ for an electron roaming on a
polygon with $n_{\rm s}$ sides:

$\hat{\bm{r}}_{k}=\sum^{n_{\rm
s}-1}_{j=0}{\bm{r}}_{j}\hat{c}^{\dagger}_{j,k}\hat{c}_{j,k}~{}.$ (S1)

Here, ${\bm{r}}_{j}=d_{r}[\cos(\gamma_{j}),\sin(\gamma_{j})]^{\rm T}$ is the
position of the $j$-th site measured from the center of the molecule, $d_{r}$
is the distance of the site from the center and $\gamma_{j}=2\pi j/n_{\rm s}$
is the angle subtended between the position vector of the first site and the
one of the $j$-th site. The magnetic moment of a single molecule is defined in
terms of the position operator as [54]

$\hat{\bm{M}}_{k}=\frac{1}{2}\hat{\bm{r}}_{k}\times\Big{(}-\frac{e}{c}\dot{\hat{\bm{r}}}_{k}\Big{)}~{}.$
(S2)

We use the Heisenberg equation of motion to find an expression for the
velocity operator $\dot{\hat{\bm{r}}}_{k}$,

$\dot{\hat{\bm{r}}}_{k}=i[\hat{\cal H},\hat{\bm{r}}_{k}]~{}.$ (S3)

Replacing this result in the definition of the magnetic moment (Eq. (S2)) we
get

$\hat{\bm{M}}_{k}=\frac{-e}{c}\frac{i}{2}\hat{\bm{r}}_{k}\times[\hat{\cal
H},\hat{\bm{r}}_{k}]~{}.$ (S4)

We now consider the magnetic moment along the ${\bm{e}}_{z}$ direction,
$\hat{M}_{z,k}$, and we use the explicit form of the total Hamiltonian
$\hat{\cal H}$ in Eq. (S4). We find

$\displaystyle\hat{M}_{z,k}$
$\displaystyle=\frac{\lambda}{B_{0}\sqrt{N}}\times$
$\displaystyle\sum_{j=0}^{n_{\rm s}-1}$
$\displaystyle\times\left[it_{j,j+1}e^{-i\frac{\lambda}{\sqrt{N}}(\hat{a}+\hat{a}^{\dagger})}\hat{c}_{j,k}^{\dagger}\hat{c}_{j+1,k}+{\rm
H.c.}\right]~{}.$ (S5)

In the presence of $N$ identical molecules, the total magnetization operator
is expressed as

$\hat{M}_{z}=\sum^{N}_{k=1}\hat{M}_{z,k}~{}.$ (S6)

In the mean-field approach described in the main text, we can trace out the
photonic degrees of freedom by projecting the total magnetization onto a
coherent state $\ket{\phi_{\alpha}}$, i.e. onto a state such that
$\hat{a}\ket{\phi_{\alpha}}=\alpha\sqrt{N}\ket{\phi_{\alpha}}$. In the
thermodynamic limit ($N\to\infty$), this mean-field procedure reduces to the
following formal replacement $\hat{a}\to\alpha\sqrt{N}$. In this limit, the
total magnetization operator becomes

$\hat{M}_{z}(\alpha)=\sum^{N}_{k=1}\hat{M}_{z,k}(\alpha)~{},$ (S7)

where

$\hat{M}_{z,k}(\alpha)=\frac{\lambda}{B_{0}\sqrt{N}}\sum^{n_{\rm
s}-1}_{j=0}\left(it_{j,j+1}e^{-i2\lambda\alpha}\hat{c}^{\dagger}_{j,k}\hat{c}_{j+1,k}+{\rm
H.c.}\right)~{}.$ (S8)

To map this operator into a more familiar object with physical dimensions of a
magnetic moment, it is necessary to introduce an effective mass $m_{\rm eff}$
for our electrons hopping through the lattice sites: $m_{\rm
eff}=1/(2|t|A_{\rm p})$, where $A_{\rm p}$ is the area of a plaquette and $t$
is the hopping parameter that we have introduced in the main text. For
instance, for a square plaquette of side $d$ ($n_{\rm s}=4$), the effective
mass is $m_{\rm eff}=1/(2|t|d^{2})$. Introducing the effective mass into the
definition of the magnetization operator we find

$\displaystyle\hat{M}_{z}(\alpha)=-\frac{e}{2m_{\rm
eff}c}\sum_{k,j}\left(i\frac{t_{j,j+1}}{n_{s}|t|}e^{-i2\lambda\alpha}\hat{c}_{j,k}^{\dagger}\hat{c}_{j+1,k}+{\rm
H.c.}\right)~{},$ (S9)

where the prefactor has the natural form of a magnetic moment $e/(2m_{\rm
eff}c)$.

The total current operator $\hat{J}$ can be derived by starting from a
discretized form of the continuity equation for the density operator. The
local density $\hat{n}_{j,k}\equiv\hat{c}^{\dagger}_{j,k}\hat{c}_{j,k}$ obeys
the Heisenberg equation of motion

$\dot{\hat{n}}_{j,k}=i[\hat{\cal H},{\hat{n}}_{j,k}]~{}.$ (S10)

By this comparing this equation to the following discretized continuity
equation,

$\dot{\hat{n}}_{j,k}=-\frac{1}{d}\big{(}\hat{J}_{j,k}-\hat{J}_{j-1,k}\big{)}~{},$
(S11)

where $\hat{J}_{j,k}$ is the current flowing from site $j$ to site $j+1$ and
$d$ is the lattice spacing, we obtain the following expression for the local
current,

${\hat{J}}_{j,k}=d\big{(}it_{j,j+1}e^{-i\lambda(\hat{a}+\hat{a}^{\dagger})/\sqrt{N}}\hat{c}^{\dagger}_{j+1,k}\hat{c}_{j,k}-{\rm
H.c}\big{)}~{}.$ (S12)

The total current $\hat{J}$ is the sum of all local terms, i.e.
$\hat{J}=\sum_{k,j}{\hat{J}}_{j,k}$, and is given by

${\hat{J}}=d\sum_{k,j}\big{(}it_{j,j+1}e^{-i\lambda(\hat{a}+\hat{a}^{\dagger})/\sqrt{N}}\hat{c}^{\dagger}_{j+1,k}\hat{c}_{j,k}+{\rm
H.c}\big{)}~{}.$ (S13)

As usual, in the mean-field approximation, we replace
$\hat{a}\to\alpha\sqrt{N}$ obtaining a current operator acting only on the
matter degrees of freedom:

${\hat{J}}(\alpha)=d\sum_{k,j}\big{(}it_{j,j+1}e^{-i2\lambda\alpha}\hat{c}^{\dagger}_{j+1,k}\hat{c}_{j,k}+{\rm
H.c}\big{)}~{}.$ (S14)

As expected, the current operator ${\hat{J}}(\alpha)$ and the magnetization
operator $\hat{M}_{z}(\alpha)$ in Eq. (S9) are directly proportional to each
other. Specifically, $\hat{M}_{z}(\alpha)=-({e}/{c}){\hat{J}}(\alpha)[{A_{\rm
p}}/({n_{s}d})]$.

##  Appendix B The effect of an external classic field

As we have seen in the main text, the orbital magnetic response changes sign
when $\tau>0$, i.e. when one of hopping parameters in the single-molecule
Hamiltonian (1) has a different sign with respect to the others (we remind the
reader that $t<0$). We now explain how one can achieve this frustrated
condition on the hoppings by considering the action of an external classical
magnetic field. Starting from the mean-field Hamiltonian [cfr. Eq. (5) in the
main text],

$\displaystyle\frac{\hat{\mathcal{H}}_{\rm MF}(\alpha)}{N}$
$\displaystyle=\frac{\matrixelement{\phi_{\alpha}}{\mathcal{\hat{H}}}{\phi_{\alpha}}}{N}=\omega_{\rm
c}\alpha^{2}+\frac{1}{N}\sum^{N}_{k=1}\hat{h}_{{\rm e},k}(\alpha)~{},$
$\displaystyle\hat{h}_{{\rm e},k}(\alpha)$ $\displaystyle=\sum_{j=0}^{n_{\rm
s}-1}\left(e^{-i2\lambda\alpha}t_{j,j+1}\hat{c}^{\dagger}_{j,k}\hat{c}_{j+1,k}+{\rm
H.c.}\right)~{},$

and considering the special case in which all the hopping parameters are equal
to $t<0$ but for $t_{0,1}=e^{-i\Theta}t$, we see that the matter Hamiltonian
reduces to

$\displaystyle\hat{h}_{{\rm e},k}(\alpha)$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{n_{\rm
s}-1}\left(te^{-i2\lambda\alpha}\hat{c}^{\dagger}_{j,k}\hat{c}_{j+1,k}+{\rm
H.c.}\right)$ $\displaystyle+$
$\displaystyle\left(te^{-i\Theta}e^{-i2\lambda\alpha}\hat{c}^{\dagger}_{0,k}\hat{c}_{1,k}+{\rm
H.c.}\right)~{}.$

By applying the unitary transformation $\hat{c}_{0,k}\to
e^{-i\Theta}\hat{c}_{0,k}$ and $\hat{c}_{j,k}\to e^{-i\Theta(j/n_{\rm
s})}\hat{c}_{j,k}$ for $j\neq 0$, we get

$\displaystyle\hat{h}_{{\rm e},k}(\alpha)$ $\displaystyle=\sum_{j=0}^{n_{\rm
s}-1}\left(te^{-i\Theta/n_{\rm
s}}e^{-i2\lambda\alpha}\hat{c}^{\dagger}_{j,k}\hat{c}_{j+1,k}+{\rm
H.c.}\right)~{}.$ (S15)

This suggests that a classical magnetic field $\bm{B}_{\rm cl}=B_{\rm
cl}\bm{e}_{z}$ with $B_{\rm cl}=-c\Theta/(eA_{\rm p})$ can be used to the
change the sign of one of the hopping parameters, thereby paving the way for
orbital paramagnetism in our toy model and the occurrence of photon
condensation.

##  Appendix C Condensation criteria

In this Section we derive a general criterion to study the onset of photon
condensation. Here, we will consider the energy functional
$E[\alpha,\psi]\equiv\langle\psi|\hat{\mathcal{H}}_{\rm
MF}(\alpha)|\psi\rangle/N$ as a function of the photonic order parameter
$\alpha$ and an arbitrary trial wave-function $|\psi\rangle$. By analyzing the
function ${\cal E}(\alpha)={\rm
min}_{\psi}E[\alpha,\psi]=\langle\psi_{0}(\alpha)|\hat{\mathcal{H}}_{\rm
MF}(\alpha)|\psi_{0}(\alpha)\rangle/N$, the appearance of a finite order
parameter $\tilde{\alpha}$ corresponds to the instability of the normal ground
state, ${\cal E}(\alpha)\leq{\cal E}(0)$. Since we are interested only in the
onset of the phase transition, in this derivation we retain only terms up to
second order in $\alpha$.

First of all, we remind the reader that $\bar{\epsilon}_{n}(0)$ and
$\ket{\psi_{n}(0)}$ are the eigenvalues and eigenvectors of the operator
$\hat{\mathcal{H}}_{0}\equiv{\hat{\mathcal{H}}_{\rm e}}/{N}$ (see Eq. (8) in
the main text).

Now we expand the energy functional $E[\alpha,\psi]$ up to second order in
$\alpha$. Hence, we can consistently approximate the diamagnetic magnetization
with the one in absence of light, $\matrixelement{\psi}{\hat{\mathcal{M}}_{\rm
d}(0)}{\psi}\approx\matrixelement{\psi_{0}(0)}{\hat{\mathcal{M}}_{\rm
d}(0)}{\psi_{0}(0)}=\mathcal{M}^{0,0}_{\rm d}(0)$. Up to the second order in
$\alpha$ the energy functional reads

$E[\alpha,\psi]=\Omega\alpha^{2}+\matrixelement{\psi}{\hat{\mathcal{H}}_{0}}{\psi}+\matrixelement{\psi}{\hat{\mathcal{M}}_{\rm
p}(0)}{\psi}\alpha~{},$ (S16)

where $\Omega\equiv\omega_{\rm c}+\mathcal{M}^{0,0}_{\rm d}(0)/2$. By
minimizing with respect to $\psi$ we get ${\cal E}(\alpha)\equiv{\rm
min}_{\psi}E[\alpha,\psi]$. Hence, we need to calculate

${\cal E}(\alpha)=\Omega\alpha^{2}+{\rm
min}_{\psi}\\{\matrixelement{\psi}{(\hat{\mathcal{H}}_{0}+\hat{\mathcal{M}}_{\rm
p}(0)\alpha)}{\psi}\\}~{},$ (S17)

Using second-order perturbation theory we obtain

$\displaystyle{\rm
min}_{\psi}\\{\matrixelement{\psi}{(\hat{\mathcal{H}}_{0}+\hat{\mathcal{M}}_{\rm
p}(0)\alpha)}{\psi}\\}=$ $\displaystyle=$
$\displaystyle\bar{\epsilon}_{0}(0)+\frac{\chi_{\rm
p}(0)}{2}\alpha^{2}+{\mathcal{M}^{0,0}_{\rm p}(0)}\alpha~{},$

where

$\displaystyle{\mathcal{M}^{0,0}_{\rm p}(0)}$ $\displaystyle\equiv$
$\displaystyle\matrixelement{\psi_{0}}{\hat{\mathcal{M}}_{\rm
p}(0)}{\psi_{0}}~{},$ (S19) $\displaystyle{\chi_{\rm p}}(0)$
$\displaystyle\equiv$ $\displaystyle-2\sum_{n\neq
0}\frac{\absolutevalue{\matrixelement{\psi_{n}}{\hat{\mathcal{M}}_{\rm
p}(0)}{\psi_{0}}}^{2}}{\bar{\epsilon}_{n}(0)-\bar{\epsilon}_{0}(0)}\leq 0~{}.$
(S20)

At this point there are two different cases to distinguish: i) the normal
ground state pararamagnetic moment is zero, ${\mathcal{M}^{0,0}_{\rm p}(0)}=0$
or ii) the normal phase has a finite paramagnetic character
${\mathcal{M}^{0,0}_{\rm p}(0)}\neq 0$. In the case i), Eq. (C) reads

${\cal E}(\alpha)=\bar{\epsilon}_{0}(0)+\big{(}\Omega+\frac{\chi_{\rm
p}(0)}{2}\big{)}\alpha^{2}~{}.$ (S21)

It is possible to define the magnetic susceptibility $\chi_{\rm M}$ as the
concavity of the molecular energy $\bar{\epsilon}_{0}(\alpha)$

$\chi_{\rm
M}\equiv-\frac{1}{2}\partial^{2}_{\alpha}\bar{\epsilon}_{0}(\alpha)\Big{|}_{{\alpha}=0}~{}.$
(S22)

In order to calculate $\chi_{\rm M}$ it is useful to notice that
$\bar{\epsilon}_{0}(\alpha)={\cal E}(\alpha)-\omega_{\rm c}\alpha^{2}~{}$. By
comparing Eq. (S21) with Eq. (S22), we enstablish the a relation between the
magnetic susceptibility $\chi_{\rm M}$ and the following microscopic
quantities

$\displaystyle-\chi_{\rm M}+\omega_{\rm c}$ $\displaystyle=$
$\displaystyle\big{(}\Omega+\frac{\chi_{\rm p}(0)}{2}\big{)}~{},$ (S23)

or equivalently $-\chi_{\rm M}={\mathcal{M}^{0,0}_{\rm d}(0)}/2+\chi_{\rm
p}(0)/{2}$. Hence, the criterion for the occurrence of photon condensation
${\cal E}(\alpha)\leq{\cal E}(0)$ is given by

$\chi_{\rm M}\geq\omega_{\rm c}~{}.$ (S24)

Now, we focus on the case ii), where ${\mathcal{M}^{0,0}_{\rm p}(0)}\neq 0$.
Eq.(S17) can be expressed at the linear order in $\alpha$ as

${\cal E}(\alpha)=\bar{\epsilon}_{0}(0)+{\mathcal{M}^{0,0}_{\rm
p}(0)}\alpha~{}.$ (S25)

In this case, it is always possible to find a $\bar{\alpha}$ such that ${\cal
E}(\bar{\alpha})\leq{\cal E}(0)$ by choosing $\bar{\alpha}$ to have a
different sign with respect to ${\mathcal{M}^{0,0}_{\rm p}(0)}$. In this case,
the system is always unstable and displays photon condensation.
$\mathcal{M}^{0,0}_{\rm p}(0)$ can be different from zero if an additional
classic field is introduced. Fig. 1(b) reproduces both conditions (S24) and
(S25) in the case of an ensemble of squared rings, varying one hopping
parameter and the external classic field.

##  Appendix D Bosonization and polaritons

In order to study the polaritonic properties of our coupled light-matter
system, we start from Eq. (4). First of all, as stated in the main text, we
shift the cavity photon operators, $\hat{a}=\alpha\sqrt{N}+\delta\hat{a}$,
where $\delta\hat{a}$ describes zero-average fluctuations around the mean-
field solution $\alpha\sqrt{N}$. Hence, in the thermodynamic limit, and
expanding the shifted Hamiltonian up to second order in the fluctuations
$\delta\hat{a}$, we find the following Hamiltonian:

$\displaystyle\hat{\mathcal{H}}$ $\displaystyle=$ $\displaystyle\omega_{\rm
c}\left[\delta\hat{a}^{\dagger}\delta\hat{a}+\alpha\sqrt{N}\left(\delta\hat{a}+\delta\hat{a}^{\dagger}\right)+\alpha^{2}N\right]$
$\displaystyle+$ $\displaystyle\sum_{k=1}^{N}\hat{h}_{{\rm
e},k}(\alpha)+\frac{\sqrt{N}}{2}\hat{\mathcal{M}}_{\rm
p}(\alpha)\left(\delta\hat{a}+\delta\hat{a}^{\dagger}\right)$ $\displaystyle+$
$\displaystyle\frac{1}{8}\hat{\mathcal{M}}_{\rm
d}(\alpha)\left(\delta\hat{a}+\delta\hat{a}^{\dagger}\right)^{2}~{}.$

It is now useful to introduce fermionic operators, $\hat{\varphi}_{l,k}$ and
$\hat{\varphi}_{l,k}^{\dagger}$, that destroy and create an electron in the
energy eigenstate $\ket{\varphi_{l}(\alpha)}_{k}$ of the Hamiltonian
$\hat{h}_{{\rm e},k}(\alpha)$, defined in Eq. (6) of the main-text, i.e.
$\hat{h}_{{\rm
e},k}(\alpha)\ket{\varphi_{l}(\alpha)}_{k}=\epsilon_{l}(\alpha)\ket{\varphi_{l}(\alpha)}_{k}$.
Notice that $\hat{\varphi}_{l,k}$ depends on the photonic mean-field $\alpha$
and we dropped that explicit dependence for the ease of readability. Here, the
eigenvalues $\epsilon_{l}(\alpha)$ are $k$-independent, since all molecules
are identical. Expressing $\hat{\mathcal{H}}$ in terms of the operators
$\hat{\varphi}_{l,k}$, we obtain

$\displaystyle\hat{\mathcal{H}}$ $\displaystyle=$ $\displaystyle\omega_{\rm
c}\left[\delta\hat{a}^{\dagger}\delta\hat{a}+\alpha\sqrt{N}\left(\delta\hat{a}+\delta\hat{a}^{\dagger}\right)+\alpha^{2}N\right]$
$\displaystyle+$ $\displaystyle\sum_{l=0}^{n_{\rm
s}-1}\epsilon_{l}(\alpha)\sum_{k=1}^{N}\hat{\varphi}_{l,k}^{\dagger}\hat{\varphi}_{l,k}$
$\displaystyle+$
$\displaystyle\frac{1}{2\sqrt{N}}\left(\delta\hat{a}+\delta\hat{a}^{\dagger}\right)\sum_{l,m=0}^{n_{\rm
s}-1}\mathcal{M}_{\rm
p}^{l,m}(\alpha)\sum_{k=1}^{N}\hat{\varphi}_{l,k}^{\dagger}\hat{\varphi}_{m,k}$
$\displaystyle+$
$\displaystyle\frac{1}{8N}\left(\delta\hat{a}+\delta\hat{a}^{\dagger}\right)^{2}\sum_{l,m=0}^{n_{\rm
s}-1}\mathcal{M}_{\rm
d}^{l,m}(\alpha)\sum_{k=1}^{N}\hat{\varphi}_{l,k}^{\dagger}\hat{\varphi}_{m,k}~{},$

where

$\displaystyle\mathcal{M}_{\rm p}^{l,m}(\alpha)$
$\displaystyle=\sum^{N}_{k=1}\tensor{[}_{k}]{\langle}{}{\varphi_{l}(\alpha)}|{\hat{\cal
M}_{\rm p}(\alpha)}\ket{\varphi_{m}(\alpha)}_{k}$
$\displaystyle=-\frac{2B_{0}}{\sqrt{N}}\sum^{N}_{k=1}\tensor{[}_{k}]{\langle}{}{\varphi_{l}(\alpha)}|{\hat{M}_{z,k}(\alpha)}\ket{\varphi_{m}(\alpha)}_{k}$
(S28)

and

$\displaystyle\mathcal{M}_{\rm d}^{l,m}(\alpha)$
$\displaystyle=\sum^{N}_{k=1}\tensor{[}_{k}]{\langle}{}{\varphi_{l}(\alpha)}|{\hat{\cal
M}_{\rm d}(\alpha)}\ket{\varphi_{m}(\alpha)}_{k}$
$\displaystyle=-\frac{2B_{0}}{\sqrt{N}}\sum^{N}_{k=1}\tensor{[}_{k}]{\langle}{}{\varphi_{l}(\alpha)}|{\partial_{\alpha}\hat{M}_{z,k}(\alpha)}\ket{\varphi_{m}(\alpha)}_{k}~{},$
(S29)

In writing Eq. (D) we used that the matrix elements
$\tensor{[}_{k}]{\langle}{}{\varphi_{l}(\alpha)}|{\hat{M}_{z,k}(\alpha)}\ket{\varphi_{m}(\alpha)}_{k}$
and
$\tensor{[}_{k}]{\langle}{}{\varphi_{l}(\alpha)}|{\partial_{\alpha}\hat{M}_{z,k}(\alpha)}\ket{\varphi_{m}(\alpha)}_{k}$
are independent of the $k$ index. Making use of the collective notation

$\hat{\Sigma}_{l,m}\equiv\sum_{k}\hat{\varphi}_{l,k}^{\dagger}\hat{\varphi}_{m,k}~{},$
(S30)

the Hamiltonian takes a more compact form

$\displaystyle\hat{\mathcal{H}}$ $\displaystyle=$ $\displaystyle\omega_{\rm
c}\left[\delta\hat{a}^{\dagger}\delta\hat{a}+\alpha\sqrt{N}\left(\delta\hat{a}+\delta\hat{a}^{\dagger}\right)+\alpha^{2}N\right]$
$\displaystyle+$ $\displaystyle\sum_{l=0}^{n_{\rm
s}-1}\epsilon_{l}(\alpha)\hat{\Sigma}_{l,l}+\frac{1}{2\sqrt{N}}\left(\hat{a}+\hat{a}^{\dagger}\right)\sum_{l,m=0}^{n_{\rm
s}-1}\left[\mathcal{M}_{\rm p}^{l,m}(\alpha)\hat{\Sigma}_{l,m}\right]$
$\displaystyle+$
$\displaystyle\frac{1}{8N}\left(\hat{a}+\hat{a}^{\dagger}\right)^{2}\sum_{l,m=0}^{n_{\rm
s}-1}\left[\mathcal{M}_{\rm d}^{l,m}(\alpha)\hat{\Sigma}_{l,m}\right]~{}.$

By following Ref. [70], we focus on the symmetric Hilbert subspace, which is
spanned by the occupation number states defined as

$\displaystyle\ket{n_{0},m_{1},\cdots,p_{n_{\rm
s}-1}}=\frac{1}{\sqrt{n!m!\cdots p!}}\sum_{\rm
perm}\ket{\varphi_{0}(\alpha)}_{1}\cdots$ (S32) $\displaystyle\times$
$\displaystyle\ket{\varphi_{0}(\alpha)}_{n}\ket{\varphi_{1}(\alpha)}_{n+1}\cdots\ket{\varphi_{1}(\alpha)}_{n+m}\dots$
$\displaystyle\times$ $\displaystyle\ket{\varphi_{n_{\rm
s}-1}(\alpha)}_{N-p+1}\cdots\ket{\varphi_{n_{\rm s}-1}(\alpha)}_{N}~{},$

which means that $n$ molecules have a single electron that occupies the state
$\ket{\varphi_{0}(\alpha)}$, $m$ molecules have a single electron that
occupies the state $\ket{\varphi_{1}(\alpha)}$, $\ldots$, and $p$ molecules
have a single electron that occupies the state $\ket{\varphi_{n_{\rm
s}-1}(\alpha)}$, such that $N=n+m+\ldots+p$. By applying an occupation number
state on collective operators defined in Eq. (S30), we find the following
properties

$\displaystyle\hat{\Sigma}_{l,l^{\prime}}\ket{n_{l},m_{l^{\prime}},\ldots}$
$\displaystyle=$ $\displaystyle\sqrt{(n+1)m}$ (S33) $\displaystyle\times$
$\displaystyle\ket{(n+1)_{l},m_{l^{\prime}}-1,\ldots},$
$\displaystyle\hat{\Sigma}_{l,l}\ket{n_{l},m_{l^{\prime}},\ldots}$
$\displaystyle=$ $\displaystyle n\ket{n_{l},m_{l^{\prime}},\ldots},$ (S34)

where $l\neq l^{\prime}$. Moreover, the collective operators fulfill the
commutator identity
$[\hat{\Sigma}_{l,l^{\prime}},\hat{\Sigma}_{m,m^{\prime}}]=\delta_{l^{\prime},m}\hat{\Sigma}_{l,m^{\prime}}-\delta_{m^{\prime},l}\hat{\Sigma}_{m,l^{\prime}}$.
Here, we introduce $n_{\rm s}-1$ couple of bosonic creation and annihilation
operators $\hat{b}_{l}^{\dagger}$ and $\hat{b}_{l}$ such that

$\displaystyle\hat{b}_{l}\ket{n_{0},m_{l},\cdots}$ $\displaystyle=$
$\displaystyle\sqrt{m}\ket{(n+1)_{0},(m-1)_{l},\cdots}~{},$ (S35)
$\displaystyle\hat{b}_{l}^{\dagger}\hat{b}_{l}\ket{n_{0},m_{l},\cdots}$
$\displaystyle=$ $\displaystyle m\ket{n_{0},m_{l},\cdots}~{},$ (S36)

where $l=1,\ldots,n_{\rm s}-1$, and the mean-field groundstate
$\ket{\psi_{0}(\alpha)}=\ket{N_{0},0,\cdots,0}=\prod^{N}_{k=1}\ket{\varphi_{0}(\alpha)}_{k}$,
that corresponds to all molecules with a single electron in the state
$\ket{\varphi_{0}(\alpha)}$, acts as the vacuum state for any $\hat{b}_{l}$.
Creation operators $\hat{b}^{\dagger}_{l}$ applied on the vacuum state
$\ket{\psi_{0}(\alpha)}$ describe the collective matter excitations, which
represents the bright modes. By comparing Eqs. (S33)-(S34) with Eqs.
(S35)-(S36), we write the collective operators in terms of the bosonic fields
accordingly to a multilevel Holstein-Primakoff transformation [70],

$\displaystyle\hat{\Sigma}_{0,0}$ $\displaystyle=$ $\displaystyle
N-\sum_{l>0}\hat{b}_{l}^{\dagger}\hat{b}_{l}~{},$ (S37)
$\displaystyle\hat{\Sigma}_{l,0}$ $\displaystyle=$
$\displaystyle\hat{b}_{l}^{\dagger}\sqrt{N-\sum_{l^{\prime}>0}\hat{b}_{l^{\prime}}^{\dagger}\hat{b}_{l^{\prime}}}~{},$
(S38) $\displaystyle\hat{\Sigma}_{l,l^{\prime}}$ $\displaystyle=$
$\displaystyle\hat{b}_{l}^{\dagger}\hat{b}_{l^{\prime}}~{},$ (S39)

where $l,l^{\prime}>0$. In the proximity of the mean-field matter ground state
$\ket{\psi_{0}(\alpha)}$, namely for a small number of collective matter
excitations, one can approximate

$\displaystyle\hat{\Sigma}_{l,0}$ $\displaystyle\approx$
$\displaystyle\sqrt{N}\hat{b}_{l}^{\dagger}~{}.$ (S40)

By using the matter collective bosonic fields $\hat{b}_{l}^{\dagger}$ and
$\hat{b}_{l}$ to rewrite the Hamiltonian in Eq. (D), and we get rid off all
terms beyond the second order in the (light and matter) bosonic fields

$\displaystyle\hat{\mathcal{H}}$ $\displaystyle\simeq$
$\displaystyle\omega_{\rm c}\delta\hat{a}^{\dagger}\delta\hat{a}+\omega_{\rm
c}\sqrt{N}\left[\alpha+\frac{\mathcal{M}_{\rm p}^{0,0}(\alpha)}{2\omega_{\rm
c}}\right]\left(\delta\hat{a}+\delta\hat{a}^{\dagger}\right)$ (S41)
$\displaystyle+$ $\displaystyle N\omega_{\rm
c}\alpha^{2}+N\epsilon_{0}(\alpha)+\sum_{l=1}^{n_{\rm
s}-1}\left[\epsilon_{l}(\alpha)-\epsilon_{0}(\alpha)\right]\hat{b}_{l}^{\dagger}\hat{b}_{l}$
$\displaystyle+$
$\displaystyle\frac{1}{2}\left(\delta\hat{a}+\delta\hat{a}^{\dagger}\right)\sum_{l=1}^{n_{\rm
s}-1}\left[\mathcal{M}_{\rm
p}^{l,0}(\alpha)\hat{b}_{l}^{\dagger}+\mathcal{M}_{\rm
p}^{0,l}(\alpha)\hat{b}_{l}\right]$ $\displaystyle+$
$\displaystyle\frac{1}{8}\left(\delta\hat{a}+\delta\hat{a}^{\dagger}\right)^{2}\mathcal{M}_{\rm
d}^{0,0}(\alpha)~{},$

This Hamiltonian still depends upon the mean-field displacement $\alpha$. In
order to fix the value of $\alpha$, we select the value which nullifies the
derivative of the energy functional ${\cal
E}(\alpha)=\bra{\psi_{0}(\alpha)}\hat{\cal H}_{\rm
MF}(\alpha)\ket{\psi_{0}(\alpha)}/N$ with respect $\alpha$, i.e.

$\frac{d{\cal E}(\alpha)}{d\alpha}=2\omega_{\rm
c}\alpha+\frac{d\bar{\epsilon}_{0}(\alpha)}{d\alpha}~{}.$ (S42)

By using the Hellmann–Feynman theorem to compute
$({d\bar{\epsilon}_{0}(\alpha)}/{d\alpha})$,

$\displaystyle\frac{d\bar{\epsilon}_{0}(\alpha)}{d\alpha}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum^{N}_{k=1}\bra{\psi_{0}(\alpha)}\frac{d\hat{h}_{{\rm
e},k}(\alpha)}{d\alpha}\ket{\psi_{0}(\alpha)}~{},$ $\displaystyle=$
$\displaystyle\bra{\psi_{0}(\alpha)}\hat{\cal M}_{\rm
p}(\alpha)\ket{\psi_{0}(\alpha)}=\mathcal{M}_{\rm p}^{0,0}(\alpha),$

we can write a non-linear equation to determine $\tilde{\alpha}$:

$\frac{d{\cal
E}(\alpha)}{d\alpha}\bigg{|}_{\alpha=\tilde{\alpha}}=2\omega_{\rm
c}\tilde{\alpha}+\mathcal{M}_{\rm p}^{0,0}(\tilde{\alpha})=0~{}.$ (S44)

By imposing $\alpha=\tilde{\alpha}$, the linear terms in the bosonic operators
in Eq. (S41) vanish. We clearly see that Eq. (S44) coincides with Eq. (10) in
the main text.

Once we set $\alpha=\tilde{\alpha}$, the approximate Hamiltonian in Eq. (S41)
becomes a polaritonic Hamiltonian $\hat{\cal H}_{\rm pol}$ the low-energy
hybrid excitations (polaritons) on top of the mean-field ground-state
solution,

$\displaystyle\hat{\cal H}_{\rm pol}$ $\displaystyle=$
$\displaystyle\omega_{\rm
c}\delta\hat{a}^{\dagger}\delta\hat{a}+\sum_{l=1}^{n_{\rm
s}-1}[\epsilon_{l}(\tilde{\alpha})-\epsilon_{0}(\tilde{\alpha})]\hat{b}_{l}^{\dagger}\hat{b}_{l}$
(S45) $\displaystyle+$
$\displaystyle\frac{1}{2}(\delta\hat{a}+\delta\hat{a}^{\dagger})\sum_{l=1}^{n_{\rm
s}-1}\left[\mathcal{M}_{\rm
p}^{l,0}(\tilde{\alpha})\hat{b}_{l}^{\dagger}+\mathcal{M}_{\rm
p}^{0,l}(\tilde{\alpha})\hat{b}_{l}\right]$ $\displaystyle+$
$\displaystyle\frac{1}{8}\mathcal{M}_{\rm
d}^{0,0}(\tilde{\alpha})(\delta\hat{a}+\delta\hat{a}^{\dagger})^{2}~{}.$

Polaritons are linear combinations of light and matter operators,

$\hat{p}_{\nu}=X_{\nu}\delta\hat{a}+Y_{\nu}\delta\hat{a}^{\dagger}+\sum^{n_{\rm
s}-1}_{l=1}\big{(}W_{\nu,l}\hat{b}_{l}+Z_{\nu,l}\hat{b}^{\dagger}_{l}\big{)}~{}.$
(S46)

Being the polariton a proper bosonic excitation of the system, the operator
$\hat{p}_{\nu}$ fulfills the equation of motion of a harmonic ladder operator

$[\hat{\cal H}_{\rm pol},\hat{p}_{\nu}]=-\Omega_{{\rm
p},\nu}\hat{p}_{\nu}~{}.$ (S47)

Being the polariton a combination of
$\hat{a},\hat{a}^{\dagger},\hat{b}_{l},\hat{b}_{l}^{\dagger}$, in order to
calculate Eq. (S47) we need the equations of motion for these light and matter
operators. These are given by the following commutators

$\displaystyle[\hat{\cal H}_{\rm pol},\delta\hat{a}]$ $\displaystyle=$
$\displaystyle-\omega_{\rm c}\delta\hat{a}-\frac{1}{4}\mathcal{M}_{\rm
d}^{0,0}(\tilde{\alpha})(\delta\hat{a}+\delta\hat{a}^{\dagger})$ (S48)
$\displaystyle-$ $\displaystyle\frac{1}{2}\sum^{n_{\rm
s}-1}_{l=1}(\mathcal{M}_{\rm
p}^{l,0}(\tilde{\alpha})\hat{b}_{l}^{\dagger}+\mathcal{M}_{\rm
p}^{0,l}(\tilde{\alpha})\hat{b}_{l})~{},$ $\displaystyle[\hat{\cal H}_{\rm
pol},\delta\hat{a}^{\dagger}]$ $\displaystyle=$ $\displaystyle\omega_{\rm
c}\delta\hat{a}^{\dagger}+\frac{1}{4}\mathcal{M}_{\rm
d}^{0,0}(\tilde{\alpha})(\delta\hat{a}+\delta\hat{a}^{\dagger})$ (S49)
$\displaystyle+$ $\displaystyle\frac{1}{2}\sum^{n_{\rm
s}-1}_{l=1}(\mathcal{M}_{\rm
p}^{l,0}(\tilde{\alpha})\hat{b}_{l}^{\dagger}+\mathcal{M}_{\rm
p}^{0,l}(\tilde{\alpha})\hat{b}_{l})~{},$ $[\hat{\cal H}_{\rm
pol},\hat{b}_{l}]=-[\epsilon_{l}(\tilde{\alpha})-\epsilon_{0}(\tilde{\alpha})]\hat{b}_{l}-\frac{1}{2}\mathcal{M}_{\rm
p}^{0,l}(\tilde{\alpha})(\delta\hat{a}+\delta\hat{a}^{\dagger})~{},$ (S50)
$[\hat{\cal H}_{\rm
pol},\hat{b}_{l}^{\dagger}]=\\-[\epsilon_{l}(\tilde{\alpha})-\epsilon_{0}(\tilde{\alpha})]\hat{b}_{l}^{\dagger}+\frac{1}{2}\mathcal{M}_{\rm
p}^{l,0}(\tilde{\alpha})(\delta\hat{a}+\delta\hat{a}^{\dagger})~{},$ (S51)

By using these commutators and Eq. (S46), the equation of motion expressed in
Eq. (S47) can be mapped in the following eigenvalue problem

$\Xi{\bm{v}}_{\nu}=\Omega_{\rm p,\nu}{\bm{v}}_{\nu}~{},$ (S52)

where ${\bm{v}}_{\nu}=(X_{\nu},Y_{\nu},{\bm{W}}_{\nu},{\bm{Z}}_{\nu})^{\rm
T}$, and

$\Xi=\begin{pmatrix}\omega_{\rm c}+\frac{1}{4}\mathcal{M}_{\rm
d}^{0,0}(\tilde{\alpha})&-\frac{1}{4}\mathcal{M}_{\rm
d}^{0,0}(\tilde{\alpha})&{\bm{g}}(\tilde{\alpha})&-{\bm{g}}^{*}(\tilde{\alpha})\\\
\frac{1}{4}\mathcal{M}_{\rm d}^{0,0}(\tilde{\alpha})&-\omega_{\rm
c}-\frac{1}{4}\mathcal{M}_{\rm
d}^{0,0}(\tilde{\alpha})&{\bm{g}}(\tilde{\alpha})&-{\bm{g}}^{*}(\tilde{\alpha})\\\
{{\bm{g}}^{*}}^{\rm T}(\tilde{\alpha})&-{{\bm{g}}^{*}}^{\rm
T}(\tilde{\alpha})&{\bm{\Omega}}(\tilde{\alpha})&0\\\ {\bm{g}}^{\rm
T}(\tilde{\alpha})&-{\bm{g}}^{\rm
T}(\tilde{\alpha})&0&-{\bm{\Omega}}(\tilde{\alpha})\\\ \end{pmatrix}~{},$
(S53)

where ${\bm{g}}(\tilde{\alpha})=({1}/{2})[\mathcal{M}_{\rm
p}^{1,0}(\tilde{\alpha}),\dots,\mathcal{M}_{\rm p}^{n_{\rm
s}-1,0}(\tilde{\alpha})]$ and ${\bm{\Omega}}(\tilde{\alpha})={\rm
diag}[\epsilon_{1}(\tilde{\alpha})-\epsilon_{0}(\tilde{\alpha}),\dots,\epsilon_{n_{\rm
s}-1}(\tilde{\alpha})-\epsilon_{0}(\tilde{\alpha})]$.

Eigenvalues of the Hopfield matrix $\Xi$ can be determined by the roots of the
determinant ${\cal D}(\Omega_{{\rm p},\nu})$

${\cal D}(\Omega_{{\rm p},\nu})={\rm Det}(\Omega_{{\rm
p},\nu}\openone-\Xi)=0~{}.$ (S54)

The determinant ${\cal D}(\Omega_{{\rm p},\nu})$ can be calculated by using
the following algebraic property of the block matrices

${\rm Det}\begin{pmatrix}A&B\\\ C&D\end{pmatrix}={\rm Det}(D){\rm
Det}(A-BD^{-1}C)~{},$

which leads to

$\displaystyle{\cal D}(\Omega_{{\rm p},\nu})$ $\displaystyle=$
$\displaystyle\prod^{n_{\rm s}-1}_{l=1}\\{\Omega_{{\rm
p},\nu}^{2}-[\epsilon_{l}(\tilde{\alpha})-\epsilon_{0}(\tilde{\alpha})]^{2}\\}$
$\displaystyle\times$ $\displaystyle{\rm Det}\begin{pmatrix}\omega_{\rm
c}+[\frac{1}{4}\mathcal{M}_{\rm
d}^{0,0}(\tilde{\alpha})+\frac{1}{4}\chi^{(\tilde{\alpha})}_{\rm
p}(\Omega_{{\rm p},\nu})]-\Omega_{{\rm p},\nu}&-[\frac{1}{4}\mathcal{M}_{\rm
d}^{0,0}(\tilde{\alpha})+\frac{1}{4}\chi^{(\tilde{\alpha})}_{\rm
p}(\Omega_{{\rm p},\nu})]\\\ [\frac{1}{4}\mathcal{M}_{\rm
d}^{0,0}(\tilde{\alpha})+\frac{1}{4}\chi^{(\tilde{\alpha})}_{\rm
p}(\Omega_{{\rm p},\nu})]&-\omega_{\rm c}-[\frac{1}{4}\mathcal{M}_{\rm
d}^{0,0}(\tilde{\alpha})+\frac{1}{4}\chi^{(\tilde{\alpha})}_{\rm
p}(\Omega_{{\rm p},\nu})]-\Omega_{{\rm p},\nu}\end{pmatrix}~{},$

where

$\displaystyle\chi^{(\tilde{\alpha})}_{\rm p}(\omega)$ $\displaystyle=$
$\displaystyle 2\sum^{n_{\rm s}-1}_{l=1}\frac{\absolutevalue{\mathcal{M}_{\rm
p}^{l,0}(\tilde{\alpha})}^{2}[\epsilon_{l}(\tilde{\alpha})-\epsilon_{0}(\tilde{\alpha})]}{\omega^{2}-[\epsilon_{l}(\tilde{\alpha})-\epsilon_{0}(\tilde{\alpha})]^{2}}~{}.$
(S56)

The static limit of $\chi^{(\tilde{\alpha})}_{\rm p}(\omega)$ generalizes the
static paramagnetic susceptibility for a finite value of $\tilde{\alpha}$,
such that $\chi^{(0)}_{\rm p}(0)=\chi_{\rm p}(0)$, where $\chi_{\rm p}(0)$ has
been defined in the main text as

$\chi_{\rm p}(0)=-2\sum_{n\neq
0}\frac{\absolutevalue{\matrixelement{\psi_{n}(0)}{\hat{\mathcal{M}}_{\rm
p}(0)}{\psi_{0}(0)}}^{2}}{\bar{\epsilon}_{n}(0)-\bar{\epsilon}_{0}(0)}~{}.$
(S57)

To verify this generalization, we note that the matrix element
$\matrixelement{\psi_{n}(0)}{\hat{\mathcal{M}}_{\rm p}(0)}{\psi_{0}(0)}$
connects the mean-field many-body ground state $\ket{\psi_{0}(0)}$ with a
mean-field many-body state $\ket{\psi_{n}(0)}$ with only an excited molecule,
namely the $N$-tuple $\\{l_{1,n},\ldots,l_{k,n},\ldots,l_{N,n}\\}$ associated
with the many-body excited state $\ket{\psi_{n}(0)}$ consists of $N-1$ zeroes
and a single non-zero element $l_{k^{\star},n}=l^{\star}_{n}>0$.

Thus, the energy difference between the many-body excited state and the many-
body ground state can be expressed in terms of the single-particle energy
difference between excited states and the ground state as
$\bar{\epsilon}_{n}(0)-\bar{\epsilon}_{0}(0)=[\epsilon_{l^{\star}_{n}}(0)-\epsilon_{0}(0)]/N$.
Moreover, since we are dealing with $N$ identical molecules, there are $N$
different excited states (with a single excitation) associated with the gap
$\bar{\epsilon}_{n}(0)-\bar{\epsilon}_{0}(0)$. Thanks to these remarks and by
means of Eq. (D), the static paramagnetic susceptibility (in Eq. (S57)) can be
expressed as

$\displaystyle\chi_{\rm p}(0)=-2\sum^{n_{\rm
s}-1}_{l=1}\frac{\absolutevalue{\mathcal{M}_{\rm
p}^{l,0}(0)}^{2}}{[\epsilon_{l}({0})-\epsilon_{0}(0)]}~{}.$ (S58)

Comparing this expression with Eq. (S56) evaluated at $\omega=0$, we conclude
that

$\displaystyle\chi_{\rm p}^{(0)}(0)=\chi_{\rm p}(0)~{}.$ (S59)

With few algebraic manipulations, the determinant can be expressed as

$\displaystyle{\cal D}(\Omega_{{\rm p},\nu})$ $\displaystyle=$
$\displaystyle\prod^{n_{\rm s}-1}_{l=1}\\{\Omega_{{\rm
p},\nu}^{2}-[\epsilon_{l}(\tilde{\alpha})-\epsilon_{0}(\tilde{\alpha})]^{2}\\}\bigg{\\{}\Omega_{{\rm
p},\nu}^{2}+$ $\displaystyle-$ $\displaystyle\omega_{\rm c}\left[\omega_{\rm
c}+\left(\frac{1}{2}\mathcal{M}_{\rm
d}^{0,0}(\tilde{\alpha})+\frac{1}{2}\chi^{(\tilde{\alpha})}_{\rm
p}(\Omega_{{\rm p},\nu})\right)\right]\bigg{\\}}~{}.$

Hence, polaritonic eigenenergies are given by the following non-linear
equation

$\Omega_{{\rm p},\nu}^{2}-\omega_{\rm c}\left\\{\omega_{\rm
c}+\left[\frac{1}{2}\mathcal{M}_{\rm
d}^{0,0}(\tilde{\alpha})+\frac{1}{2}\chi^{(\tilde{\alpha})}_{\rm
p}(\Omega_{{\rm p},\nu})\right]\right\\}=0~{},$ (S61)

where we remind that $\chi^{(\tilde{\alpha})}_{\rm p}(\Omega_{{\rm p},\nu})$
is a function of $\tilde{\alpha}$. The onset of the superradiant phase
transition corresponds to a softening of a polariton, i.e. a polariton with
zero energy. Eq. (S61) for $\Omega_{{\rm p},\nu}=0$ reduces to

$\omega_{\rm c}+\frac{1}{2}\mathcal{M}_{\rm
d}^{0,0}(\tilde{\alpha})+\frac{1}{2}\chi^{(\tilde{\alpha})}_{\rm p}(0)=0~{}.$
(S62)

We note that for $\tilde{\alpha}=0$, employing Eq. (S59), the previous
equation coincides with the instability criterion shown in Eq. (S24).

##  Appendix E Microscopic theory of Van Vleck paramagnetism

In this Section we briefly remind the reader about the microscopic theory of
Van Vleck paramagnetism.

We consider a generic many-electron Hamiltonian of the form

$\hat{\cal{H}}_{0}=\sum_{k}\Big{[}\frac{\hat{\bm{p}}^{2}_{k}}{2m}+V(\hat{\bm{r}}_{k})\Big{]}+\frac{1}{2}\sum_{k\neq
k^{\prime}}v(\hat{\bm{r}}_{k}-\hat{\bm{r}}_{k}^{\prime})~{},$ (S63)

where $m$ is the electronic mass, $\hat{\bm{p}}_{k}$ is the momentum of the
$k-$ electron, $V(r)$ is a external potential and
$v(\hat{\bm{r}}_{k}-\hat{\bm{r}}_{k}^{\prime})$ is the electron-electron
interaction. The coupling to a uniform magnetic field ${\bm{B}}$ is made as
customary via the minimal coupling substitution, i.e.
$\hat{\bm{p}}_{k}\to\hat{\bm{p}}_{k}+(e/c)\bm{A}(\bm{r}_{k})$.

In the symmetric gauge, the vector potential can be expressed as follows

$\bm{A}(\bm{r})=\frac{1}{2}B\hat{\bm{e}}_{z}\times\bm{r},$ (S64)

where $\bm{B}=B{\bm{e}}_{z}$ is the magnetic field. Notice that the symmetric-
gauge vector potential obeys also the Coulomb condition, i.e.
$\bm{\nabla}\cdot\bm{A}(\bm{r})=0$. Carrying out the minimal coupling in the
symmetric gauge, we find the following Hamiltonian of the many-electron
systems coupled to a uniform external magnetic field $\bm{B}$:

$\hat{\cal H}=\hat{\cal
H}_{0}+\frac{e}{2mc}\bm{B}\cdot\hat{\bm{L}}+\frac{e^{2}}{8mc^{2}}B^{2}\sum_{k=1}^{N}|\hat{\bm{r}}_{k}|^{2}~{},$
(S65)

where $\hat{\bm{L}}=\sum_{k=1}^{N}(\hat{\bm{r}}_{k}\times\hat{\bm{p}}_{k})$ is
the (paramagnetic) angular momentum operator and the last term is the
magnetostatic energy. Note that, accordingly to the main text, we have
neglected the Zeeman coupling.

The magnetization $\hat{\bm{M}}$ is related to the angular momentum and
magnetic field by the following relation:

$\displaystyle\hat{\bm{M}}$
$\displaystyle\equiv-\frac{ei}{2c}\sum_{k=1}^{N}\hat{\bm{r}}_{k}\times[\mathcal{H},\hat{\bm{r}}_{k}]$
$\displaystyle=-\frac{e}{2mc}\hat{\bm{L}}-\frac{e^{2}}{2mc^{2}}\sum_{k=1}^{N}\hat{\bm{r}}_{k}\times\bm{A}(\bm{r}_{k})$
$\displaystyle=-\frac{e}{2mc}\hat{\bm{L}}-\frac{e^{2}B}{4mc^{2}}\Big{(}\sum_{k=1}^{N}|\hat{\bm{r}}_{k}|^{2}\Big{)}\hat{\bm{z}}$
$\displaystyle\equiv\hat{\bm{M}}_{\rm p}+\hat{\bm{M}}_{\rm d}~{},$ (S66)

where the first (second) term is the paramagnetic (diamagnetic) contribution.

In order to find the orbital magnetic susceptibility, we need to study the
energy variation under the applied magnetic field. This splits in two terms:
$\Delta E=\Delta E_{\rm P}+\Delta E_{\rm D}$. Introducing the exact
eigenstates and eigenvalues of the Hamiltonian $\hat{\cal H}_{0}$, we find:

$\displaystyle\Delta E_{\rm D}$
$\displaystyle\equiv\frac{e^{2}}{8mc^{2}}B^{2}\langle\psi_{0}|\sum_{k=1}^{N}|\hat{\bm{r}}_{k}|^{2}|\psi_{0}\rangle$
$\displaystyle=-\frac{1}{2}\langle\psi_{0}|\hat{{M}}_{{\rm
d},z}|\psi_{0}\rangle B^{2}~{},$ (S67)

where $\hat{{M}}_{{\rm d},z}$ is the $\bm{e}_{z}$ component of the vector
$\hat{\bm{M}}_{\rm d}$ and

$\displaystyle\Delta E_{\rm P}$
$\displaystyle\equiv\Big{(}\frac{e}{2mc}\Big{)}^{2}\sum_{n\neq
0}\frac{|\langle\psi_{n}|\bm{B}\cdot\hat{\bm{L}}|\psi_{0}\rangle|^{2}}{E_{n}-E_{0}}$
$\displaystyle=\sum_{n\neq
0}\frac{|\langle\psi_{n}|\bm{B}\cdot\hat{\bm{M}}_{\rm
p}|\psi_{0}\rangle|^{2}}{E_{n}-E_{0}}~{}.$ (S68)

$\Delta E_{\rm B}$ is the magnetic contribution, while by construction $\Delta
E_{\rm D}>0$ and $\Delta E_{\rm P}<0$ i.e. they are a diamagnetic and a
paramagnetic contribution, respectively. If the first excited state $n=1$ is
nearly degenerate with the ground state $E_{1}-E_{0}\approx 0$, the
paramagnetic contribution is dominant and the system has a paramagnetic
response (known as Van Vleck paramagnetism).

The paramagnetic contribution can be recast as $\Delta E_{\rm
P}=-(1/2){\chi_{\rm p}}B^{2}$ where $\chi_{\rm p}$ is the magnetization
response function,

$\displaystyle\chi_{{\rm p}}$ $\displaystyle\equiv$
$\displaystyle-2\sum_{n\neq
0}\frac{|\langle\psi_{n}|\hat{\bm{z}}\cdot\hat{\bm{M}}_{\rm
p}|\psi_{0}\rangle|^{2}}{E_{n}-E_{0}},$ (S69) $\displaystyle=$
$\displaystyle-2\Big{(}\frac{e}{2mc}\Big{)}^{2}\sum_{n\neq
0}\frac{|\langle\psi_{n}|\hat{\bm{z}}\cdot\hat{\bm{L}}|\psi_{0}\rangle|^{2}}{E_{n}-E_{0}},$
(S70)

Thus, orbital paramagnetism is governed by the angular momentum-angular
momentum response function. From the expression for the response function in
Eq. (S70), it can be seen that non-zero orbital paramagnetic response occurs
only if the system does not have rotational invariance around the
$\hat{\bm{z}}$ axis. Indeed, if $[\hat{\cal
H}_{0},\hat{\bm{z}}\cdot\hat{\bm{L}}]=0$, one can choose a common eigenstate
basis of the energy and the angular momentum projection
$\hat{\bm{z}}\cdot{\bm{L}}$ and the matrix element
$\bra{\psi_{0}}\hat{\bm{z}}\cdot\hat{\bm{L}}\ket{\psi_{n}}$ vanishes,
precluding the presence of Van Vleck paramagnetism in rotationally invariant
systems, such as closed shell atoms.