Photon condensation, Van Vleck paramagnetism, and chiral cavities

We introduce a gauge-invariant model of planar, square molecules coupled to a quantized spatially-varying cavity electromagnetic vector potential A(r). Specifically, we choose a temporally chiral cavity hosting a uniform magnetic field B, as this is the simplest instance in which a transverse spatially-varying A(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.

On a seemingly disconnected path, the impact of vacuum cavity fields on the chemical and physical properties of molecules has been demonstrated [37][38][39][40].A great deal of interest emerged after seminal experiments unveiled how photochemical reaction rates can be modified within cavities [41,42].This blooming field is nowadays known as polaritonic chemistry [40].Fundamental theoretical work in this field [43][44][45][46][47][48] is carried out in the framework of the electrical dipole approximation, whereby molecular transitions couple solely with the cavity electric field, thereby neglecting magnetic effects.The need to transcend the electrical dipole approximation in polaritonic chemistry is beginning to emerge.For example, the authors of Ref. [49] go beyond it in order to distinguish arXiv:2302.09964v3[cond-mat.mes-hall]22 Mar 2024 between two enantiomers in chiral cavities.
The key question we try and answer in this Letter is the following: can photon condensation-which, so far, has been studied only in itinerant electron systems [25][26][27]-occur in the realm of polaritonic chemistry?In this Letter, we show that incorporating the effects of magnetic coupling, photon condensation can occur in a system of molecules, which does not display extended Bloch states and an itinerant character.Two conditions need to be 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" [50].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 A(r) is at play.Since a spatially-varying A leads to a finite magnetic field B(r) = ∇ r ×A(r), the simplest cavity one can consider is a cavity hosting a uniform B field.Our cavity falls into the category of chiral cavities [51] with chirality of temporal (rather than structural) character.

II. MODEL
We consider a single, spinless electron [52] hopping between the n s sites of a planar plaquette, such as a square or triangular plaquette, lying in the x-ŷ plane.(Here,x, ŷ, and ẑ are unit vectors in the x, y, and z directions, respectively.)This system will be considered as a toy model of a "molecule" since the plaquette contains a small number of sites (n s = 4 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, Ĥe = N k=1 ĥe,k , where is the single-molecule Hamiltonian.Here, ĉ † j,k (ĉ j,k ) creates (destroys) an electron on the j-th site of the k-th plaquette, t j,j+1 = δ j,0 t ′ + (1 − δ j,0 )t is the tunnelling amplitude between neighboring sites of a plaquette, and where τ ∈ [−t, t] (with t>0 is assumed to be real, without any loss of generality) and Θ ∈ [0, π[ are the tunneling "intensity" and phase, respectively.We emphasize that we have taken t ′ ̸ = t to simulate different physical conditions.For Θ ̸ = 0 (Θ ∈]0, π[, yielding therefore Im[t ′ ] ̸ = 0), the electronic Hamiltonian is not invariant under time-reversal symmetry (TRS), and the plaquette supports a circulating ground-state current.Viceversa, when Θ = 0 the system is invariant under TRS and no ground-state currents flow.As discussed in Appendix C, Hamiltonians like the one in Eq. ( 1), which break TRS, can be physically obtained by applying an external, classical magnetic field along the ẑ direction.Note that, counterintuitively, also the particular case Θ = 0 and t ′ = −t (which respects TRS) can be obtained by applying a particular classical transverse static magnetic field B cl to a system of square molecules that, in the absence of such field, have four identical hopping parameters t j,j+1 = t ∀j, as detailed in Appendix C. Note that, despite the presence of B cl , the Hamiltonian in this case respects TRS.This occurs for a special magnetic field, B cl = ±B π ẑ, with B π = cπ/(eA p ). Consider, indeed, the case in which one is in the presence of B cl = ±B π ẑ.
In both cases, it is easy to check that one gets Θ = 0 and t ′ = −t.Now, applying TRS has the net result of changing the classical field from We conclude that at the fields B cl = ±B π ẑ one gets a Hamiltonian that is invariant under TRS.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 where ω c is the cavity photon energy (ℏ = 1), and the Peierls phase θ j,j+1 = (−e/c) rj+1 rj A(r) • d 2 r is necessary to minimally couple the matter degrees of freedom living on the plaquette to the cavity field [52].Here, for the sake of simplicity, we consider a system composed of planar plaquettes with all the same orientation, so that the tunneling coefficients are independent of the specific k molecule.In principle, a system of molecules with the same orientation can be realized by growing a molecular crystal (see, e.g., Ref. [53]).
In the symmetric gauge, the vector potential of the single photon mode is A(r) = −By x/2 + Bx ŷ/2, with ∇ r × A(r) = B = B ẑ. Quantization of the cavity field is carried out in the usual manner by promoting B from a cnumber to a bosonic quantum operator B → B 0 (â + â † ), where â † (â) creates (destroys) a cavity photon.The final Hamiltonian describing light-matter interactions in the cavity is: where λ = −2π[Φ/(Φ 0 n s )] √ N is a dimensionless lightmatter coupling constant, which is proportional to the ratio between the magnetic flux Φ ≡ B 0 A p piercing the plaquette of area A p and the flux quantum Φ 0 = 2πc/e.For example, for a square plaquette (n s = 4) of side d, A p = d 2 .Note that i) in the thermodynamic N → ∞ limit, λ is independent of N , since B 0 scales as ∼ 1/ √ N to make sure that the magnetic field B = B 0 (â + â † ) is an intensive quantity (i.e. it does not scale with N ); ii) the physical flux is Φ = Φ(â + â † ).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 electrical dipole approximation [24,54,55]; iii) For spinless fermions, TRS operates as a complex conjugation on electronic operators.For photonic operators, TRS changes the sign, i.e. â → −â, in order to reverse the direction of the magnetic field.If all hopping terms t j,j+1 are real, then the total Hamiltonian given in Eq. ( 4) is invariant under TRS.

III. MEAN-FIELD THEORY AND PHOTON CONDENSATION CRITERION
To the end of studying the possible emergence of photon condensation, we approximate the ground-state of Ĥ as a product state of the form |Ψ⟩ = |ψ⟩ |ϕ⟩, where |ψ⟩ and |ϕ⟩ are the matter and light quantum states [56][57][58][59][60], respectively.In the thermodynamic limit, the photonic state can be considered a coherent state â |ϕ α ⟩ = α √ N |ϕ α ⟩.Photon condensation occurs when the photonic order parameter α ≡ ⟨ϕ α |â|ϕ α ⟩/ √ N acquires a non-zero value in the thermodynamic limit.From now on, we will take α ∈ R, without loss of generality.We hasten to emphasize that our order parameter α is not plagued by any gauge ambiguity [61] since it physically corresponds to the magnetic flux-defined as the expectation value of Φ over the ground state |ϕ α ⟩, i.e. 2αΦ √ N = ⟨ϕ α | Φ|ϕ α ⟩-which is measurable e.g.via SQUID or NV-center magnetometry [62,63].
At the onset of the phase transition, α is a small parameter and the mean-field Hamiltonian can be therefore expanded in a power series of α, retaining only terms up to O(α 2 ).Using the magnetization operator derived in Appendix B, we get ĤMF (α) where Ĥ0 = Ĥe /N , and we have introduced the paramag- contributions to the magnetic moment operator.We emphasize that since B 0 ∼ 1/ √ N and Mz (α) ∼ N , both Mp (α) and Md (α) are intensive quantities.In Appendix B we show that Mp (α) is proportional to the paramagnetic current operator.
As detailed in Appendix D, 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, ⟨ψ 0 (0)| Mp (0)|ψ 0 (0)⟩ ̸ = 0-where |ψ 0 (0)⟩ is the ground state of Ĥ0 with eigenvalue ε0 (0)-the coupled light-matter system is always in a photon condensate state for any value of the lightmatter coupling λ, and the real ground state energy shift is linear in α.The physical reason for this phenomenon is that, since ⟨ψ 0 (0)| Mp (0)|ψ 0 (0)⟩ ̸ = 0, the matter ground state |ψ 0 (0)⟩ carries a persistent current, which, in turn, creates a non-zero magnetic field in the cavity.This selfgenerated 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 Θ (see Appendix C), explicitly breaking TRS.From now on, we will focus only on Θ = 0; (ii) Conversely, as further discussed in Appendix D, if there are no circulating currents in the uncoupled matter system, i.e. if Θ = 0 and ⟨ψ 0 (0)| Mp (0)|ψ 0 (0)⟩ = 0,
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 χ 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 spatiallyuniform fluctuating B field can be easily encoded in the right-hand side of the inequality (9), changing ω c into a more complicated electromagnetic form factor.Our criterion (9) provides guidance in the experimental search for photon condensation in polaritonic chemistry [40], emphasizing that the quest for this exotic state of matter should focus on the combination between molecular systems with a positive orbital magnetic susceptibility χ M and chiral cavities.

IV. VARIATIONAL THEORY OF THE PHOTON CONDENSATE STATE
The perturbative approach for α ≪ 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 α and the spectra ε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 α by imposing that the optimal value α yields a minimum of the mean-field ground-state energy functional E(α) ≡ ⟨ψ 0 (α)| ĤMF (α) |ψ 0 (α)⟩ /N .Hence, we determine α by imposing that ∂ α E(α)| α= α = 0: The solution of Eq. ( 10) not only determines whether the system is in a normal (α = 0) or photon condensate (α ̸ = 0) phase but also yields α as a function of the microscopic parameters τ , Θ, and λ.If Θ = 0, the photon condensate state ( α ̸ = 0) spontaneously breaks TRS due the appearance of a finite magnetic field and circulating currents.Conversely, in the case Θ ̸ = 0, the system is not invariant under TRS to begin with.Results for α as a function of τ /ω c and Θ for λ = 0.3 have been reported in Fig. 1.
Figure 2(a) shows the molecular spectra ϵ l (α) as functions of τ /ω c .The dashed lines describe the molecular spectra in the absence of cavity (i.e. for λ = 0), while the solid lines describe the case of a finite light-matter coupling (λ = 0.3).The vertical dotted line marks the critical value of τ beyond which a transition to a photon condensate state occurs.We clearly see that at this value of τ the spectra are largely affected by the cavity.Fig. 2(b) shows the order parameter α (solid blue line) and the orbital magnetic susceptibility χ M (red dashed line) as functions of τ /ω c .In agreement with Eq. ( 9), photon condensation (i.e.α ̸ = 0) occurs when χ M /ω c > 1. Diamagnetism (paramagnetism) corresponds to χ M < 0 (χ M > 0).Note that χ p (0) → −∞ (therefore yielding χ M → +∞) due to the degeneracy ϵ 1 (0) = ϵ 0 (0) shown in Fig. 2(a) (see blue and orange dashed lines at τ = −ω c ).

V. POLARITONS
Measuring directly the molecular spectrum ϵ l (α) 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 [65][66][67][68].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 â → α√ N + δâ, where δâ describes a zero-average fluctuation around the mean-field solution α√ N , which has been described above.We then introduce the collective "bright mode" creation operator [69][70][71][72][73][74], with l > 0, where φ † l,k ( φ0,k ) creates (destroys) an electron in the eigenstate |φ l (α)⟩ k (|φ 0 (α)⟩ 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 ϵ l (α).The associate transition energy is ϵ l (α) − ϵ 0 (α).In the thermodynamic N → ∞, b † l behaves as a quasibosonic operator [69][70][71][72][73][74], i.e. [ bm , b † l ] ≈ δ m,l .Further details are given in Appendix F. Hence, in the thermodynamic limit, and expanding the shifted Hamiltonian up to the second order in the photonic fluctuations δâ and in the bright mode bosonic operators b † l , we can write an approximate polaritonic Hamiltonian Ĥpol , describing the lowest excited states of the coupled cavity-molecular system: where Finally, the polariton frequencies can be derived by diagonalizing the Hopfield matrix Ξ, i.e. the matrix that represents the quadratic polaritonic Hamiltonian in Eq. (12).Further details are reported in Appendix F.
Figure 3 shows the four polariton frequencies Ω p as functions of λ and for a fixed value of τ .We clearly see that, at λ = λ c , the lowest polariton mode softens, signaling that the transition to a photon condensate state is a second-order quantum phase transition.When λ exceeds a critical value λ c , α increases from zero to a finite value.Physically, this means that at λ > λ c , a magnetic flux appears spontaneously.
The physics discussed so far does not require large values of the light-matter coupling λ.Although reaching the ultra-strong coupling regime [39,40,75,76] 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 λ → 0 limit.In Appendix F, we show that our model displays such a pleasant feature, provided that one chooses τ = −ω c .Indeed, for photon condensation to occur is more important to hunt for molecular systems with a large value of χ M > 0 and design cavities with a suitable electromagnetic vacuum structure so that χ M /ω 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 [25][26][27]) provided that magnetic effects beyond the electrical dipole approximation are taken into account.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 χ 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 complex chiral cavities [51], which can be studied with recently developed ab initio numerical approaches [43].Finally, ultra-strong magnetic coupling between magnons and a planar superconducting resonator [77] has been recently demonstrated [53].This could be a promising platform to test our findings.Results similar to ours, which show the importance of a cavity with a significant magnetic component, have been also found in a two-leg ladder model: see Ref. [78].

VI. 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 "Re-search 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

Appendix A: Derivation of the effective electronic Hamiltonian
A (single) molecule is composed by mutually interacting electrons and nuclei.Their Hamiltonian Ĥ includes their respective kinetic energies and all Coulomb interactions among them (electron-electron, electron-nucleus, and nucleus-nucleus).In the non-relativistic limit, this Hamiltonian reads: where: In these equations, r i ( pi ) and R I ( PI ) are the position (momentum) of the i-th electron and the I-th ion, respectively.
Given the large difference in mass between electrons and nuclei, M I ≫ m, the dynamics of electrons in a molecule is substantially faster than that of the nuclei.Due to this separation of time scale, we can invoke the Born-Oppenheimer approximation, which assumes that the nuclear motion and the electronic motion can be decoupled and, therefore, the total wavefunction of the system Ψ(r, R) is assumed to be a product of an electronic wavefunction ψ(r; R), which depends on the nuclear positions, and a nuclear wavefunction χ(R), Within this approximation, when we deal with electrons we can ignore the motion of the nuclei, and assume that they see a fixed arrangement of charge.Hence, the electronic wavefunction satisfies the resulting electronic Schrödinger equation: where E e (R) is the electronic energy corresponding to a given nuclear configuration R. The function ψ(r; R) can be thought of as the ground state solution of the electronic Schrödinger equation for fixed nuclear positions.In this expression, the electronic Hamiltonian Ĥe (R) is only parametrically dependent on the static nuclear configuration defined by R, which can be treated as classical parameters.
Lastly, we project the electronic Hamiltonian Ĥe (R) on a basis of localized wavefunctions in order to get a simplified tight binding Hamiltonian.Given Ĥe (R), we can project it onto the localized basis |j⟩ by computing the matrix elements of Ĥe (R) in this basis: Here, H j,j ′ represents the matrix element of the electronic Hamiltonian Ĥe (R) projected on localized states |n⟩ and |m⟩.If the states |j⟩ and |j ′ ⟩ are well-localized, H j,j ′ would be zero for most pairs (j, j ′ ), leading to a sparse Hamiltonian matrix -a key feature of tightbinding models.
In general, the elements H j,j ′ will generally include terms related to the on-site energy of a particle in state |j⟩, and terms related to the hopping of a particle from state |j⟩ to state |j ′ ⟩.Hence, we obtain a tight-binding Hamiltonian in real space for electrons: Here ϵ j is the on-site energy of an electron on site j, t j,j ′ is the hopping amplitude for an electron to move from the site j ′ to site j.Both ϵ j and t j,j ′ can be determined by computing the appropriate matrix elements of the original Hamiltonian, as mentioned above.In what follows, we assume ϵ j = 0.
In second quantization, we can rewrite Eq. (A4) as: where ĉ † j and ĉj are the creation and annihilation operators for an electron at site j, respectively.Eq. ( 1) of the main text is a particular instance of Eq. (A5), providing the addition of a further index k to denote different molecules.

Appendix B: 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 rk , associated with the k-th molecule, can be written in terms of electron creation and annihilation operators, ĉ † j,k and ĉj,k for an electron roaming on a polygon with n s sides: Here, r j = d r [cos(γ j ), sin(γ j )] 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 γ j = 2πj/n 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 [79] Mk = We now consider the magnetic moment along the e z direction, Mz,k , and we use the explicit form of the total Hamiltonian Ĥ in Eq. (B4).We find In the presence of N identical molecules, the total magnetization operator is expressed as 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 |ϕ α ⟩, i.e. onto a state such that â |ϕ α ⟩ = α √ N |ϕ α ⟩.In the thermodynamic limit (N → ∞), this mean-field procedure reduces to the following formal replacement â → α √ N .In this limit, the total magnetization operator becomes where ) To map this operator into a more familiar object with the physical dimensions of a magnetic moment, it is necessary to introduce an effective mass m eff for our electrons hopping through the lattice sites: m eff = 1/(2|t|A p ), where A 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 s = 4), the effective mass is m eff = 1/(2|t|d 2 ).Introducing the effective mass into the definition of the magnetization operator we find where the prefactor has the natural form of a magnetic moment e/(2m eff c).
The total current operator Ĵ can be derived by starting from a discretized form of the continuity equation for the density operator.The local density nj,k ≡ ĉ † j,k ĉj,k obeys the Heisenberg equation of motion By this comparing this equation to the following discretized continuity equation, where Ĵ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, The total current Ĵ is the sum of all local terms, i.e.Ĵ = k,j Ĵj,k , and is given by (B13) As usual, in the mean-field approximation, we replace â → α √ N obtaining a current operator acting only on the matter degrees of freedom: As expected, the current operator Ĵ(α) and the magnetization operator Mz (α) in Eq. (B9) are directly proportional to each other.Specifically, Mz (α) = −(e/c) Ĵ(α)[A p /(n s d)].
Appendix C: The effect of an external classic field As we have seen in the main text, the orbital magnetic response changes sign when τ < 0, i.e. when one of the 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], ĤMF (α) and considering the special case in which all the hopping parameters are equal to t > 0 but for t 0,1 = e −iΘ t, we see that the matter Hamiltonian reduces to By applying the unitary transformation ĉ0,k → e −iΘ ĉ0,k and ĉj,k → e −iΘ(j/ns) ĉj,k for j ̸ = 0, we get This suggests that a classical magnetic field B cl = B cl e z with B cl = −cΘ/(eA p ) can be used to 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.The phase Θ in (C1) can be straightforwardly obtained by applying the Peierls substitution which describes the orbital effect of B cl .Using the vector potential in Landau gauge A(r) = (0, B cl x, 0) T which generates the static magnetic field B cl , in a square plaquette with the four sites located at (0, 0) T , (a, 0) T , (0, a) T , and (a, a) T , one finds that the phase of the link between the sites j and j = +1, θ j,j+1 = (−e/c) rj+1 rj A(r) • d 2 r, is not zero only for the link that connects (a, 0) T and (a, a) T .
The Hamiltonian in Eq. (C2) can be diagonalized by performing a discrete Fourier transformation (DFT) on the index j of the creation and annihilation operators.Let us denote the transformed operators as ĉq,k , which is a function of momentum q and molecule index k.The DFT will be as follows: and its Hermitian conjugate For a system with n s sites, the allowed momentum states are given by q = 0, 1, 2, ..., n s − 1. Substituting the DFT of ĉj,k and ĉ † j,k into the Hamiltonian ĥe,k (α), and performing the sum over j leads to ĥe,k (α) = −2t q cos (Θ/n s + 2λα + 2πq/n s ) ĉ † q,k ĉq,k .

(C5)
The ground state of this Hamiltonian, when it contains a single electron, is denoted by ĉ † q,k |0⟩, where q denotes the momentum corresponding to the lowest energy.In this notation, |0⟩ represents the vacuum state, the state in which there are no electrons.
The energy of this ground state is given by ϵ 0 (α) = min q [−2t cos (Θ/n s + 2λα + 2πq/n s )].The total energy of the system, E(α), is calculated by summing the electronic energy and the energy of the cavity, expressed as E(α) = ϵ 0 (α)+ω c α 2 .The minimum total energy is found by setting the derivative of the total energy with respect to α to zero, at α = α, leading to the following condition: An important observation is that, when Θ = 0, α = 0 is the only solution.This means that if all the hopping parameters are taken to be equal, the phenomenon of photon condensation does not occur in this model.
where m is the electronic mass, pk is the momentum of the k− electron, V (r) is a external potential and v(r k − r′ k ) is the electron-electron interaction.The coupling to a uniform magnetic field B is made as customary via the minimal coupling substitution, i.e. pk → pk +(e/c)A(r k ).
In the symmetric gauge, the vector potential can be expressed as follows where B = Be z is the magnetic field.Notice that the symmetric-gauge vector potential obeys also the Coulomb condition, i.e. ∇ • A(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 B: where L = N k=1 (r k × pk ) 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 M is related to the angular momentum and magnetic field by the following relation: 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 into two terms: ∆E = ∆E P + ∆E D .Introducing the exact eigenstates and eigenvalues of the Hamiltonian Ĥ0 , we find: where Md,z is the e z component of the vector Md and ∆E B is the magnetic contribution, while by construction ∆E D > 0 and ∆E 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 ≈ 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 ∆E P = −(1/2)χ p B 2 where χ p is the magnetization response function, Thus, orbital paramagnetism is governed by the angular momentum-angular momentum response function.From the expression for the response function in Eq. (E8), it can be seen that non-zero orbital paramagnetic response occurs only if the system does not have rotational invariance around the ẑ axis.Indeed, if [ Ĥ0 , ẑ • L] = 0, one can choose a common eigenstate basis of the energy and the angular momentum projection ẑ • L and the matrix element ⟨ψ 0 | ẑ • L |ψ n ⟩ vanishes, precluding the presence of Van Vleck paramagnetism in rotationally invariant systems, such as closed shell atoms.

Appendix F: 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, â = α √ N + δâ, where δâ describes zero-average fluctuations around the mean-field solution α √ N .Hence, in the thermodynamic limit, and expanding the shifted Hamiltonian up to second order in the fluctuations δâ, we find the following Hamiltonian: It is now useful to introduce fermionic operators, φl,k and φ † l,k , that destroy and create an electron in the energy eigenstate |φ l (α)⟩ k of the Hamiltonian ĥe,k (α), defined in Eq. ( 6) of the main-text, i.e. ĥe,k (α) |φ l (α)⟩ k = ϵ l (α) |φ l (α)⟩ k .Notice that φl,k depends on the photonic mean-field α and we dropped that explicit dependence for the ease of readability.Here, the eigenvalues ϵ l (α) are k-independent, since all molecules are identical.Ex-pressing Ĥ in terms of the operators φl,k , we obtain In writing Eq. (F2) we used that the matrix elements the Hamiltonian takes a more compact form By following Ref.[80], we focus on the symmetric Hilbert subspace, which is spanned by the occupation number states defined as  We note that for α = 0, employing Eq. (F34), the previous equation coincides with the instability criterion shown in Eq. (D9).
In Fig. 4 we show the polariton eigenenergies, by diagonalizing the Hopfield matrix expressed in Eq. (F28).In this figure, we illustrate the same quantities as in the analog figure of the main text but for a larger value of τ , i.e. for τ = −ω c .In this case, polariton softening occurs in the weak-coupling λ → 0 limit.The reason is easy to understand.The criterion (D9) for the occurrence of photon condensation in molecular systems we derived in this Letter depends on the intrinsic orbital magnetic response χ M of the molecular system, i.e. on the dependence of χ M on the microscopy of the molecular Hamiltonian ĥe,k in the absence of light-matter interactions.For τ = −ω c , the susceptibility χ p (0) diverges due to a zero in the denominator of a term in the sum given in Eq. (D5).This divergence arises from the double degeneracy of the ground state of the molecular system

FIG. 1 .
FIG. 1. (Color online) (a) Pictorial representation of the setup where many planar square molecules interact with a cavity's single-mode, quantized magnetic field, B. The field is perpendicular to the molecular plane.Due to the small spatial extension of the molecule cloud compared to the field's wavelength, the magnetic field is assumed constant and the electric field negligible.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 = ωc and λ = 0.3.The photon condensate order parameter α (color scale) is plotted as a function of the two microscopic parameters τ ∈ [−t, t] and Θ ∈ [0, π[.

FIG. 2 .
FIG. 2. (Color online) Photon condensation of molecules in a chiral cavity.Results obtained by setting Θ = 0, λ = 0.3, and t = ωc.(a) Eigenvalues ϵ l ( α) of the single-molecule Hamiltonian ĥe,k plotted as functions of τ /ωc.Solid lines denote results for λ = 0.3.Dashed lines denote results in the absence of the cavity.The vertical dotted line marks the critical value of τ at which a transition to a photon condensate state occurs.(b) The photon condensate order parameter α is plotted as a function of τ .When τ reaches a critical value, α becomes finite, signaling a quantum phase transition to a photon condensate state.

FIG. 3 .
FIG. 3. (Color online) Polariton softening at finite λ. Results in this figure have been obtained by setting with Θ = 0, τ = −0.5ωc,and t = ωc.The four polariton energies Ωp are plotted as functions of λ.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.
We use the Heisenberg equation of motion to find an expression for the velocity operator ṙk ,ṙk = i[ Ĥ, rk ] .(B3)Replacing this result in the definition of the magnetic moment (Eq.(B2)