Cavity-induced bifurcation in classical rate theory

We show how coupling an ensemble of bistable systems to a common cavity field affects the collective stochastic behavior of this ensemble. In particular, the cavity provides an effective interaction between the systems, and parametrically modifies the transition rates between the metastable states. We predict that the cavity induces a phase transition at a critical temperature which depends linearly on the number of systems. It shows up as a spontaneous symmetry breaking where the stationary states of the bistable system bifurcate. We observe that the transition rates slow down independently of the phase transition, but the rate modification vanishes for alternating signs of the system–cavity couplings, corresponding to a disordered ensemble of dipoles. Our results are of particular relevance in polaritonic chemistry where the presence of a cavity has been suggested to affect chemical reactions


Introduction
Classical rate theory provides a stochastic framework for investigations in natural and social sciences [1][2][3][4][5][6].It has been used to describe a wide variety of phenomena in different fields: the spreading of diseases in epidemiology [7,8], the growth and decline of companies in economics [9] and of populations in ecology [10], and disintegration of radioactive nuclei in atomic physics, to name a few.In chemistry, rate equations provide a quantitative way of understanding chemical reactions [11].
There have been several recent indications of chemical reactions being altered by the formation of vibrational polaritons, hybrid excitations formed by the coupling of molecular vibrations to vacuum electromagnetic field of an optical cavity [12][13][14].The vacuum field has been reported to work both as a catalyst [15,16] and an inhibitor [17][18][19][20] in different reactions.However, the current theoretical understanding seems to contradict these experimental results.The conventional approach that concentrates on the cavity-induced modification of the individual reaction rates within the transition state theory finds that all the polaritonic effects should disappear when the number of molecules forming the polariton increases [21][22][23][24][25]. Furthermore, very recently, there have been reports of failed replication experiments [26,27], adding to the confusion.
Regardless of the current experimental status, vacuum-modified chemistry poses a fundamental challenge on our understanding of light-matter interactions.Its difficulty arises from the notion of collectivity.The cavity fields interact with a large number of molecules which forms an interacting many-body system.The collectivity of the light-matter system leads to novel optical properties [28], including the realization of hybrid Bose-Einstein condensates [29,30] and sub-/superradiant quantum emitters [31], as well as novel interactions whose role in materials design is actively studied [32][33][34][35][36].
In this article, we predict a phase transition induced by the presence of the cavity.This transition is caused by the indirect interactions within the matter system, e.g.molecules, mediated by the vacuum field.We show that such an interaction can be taken into account in the classical rate theory as transition rates that depend parametrically on the state of the molecules.We find an order parameter for the light-matter system and show that, below a certain critical temperature which depends linearly on the number of molecules, it bifurcates so that it can have multiple solutions in the stationary state.The phase transition aligns the molecular states depending on their coupling to the cavity.If the molecules couple equally, they will spontaneously prefer being in the same state.When this is not the case, there is a phase Figure 1: Interactions between a harmonic cavity mode x and N identical modes q i , described by a bistable potential V (q).In the harmonic approximation, ∂ 2  q V (a/2) = ω 2 R and ∂ 2 q V (−a/2) = ω 2 L define the relevant frequencies.
separation based on the light-matter coupling constants.Even before the phase transition, the cavity-mediated interaction may slow down the effective transition rates, suggesting that the collectivity of the many molecule system must be considered on the level of the rate equations, not only on the level of the rates themselves.

Rate theory
Let us consider a one-dimensional and bistable system described by a potential V (q) as in Fig. 1.
Here, q is a position quadrature which can describe, e.g., vibrational modes of a molecule in its electronic ground state.Classical rate theory then relates the probability p to find the system in one of the wells to the transition rates Γ µν between the wells.We adopt a terminology of left and right wells corresponding to the minima of V (q) at q < 0 and q > 0, respectively.The state of N bistable systems can be described with the vector (s 1 , s 2 , . . ., s N ) where s j = sign q j .Thus, if p j is the probability to find the jth system in the left well (s j = −1), its rate equation 1reads as d dt The macroscopic behaviour is obtained from Eq. ( 1) by noting that the expectation value of the number of systems in the left well is given by N L = j p j .If the transition rates are independent of the system j, we find d dt In the following, we show that this approach has to be extended when the molecules have varying couplings to the cavity.
Often, the transition rates are assumed to be independent of p j and time, i.e., constants.Physically, this is motivated by separation of time scales.We assume the transitions to happen so rarely that the molecules manage to relax and their environment lose any information before the next transition [11].We consider the thermal activation regime, the ratio between the potential energy barrier and the temperature k B T controls the rate of transitions and should then be sufficiently large.This is in contrast with the quantum tunneling regime where the energy barrier height is compared to the resonant frequency inside the potential well [25].In such a case, to comply with thermodynamics, the rates must obey the detailed balance Γ LR /Γ RL = e βσ where β = 1/k B T is the inverse temperature and σ = V (−a/2) − V (a/2) corresponds to the potential bias between the different states (Fig. 1).The rate equation then describes thermalization to the Boltzmann distribution with p j = N L /N = 1/(1 + e βσ ).
The coupling of the molecules to the vacuum field leads to a state-dependent modification of the rates in Eq. ( 1).This is caused by the effective interaction between the different molecules mediated by the cavity.Before deriving the modification explicitly, we write this statement formally as Γ

State-dependent rate modification
Consider N systems with coordinates q j coupled to an external collective mode as in Fig. 1.In what follows, we call those N systems "molecules" and the collective mode is a "cavity", in analogy to the systems studied in polaritonic chemistry.However, this approach also works, for example, in the case of a large number of superconducting flux qubits coupled to a single microwave cavity, or bistable atoms in a cold-atom arrangement, coupled to a common cavity mode.These potential realizations are discussed in Appendix G.
We derive the state-dependent transition rates from the detailed balance.That is, if the stationary state of the total system thermalizes in the presence of the cavity, the rates for the molecule j obey when V tot is the potential energy for a given state of the molecule-cavity system.The quantity in square brackets represents the energy cost of moving one molecule from right to left well.Motivated by quantum electrodynamics (QED) in the dipole gauge [23,[37][38][39][40][41], we introduce a single cavity mode of frequency ω c and its electric field quadrature x and choose the classical potential energy of the system2 to be similar to that producing polariton excitations, where d j (q j ) represents the molecular dipole moment component in the direction of the cavity mode's polarization vector.In the Born-Oppenheimer approximation, d j (q j ) is proportional to the expectation value of the dipole operator in the electronic ground state with a given value of q j .We note that the energy does not include the so-called dipole self-energy terms, proportional to [d j (q j )] 2 in the long-wavelength limit, similar to other recent works in strong coupling regime [20][21][22]41].Moreover, the model containing only a single cavity mode is chosen for simplicity.Generalizing the results to the case of several cavity modes is relatively straightforward and does not lead to qualitative changes in the results.Quantitatively, this would amount to calculating the Casimir-Polder-van der Waals [42] interaction between the dipoles in a cavity.Note that such an interaction is present also without the presence of a cavity: we assume that without the cavity the interaction is to be too small to lead to the collective effects discussed here.We focus on a linear dipole moment d j (q j ) = ω c ω m g j q j to gain insight on the effect of the cavity on the transition rates.Here, ω 2 m = (ω 2 L + ω 2 R )/2 represents the mean frequency of the potential V (q), and the square root term normalizes the light-matter coupling constant g j .These choices allow for mapping this model potential to the Dicke Hamiltonian in which case g j follows directly from the QED derivation [43].However, our approach works for an arbitrary dipole function d j (q j ) with the linear approximation d j (q j ) ≈ d j (±a/2) + d ′ j (±a/2)[q j ∓ a/2] near the minima of V (q j ).The important parameters are then the transition dipole moments d ′ j (±a/2), responsible for coupling molecules in left and right wells to light, and the difference in the dipole moments d j (a/2) − d j (−a/2) between the metastable states.For a linear dipole moment, these three parameters are determined by g j , simplifying the subsequent analysis.
The light-matter coupling g j depends on both the position of the molecule within the cavity and the alignment of its dipole moment with respect to the cavity polarization vector [43][44][45].Importantly, g j can be either positive or negative which has physical consequences for the many molecule system.The strength of the collective coupling is indicated by the Rabi splitting frequency Ω ∝ j g 2 j , defined as the difference of the eigenfrequencies of the total potential on resonance ω L = ω R = ω c [46,47].We assume that g j ≪ ω c , ω m while Ω is a small fraction of ω c and ω m due to N ≫ 1.
The effective interaction mediated by the cavity can be understood by a force analysis.For a system in a stationary state, displacing the mode q j by δq j causes a force proportional to −g j δq j to the cavity mode x.This leads to a displacement δx proportional to the force to the cavity mode.The force applied on another mode q k is then proportional to −g k δx which reads as g k g j δq j in terms of the original displacement on mode q j .Whenever the coupling constants of two different molecules share the same sign, the force between them is attractive and for opposite signs repulsive.The cavity-induced interaction resembles the phonon-mediated interaction between electrons in superconductivity [48].
We solve the stationary states of V tot (x, q 1 , . . ., q N ) within the harmonic approximation which holds when g j ≪ ω c , ω m and, thus, the cavity-induced change to the local minimum points of V (q j ) is small.Technically, we first solve the values of q j 's and x from the conditions ∂ x V tot = ∂ q j V tot = 0 by treating the signs of q j as variables, representing whether the molecule j is in the left or the right well.We find that all the stationary points may be expressed by using the collective order parameter where g 2 = i g 2 i /N represents the average over the molecules and s i = sign(q i ).If there is no variation in the coupling constants g i , ∆ becomes the normalized difference between the number of molecules in the right and the left well, ∆ → δN /N = (N R − N L )/N .Finally, we calculate the value of V tot (x, q 1 , . . ., q N ) using Eq. ( 4) at these stationary points.The details of this algebraic calculation with an arbitrary dipole moment function d j (q j ) are relegated to Appendix A. We note here that the order parameter ∆ has the same form if one replaces g i with state-dependent light-matter couplings g i (s i ) arising from transition dipole moments d ′ (±a/2), but the effects reported below disappear if the dipole moment is the same on both states, i.e., d i (a/2) − d i (−a/2) = 0.
We find that the stationary states of the full N molecule system have the potential energy The first term is caused by the effective interaction via the cavity mode while the second term describes the potential bias of V (q).The magnitude of the interaction is determined by the product of an effective potential barrier E b = 1 2 ω 2 m (a/2) 23 and Here, ω(s i ) equals ω L for s i = −1 and ω R for s i = +1.The coefficient P can be related to the Rabi splitting Ω using N g 2 ∝ Ω 2 .We assume that Ω is a small fraction of ω c ω m and P ≪ 1.Consequently, P may be treated as a state-independent constant since its variation is proportional to P 2 as shown in Appendix B.
Using the detailed balance (3) and Eq. ( 6), we find in the lowest order of P that where α = 2Pβ E b plays the role of a control parameter.It should be noted that, even though P ≪ 1, the effective potential barrier E b can be much larger than the temperature so that α can be of the order of unity.Equation ( 8) expresses the relative change of the rates due to the change of the number of molecules in the left and right wells.To find the individual rates, we use the fact that reversing all coordinates also interchanges the rates.That is, we impose that r j RL must be equal to r j LR if we at the same time transform σ → −σ, ω L/R → ω R/L , and q i → −q i for every i.The control parameter α retains its value in this parity transformation.By using this equality in Eq. ( 8) and expanding the exponent of r j LR in the powers of ∆, we find that where f is a function of ∆ that stays invariant in the parity transformation.A few examples of such functions are ∆ 2 and σ∆.C is a constant that describes a state-independent cavityinduced modification that may be acquired in, e.g., the transition state theory.Here, we assume that such modifications are small and C = 1.The function f , although not forbidden by the symmetry, we neglect as there is no clear physical origin for such terms.That is, we set f = 0.In this case, r j RL = 1/r j LR .Next, we use Eq. ( 8) to evaluate the stationary state of the macroscopic system.We note that our approach is consistent with thermodynamics by construction, and hence the stationary properties could also be found by treating the light-matter system as a canonical ensemble and solving its thermal distribution.The separate rate modifications (9) allow for a continuous-time Markov chain simulation of Eq. ( 1) from which dynamics can also be investigated and compared to an analytical approach. 4We further elaborate the numerical method in Appendix E.

Cavity-induced phase transition and collective dynamics 4.1 Equal light-matter coupling strengths
We illustrate some qualitative consequences of the parametric rate modifications by taking a simplifying limit of identical couplings as in [46]; we set g j = g 0 for all j.In this case, the  2) exhibiting bifurcation for equal couplings.The dotted line represents an unstable solution whereas the solid lines are stable; the average tends to one of the stable solutions after a long time.(b) Phase separation of a polaritonic system that is divided into two equally large partitions (N 1 = N 2 = N /2) with couplings related by g 1 = −2g 2 .The solid lines represent the dynamics of the partitions alone whereas the dashed line describes the dynamics of the full system.Here, the potential V (q) is symmetric with Γ 0 LR = Γ 0 RL and σ = 0.The system has undergone the phase transition as α = 1.3.transition rates become independent of the molecule in question.The set of rate equations can be mapped to a master equation describing the probability to find exactly N L molecules in the left well as shown in Appendix D. Interestingly, a similar state-dependent master equation has been used to model the formation of time crystals in magneto-optical traps [49].We confirm numerically that such a master equation and the macroscopic rate equation (2) produce the same results as long as N ≫ 1.
We find the stationary state of the full N molecule system by inserting Eq. ( 8) into the macroscopic rate equation (2) together with d dt N L = 0.It gives The solutions for δN determine the long-time behaviour of the macroscopic system.First, let us consider σ = 0.If there is no coupling to light, the only stationary state is N L = N R = N /2 due to the left/right symmetry.For small values of α = 2Pβ E b there is no change to this stationary state.However, as soon as α > 1, we find two new stationary states while the state N L = N R becomes unstable.The stationary states are represented in Fig. 2(a) as a function of the control parameter α.As this configuration is possible only if the coupling is large enough, we call it the cavity-induced bifurcation of the stationary state.
The cavity-induced bifurcation remains in the presence of a finite bias σ.However, the critical value of α increases as σ increases.For a very large bias βσ ≫ α, all the systems eventually end up in the right well which is lower in energy.This is expected, as biasing the potential pushes the molecules towards the right well whereas the attractive interaction tries to keep the molecules in the same well.Before the bifurcation, one can also see the effect of the cavity-mediated attraction in that more molecules can be found in the right well than without the cavity.This seems as if the cavity increased the potential bias.
The mere state-dependent light-matter coupling -caused by a nonlinear d i (q i ) -may bias the system.These effects are described at length in Appendix F. Similar to a potential bias, the molecules prefer the state with the larger coupling strength even with a single stationary state, and the larger the coupling difference, the larger the critical value of α.

Effect of coupling strength distribution
If there is only a partial control of the coupling constants as is typical in polaritonic chemistry, one must include the variation of the couplings.We argue that this can but does not necessarily remove the cavity induced effects.Especially, the role of the average 〈g〉 is important.This average is independent of the optically observed Rabi splitting which is proportional to N 〈g 2 〉 [12-14, 50, 51] and hence does not depend on the sign of the individual g j 's.
The order parameter ∆ determines the individual rates in Eq. ( 8).If we consider a partition of the N molecules such that all molecules within the partition share the same coupling constant g j , they also share the same transition rates Γ j µν .Macroscopic rate equations similar to Eq. ( 2) can be derived for each partition but deriving a closed-form differential equation for ∆ is impossible.Consequently, the full macroscopic rate equation cannot be evaluated.However, the stationary states of the partitions can be solved in terms of ∆.As shown in Appendix C, this leads to a self-consistency equation whose solution can be used to find given a distribution of coupling constants g.We specifically discuss the case σ = 0.By expanding Eq. ( 11) to the third order in ∆, one finds the solutions ∆ ∝ ± α − 1 and ∆ = 0. Thus, α = 1 is the critical value for the bifurcation of the order parameter ∆.This agrees with Fig. 2 in which ∆ = δN /N .The critical temperature of the phase transition is approximately given by At the same time, Eq. ( 12) may be linearized giving δN ∝ 〈g〉∆.Whenever 〈g〉 = 0, we find δN = 0 or N L = N R .Thus, even if ∆ bifurcates, it is possible that no change in the macroscopic state is seen.The bifurcation of ∆ indicates formation of a new phase; the coupling constants g j and the corresponding modes q j become strongly correlated.This is exemplified in Fig. 2(b) in a simple case of two equally large partitions with two different couplings.The partitions have distinctly different dynamical behaviors and, in the stationary state, the molecules are found mostly in the right well for one partition and in the left well for the other.Such a correlation may be observed even if all the couplings are different.

Change to thermalization or reaction rate
To conclude the analysis of the rate theory approach, let us consider the dynamics of thermalization, i.e., reaction kinematics.To this end, consider an experiment where every molecule is initialized in the left well ("reactant") and, at a time t later, the percentage of the molecules in the right well ("product") is measured.The thermalization (or reaction) rate is initially controlled by Γ LR alone.Averaging over Eq. ( 9), this initial rate is modified from Γ 0 LR by a factor in the leading order of 1/N and to the second order in the cumulant expansion.This result, derived in Appendix B.1, is independent of the bias σ.In the inset of Fig. 3, we show the effect of α and the variance of the coupling distribution.The effective transition rate becomes always slower due to the presence of the cavity -the cavity is an inhibitor -and the effect is stronger the less varied the couplings are.The agreement between our theoretical approach and the Markov chain simulation is shown in Fig. 3.For equal couplings, Var(g) = 0, we use the macroscopic rate equation (2) directly and the results are nearly indistinguishable.Our solutions for the stationary and transient behavior also agree well with the simulation.

Phase transition from a quantum mechanical model
Previous section describes how the cavity-induced interaction between molecules leads to a bifurcation in the rate equation model, resembling a phase transition.This phenomenon has some similarities to the superradiant phase transition in the Dicke model, originally discussed by Hepp and Lieb [52].Namely, a collection of electric dipoles interacting with the electromagnetic vacuum field can be exactly mapped to an Ising model that exhibits a ferromagnetic phase transition.However, while these two phase transitions have some similarities, they are not the same as they exhibit a different set of parameters.
Consider N two-state systems having a dipole operator proportional to σ z .Let us assume that there is no energy difference between the two states.A conventional (low-energy) quantum electrodynamical Hamiltonian to describe the interaction of these dipoles with a cavity or vacuum field reads as (ħ h = 1) [43] Here, ω α is the eigenfrequency of the cavity mode α to which photons are created by the operator a † α and from which they are annihilated by a α .The light-matter coupling between these cavity modes and a two-level system j are given by g jα .
A fully-connected Ising Hamiltonian is now found by applying a unitary transform with where we have neglected a constant energy shift in the last step.
The spontaneous symmetry breaking associated with the Ising model corresponds to the alignment of the electric dipoles σ j z and the spontaneous polarization of the cavity field.For the sake of a simple argument, consider only a single relevant cavity mode α to which all the two-level systems couple equally with strength g.The effective Ising Hamiltonian becomes If embedded in a bath of temperature T , the Ising model's phase transition happens when T < T c = (N − 1) , which spontaneously aligns the two-level systems to the same state [53].
Let us compare this to the critical temperature ( 13) obtained within the rate theory analysis.It is analogous in that both contain spontaneous symmetry breaking but, in the case of a semiclassical stochastic process, the potential energy surface increases the critical temperature by a factor of 2E b /ω b compared to the Ising model.In the Ising model literature, the two-level systems describe the spins of a ferromagnet and their alignment corresponds to a finite magnetization.Here, however, the electric dipoles are aligned.As a direct consequence, the electric field of the cavity mode, proportional to a α + a † α , breaks the symmetry as well, obtains a finite expectation value, and the vacuum becomes polarized.
The existence of an equilibrium phase transition in QED involving such a superradiant or a photon condensate state has been debated over the years.Often, this is formulated in terms of no-go theorems which show that such transition cannot happen due to gauge invariance which is broken when the self-interaction term is neglected [54].However, these no-go theorems rely on strong assumptions such as that the molecules are located in a region much smaller than the wavelength of the cavity field [55] and that there are no direct intermolecular interactions [38,56].In the context of polaritonic chemistry, at least, such assumptions do not respect physical reality of the experiments.Furthermore, it has been argued that the direct Coulomb interactions between the dipoles exactly cancel the self-interaction term so that the one is left with the Dicke model as in Eq. ( 15) [38].
The spontaneous symmetry breaking is not a resonant effect, caused by tuning the cavity eigenfrequencies ω α which is often associated with polaritonics.However, it is intrinsically a collective effect.One cannot find this phase transition in a model consisting of a single two-level system.

Conclusion
To summarize, we have shown how a cavity-mediated interaction between molecules can cause a collective change in their stochastic behavior.The cavity may cause a slowdown of the thermalization rate of the full system and a bifurcation in the stationary state.The cavity effect depends on the number of molecules and the distribution of the light-matter couplings.The phase transition can even show up as a phase separation without altering the macroscopic state.The semiclassical model based on classical rate theory that underlies these results can be understood as an effective dipole-dipole interaction mediated by the electromagnetic vacuum field.This we illustrate by showing that a quantum mechanical light-matter Hamiltonian can be mapped to a fully-connected Ising model of interacting dipoles in a special case.
In polaritonic chemistry, our approach is straightforward to extend to studies of multiple cavity modes, molecular symmetries [57], and cavity-induced reaction selectivity [19].Thermodynamic extensions should also open a fruitful avenue of study: Here, the light-matter system forms a canonical ensemble but, instead of connecting to classical rate theory, one may investigate its thermodynamic properties.The phase transition should have associated thermodynamic signatures, such as latent heat, whose magnitude and significance is unknown.Since we expect that, experimentally, the molecules and their dipoles are highly disordered, these thermodynamic signatures might provide the most practical observable to measure while the reaction rates or equilibrium concentrations are unaffected.Using the approach of this paper, it seems possible to consider grand canonical ensembles to describe a much wider variety of chemical systems.Furthermore, it is not clear to us how to incorporate cavity-induced phase transition physics into conventional rate theories used in chemistry, e.g., transition-state theory.
Our work indirectly raises the question whether it is possible to replace the phenomenological description of the light-matter coupling by a microscopic QED description.That is, one might be able to derive an expression for the effective dipole-dipole interaction mediated by the cavity vacuum field in terms of the physical parameters, such as the cavity size, dielectric constant within the cavity, and fine-structure constant.It would elucidate the physical mechanism of the cavity-mediated interaction as well as explain the critical temperature (13) with only the physical parameters of the system.
Due to the ubiquitous nature of light-matter coupling model used here, our results can be used to describe other systems such as electric circuits and cold atoms in optical traps.In these systems the disorder of the effective dipoles can be controlled to a much greater degree than in polaritonic chemistry experiments.This allows for direct observation of the phase transition and its effect on thermalization.

A Calculation of the state-dependent potential energy
Here, we present how we calculate the potential energies of the classical states, starting from the total potential where V (q i ) is a bistable potential and d i (q i ) characterizes the dipole moment of the ith molecule projected on the polarization vector of the cavity mode.We point out that the following calculation is very similar to the Caldeira-Leggett model except that the typical bath of harmonic oscillators is replaced by a single mode.Consequently, instead of dissipation, integration over the cavity mode x leads to effective interaction between the molecular coordinates q i .First, we solve the stationary points dictated by the conditions ∂ x V tot = 0 = ∂ q i V tot .We work here in the harmonic approximation, assuming that the coupling to the cavity quadrature x induces a small change in the positions of the stationary points q * i which are the minimum points of the potential fulfilling ∂ q i V (q i ) q i =q * i = 0. We choose the coordinate system so that q * i = ±a/2 without the presence of the cavity field.Thus, a is the distance between the (local) minimum points and the molecular state can be encoded into s i = sign(q i ) = ±1.To characterize the displacement of the minimum points, we introduce a variable δq i = q i − a 2 s i .We can now write the potential in the harmonic approximation as where σ represents the energy difference between the left and right well and ω 2 i is a shorthand notation for which alternates between ω 2 L and ω 2 R depending on the molecular state.Similarly, we assume that the dipole moment can be linearized near q i = ±a/2 where λ 2 i (s i ) describes the (possibly state-dependent) strength of the light-matter coupling and q 0 additionally characterizes the difference in the dipole moment between right and left well.The coupling constants λ 2 i (s i ) can also vary from molecule to molecule, for instance, due to disorder in the polarization vectors and positions.Here, we choose λ i (s i ) to have the units of frequencies ω i , but it should be noted that λ 2 i (s i ) can be negative.The results of the main paper are obtained by replacing q 0 = a and λ 2 i (s i ) → ω m ω c g i where g i is independent of the molecular state.The presence of the square root term is due to the choice of coordinates q i and x so that, for a harmonic potential V (q i ), we retain the light-matter interaction with constants g i as g i (a † + a)(b † i + b i ) after quantization as in the Dicke model.As the derivative of the sign function s i = sign(q i ) vanishes everywhere except at q i = 0, and q i obtains values only near ±a/2 by assumption, we find This set of equations can be solved by noting that Eq. (A.5a) depends on Q = i λ 2 i (s i )δq i and, at the same time, we can derive from Eq. (A.5b) an equation for the collective variable Q. Inserting the solution of Q in terms of x to Eq. (A.5a) gives In the latter equality, we introduce the useful notations from the main text -which we gather here for convenience (A.8) Though, here, the average coupling λ 4  s still depends on the states of the molecules which is indicated by the subscript s.The solution of x allows for the solutions of individual δq i which are simply δq i = − Note that the solution is consistent with the harmonic approximation as long as the dimensionless constant P is well below unity and Finding the potential minima is now straightforward.We insert the obtained local minimum points into the expression of total potential V tot .The potential energies found in this way still depend on the state of the system in the sense that the variables s i = sign(q i ) are not fixed.Within the harmonic approximation it holds that |δq i | < a/2 and it is clear that the signs of q i indicate whether the system is in the left (q i < 0, s i = −1) or right (q i > 0, s i = +1) well.We find which matches with the result of the main text when 2 is identified in the latter term.Since E b has the dimension of energy, we call it the effective potential barrier.In this article, for the most part, we focus on dipole moment functions d i (q i ) that can be written as a linear function over the whole range of q i .Physically, this means that both molecular states are equally coupled to the cavity.The condition is easily relaxed to investigate a multitude of possible polaritonic systems -we briefly discuss one example in Sec.F. In terms of the model at hand, it means that d i (q i ) = λ 2 i q i where λ 2 i is a constant for each molecule i.Thus, the average λ 4 is state independent and, as seen in the next section, P can be regarded as a constant in the lowest order of the light-matter coupling.In addition, if also ω L = ω R , E b simply describes the potential energy V (q = 0) within the harmonic approximation since q 0 = a.

B Detailed balance and the transition rate modification
The detailed balance states that the transition rates must obey in order for the system to thermalize to the Boltzmann distribution.We assume here that the light-matter coupling is state independent and hence λ 2 i (s i ) = λ 2 i .First, we calculate the difference in P = P(s 1 , . . ., s N ) in Eq. (A.8) due to one molecular state change.It is useful to denote K = 1 N i . Thus, we formally can expand P around K = 1 which leads to Since we assume that N λ 4 /(ω 2 c ω 2 m ) is small and P 0 ≪ 1, P ≈ P 0 + P 2 0 (K − 1) to a good approximation.The higher powers of K also scale in the powers of 1/N .Finding the difference is now straightforward and results in δP j = P(. . ., s j−1 , −1, s j+1 , . . . ) − P(. . ., s j−1 , +1, s j+1 , . . .
We neglect the second order contribution of P 0 in our analysis and, thus, the state dependency of P ≈ P 0 as δP j ≈ 0. Since P can be regarded as a constant in the calculation of the difference in potential energies, the calculation simplifies greatly and we obtain However, we still have to connect this to what would be observed in experiments which deal with macroscopic numbers of N .With the value of δV tot , the detailed balance states that Here, one can readily identify the ratio of the rate modification factors which is the latter exponent.Now, we can appeal to symmetry: if we interchange left and right everywhere in the potential V (q), the rate modifications should remain the same.Thus, where {q i } refers to the set of all q i 's.One can confirm that P remains invariant in this transformation (the relevant quantity in ).Following the line of deduction as in the main text, we can thus write and r j RL = 1/r j LR .We may also replace the sum in this expression with ∆ as the error made is of the order of 1/N .

B.1 Cumulant expansion
Next, we treat the coupling constants λ 2 j as independent random variables -instead of considering them to have some fixed values -and average the detailed balance relation over all realizations.Furthermore, we assume that the molecular state s j is independent of λ 2 j which does not necessarily hold for a dynamical system but may be assumed to hold for the initial state, for instance.The average of the exponent gives rise to the cumulant expansion.To calculate the cumulants, we first calculate the raw moments of the term in the exponent.Taking only the leading order in 1/N , we find Since all cumulants may be expressed using only raw moments in a specific polynomial manner [1], we can deduce from this result that the nth cumulant K n of the term in the exponent of r j LR is That is, we have related the cumulants of the sum to the cumulant of its individual elements.
Using the cumulant expansion to second order, we have that Here, the first term follow from the average and the second from the variance.We emphasize that the approximation is valid when N ≫ 1.The cumulant expansion is otherwise exact if the distribution of λ 2 j 's is Gaussian; only the first two cumulants may be finite for a Gaussian distribution.

C Derivation of the self-consistency equation
L is the number of systems in the left well.Following the main text, we write for each subsystem k the macroscopic equation We assume that each subsystem has a stationary value.We may then solve, using Eq.(B.5), Inserting this solution to the definition of ∆ and δN leads to In both equations, we identified the structure of k (N k /N ) f (c k ) as an average over the distribution of c's which follows from the fact that there are exactly N k values of c k .In the limit N → ∞, we may construct a continuous distribution of the coupling constants which we use here as an approximation for systems of finite size.
It should be noted that we treat ∆ and δN as numbers, not stochastic variables, even though the underlying stochastic processes determine their realizations and the possible values of ∆ are constrained by the distribution of the coupling constants.
In this derivation, we have neglected the possible variation of α with the macroscopic state.This arises from the variation of P, as discussed in Appendix B. In principle, this is remedied by including an equation for α = α(∆) and solving all three equations self consistently.It should be noted that the variation of P with the macroscopic state is a second-order effect in P and, thus, one should also calculate δV tot to the same order.

D Relation to the full master equation of the N molecule system
If all the coupling constants g j are equal, the potential difference δV tot is independent of the system index j.In this case, the full N molecule system is readily described by the probability P(N L ) to find N L systems in the left well.The master equation for these probabilities reads as where µ M = M Γ 0 LR r LR (M ) describes the removal rate of particles from the left well to the right well and ν M = (N − M )Γ 0 RL r RL (M ) the inverse process.Now, the expectation value of N L is given by 〈N L 〉 = N L N L P(N L ).Using this definition, one can derive from the master equation an equation for the average as This matches exactly to the macroscopic rate equation in Eq. ( 2) if , and the resulting rate modification factors may be understood as a function of the mean value of N L only.Due to the non-linearity of the rate modifications, this can be used as a valid approximation only when the fluctuations of N L are small compared to some characteristic scale given by the rate modifications.This assumption can be checked numerically.First, we note that the transition rates are independent of time.The master equation (D.1) can thus be written in a matrix form as ∂ t ⃗ P = W ⃗ P where W is independent of time.Given an initial state ⃗ P(t = 0), we formally have ⃗ P(t) = e W t ⃗ P(t = 0).The numerical problem is then to calculate the matrix exponential e W t .Fig. 4 shows that the macroscopic rate equation and the master equation produces the same mean behavior for N L as N is large enough.
We note that even though the rate equations are state dependent and consequently nonlinear, the macroscopic behaviour fits well to the typical exponential behaviour with a suitable fitting constant γ.This observation seems to hold even when the variation of couplings is included which is possible by the numerical method presented in the next section.

E Numerical Markov chain approach
Here, we shortly describe the Markov chain algorithm we use in this work.The main task is to simulate a set of rate equations where Γ j µν depend on the state of the system.In our formulation, these states are s j = sign q j so that s j (t) ∈ {−1, +1}.The rate equation describes that in a small time step ∆t, the probability to change the system from s j = −1 to s j = +1 is given by Γ j LR ∆t while the converse process has the probability Γ j RL ∆t.The algorithm is as follows: 1. Initialize the states s j at time t if necessary.

For each
+1 using the system state at time t.
4. For each j ∈ {1, 2, . . .N } compare u j and A j .If u j < A j , set the system state at time t + ∆t to be s j (t + ∆t) = −s j (t).Otherwise, set s j (t + ∆t) = s j (t).
5. Return to step 1 with t → t + ∆t.
Thus, when an appropriate amount of time steps is taken, we have a list of system states s j (t) at the chosen time points.These states can be used to calculate e.g.∆ or δN .This is exemplified in Fig. 5. Practically, we choose N to be large enough so that there are no spurious 1/N effects and correspondence with the master equation approach holds.The time step we choose so that Γ 0 µν ∆t is below or of the order of 1 percent.The initial state is chosen according to how many systems are wanted to be found in the left well and then the states are shuffled.
Here, the coupling constants are sampled from a Gaussian distribution with zero mean.a) The time evolution of the quantities δN and ∆.As discussed in the main text, ∆ has a solution other than zero when α > 1 while δN fluctuates around zero. b) As a visual representation, the states "left" and "right" are here represented by blue and orange, respectively, for each time step.Here, the different systems are indexed in the ascending order in the coupling constant (i.e., g 1 ≤ g j ≤ g N ).The correlation becomes visible; the systems with largest values of coupling are predominantly found in the left well after a transient period.

F Effect of nonlinear dipole moment in a symmetric potential
We shortly describe how the physics changes when the assumption of equal light-matter coupling strength in both wells is lifted.Now, we have to specify two couplings and λ 2 i (s i = +1) = λ 2 iR for each molecule in Eq. (A.4) based on its state.Furthermore, the average coupling on different states may differ, that is For simplicity, we assume here a symmetric potential with ω L = ω R and σ = 0.
Physically, breaking the (anti)symmetry of the dipole moment leads to a preference or a lower energy on the state that is coupled more strongly to light.This is relevant for engineering state or reaction selectivity by vacuum fields.The bifurcation behavior presented in the main text also changes.Here, we investigate this numerically with the help of the Monte Carlo algorithm presented in the previous section.
The main mathematical difficulty in the case of a state-dependent light-matter coupling is that the quantity P is no longer a constant.For a linear dipole moment, this dimensionless parameter describes the strength of the Rabi splitting.Now, assuming relatively small couplings and thus taking only the lowest order in 〈λ 4 〉 s /(ω c ω m ), it is useful to write P in terms of an auxiliary constant λ4 Recall that λ 4 s depends on the state of the molecules as well.Now, P can characterize the total coupling strength if we choose, for instance, λ4 = ( The calculation of the energy difference between left and right states δV tot follows as in Appendix B. In the lowest order of 1/N , we find where α = 2 Pβ E b now functions as a control parameter.It is noteworthy that the nonlinearity of the dipole moment is connected with a ∆ 2 -term in the energy difference.This term shows that the different light-matter coupling constants generate an effective energy gradient to the molecular states.Furthermore, the collective coupling always vanishes if the dipole moment is fully symmetric so that λ 2 jR = −λ 2 j L .Here, we concentrate on a numerical example, even though one could proceed with a similar kind of argumentation as in Sec.B. Let us choose a light-matter coupling that vanishes on one state and is finite on the other.How does the polaritonic system relax, if all the molecules are initially in the left state, and what is the stationary state?To further connect this setup to Fig. 3 in the main text, we choose λ4 so that it is the larger one of λ 4 R and λ 4  L , choose the couplings from a normal distribution with Var(g) = 0.25 〈g〉 2 (in the description λ 2 i = ω m ω c g i ), and set α = 1.3.That is, in the case in which only the molecules on the left state couple to light, the situation is initially the same as in Fig. 3 (orange line).
The result of the numerical simulation is plotted in Fig. 6.As expected, the initial evolution of the macroscopic state is very similar in the case where the left state couples to light as to when both states are coupled to light.The stationary state changes only in a minor way.
Much more interesting is the case where the coupling to light is on the right state only.Initially, there is no light-matter coupling and the polaritonic system starts to thermalize towards an even split between the states.This initial rate is similar to that without any coupling.The stationary state however changes notably as more of the molecules prefer the right state.This again shows the effective attractive interaction mediated by the cavity.Finally, because of symmetry in the potential and the mirror symmetry of the nonlinear dipole moment, the stationary states of these two plotted scenarios are related.One gets from one to the other by simply interchanging the labels for left and right.
Finally, we note that Eq. (F.2) allows solving the stationary states analytically.This is the most straightforward in the case in which there is no variation in the left and right coupling constants, that is, Focusing still on the symmetric potential, the relevant parameter is the ratio λ 2 R /λ 2 L = g R /g L .The effect on the bifurcation diagram, comparing to Fig. 2(a) in the main text, is similar to having a bias in the potential as shown in Fig. 7.That is, the critical value of α at which the bifurcation happens increases as the ratio g R /g L deviates from unity.Similar to the effect of bias, there is a measurable effect before the critical value as more molecules are found in the state with the larger coupling.

G Realizations of cavity-induced bifurcation G.1 Polaritonic chemistry
Let us discuss the validity range of the approach in the main text in the context of polaritonic chemistry.As we derive the rate constant modification from detailed balance, our results disregard quantum effects.That is, the temperature must be above a certain threshold frequency that is related to ω c , ω L , ω R [11].This limits the systems mostly to vibrational strong coupling where these frequencies are typically of the order of 100 meV and not in the range of several eV as in electronic strong coupling (at room temperature k B T ∼ 25 meV).More notably, we treat the molecules as one-dimensional systems and assume a constant bilinear coupling to the (single mode) cavity.We believe that even though this choice neglects many details of molecules used in recent experiments, it gives an important analytical insight which can be improved upon with specific molecular models.
We note that the potential V (q) with σ = 0 and ω L = ω R is similar to the ground state potential of the Shin-Metiu model describing proton-coupled electron transfer [58].Recently, it has been used in other works in polaritonic chemistry [40,41].In these works, the typical barrier energy is of the order of 1 eV; although it depends on the exact choice of model parameters.Since the effective potential barrier E b is typically larger than the barrier energy, our order of magnitude estimate is that E b /(k B T ) may range from tens to hundreds.
To identify the strong coupling, the Rabi splitting Ω must be large enough to be observed.This means that Ω ≳ (κ + γ)/2 where κ and γ refer to the dissipation rates of the cavity and the molecular mode, respectively.We have assumed the collective strong coupling regime where Ω/ω m ∼ 0.1 leading to P ∼ 0.01.Thus, all the results are given in the first order of P and higher order corrections are yet to be found.Experimentally, it has been shown that it is possible to reach the so-called ultrastrong coupling regime in vibrational strong coupling; Ω/ω m ≈ 0.24 in Ref. [59].
Therefore, the control parameter α = 2Pβ E b can be of the order of unity in experiments of polaritonic chemistry.For example, if P = 0.01 and β E b = 100, we have α = 2.In the main text, we find that the cavity-induced phase transition happens when α > 1 for systems with two metastable states at the same energy.
We emphasize the fact that the coupling constants can differ from molecule to molecule which, somewhat surprisingly, is often neglected in the seminal theoretical works [46].To further complicate the picture, the couplings may truly be time dependent if the molecules can diffuse within the cavity.However, such time dependent variation is much faster than the typical reaction rates, and thus, we can deal with only the averages of coupling constants.There are at least two types of randomness in the couplings: 1) "orientation disorder" due to the varying directions of transition dipole moments with respect to the cavity polarization vector(s), and 2) "position disorder" due to the varying positions of the molecules within the cavity.These disorders have been shown and studied in several experiments, for instance, in plasmonic picocavities [44] and in optical cavities [45].In the main text, we find that the average of the coupling constants determines the macroscopic properties of the N molecule polaritonic system, for instance the stationary state and the rate modification.The detailed understanding of these disorders is thus important.
It is speculated in Ref. [12] that the orientation disorder must be negligible in order to see polaritonic chemistry.This indeed corresponds to the rate modification analysis in the main text, as assuming that the transition dipole moments are distributed isotropically would lead to the average of the coupling constants being zero and leaving the rate unchanged.However, it is not clear to us what processes would cause the alignment of molecules in the experiments.There are fabrication techniques to construct such aligned molecular assemblies [60,61] but they represent a departure from the microfluidistic cavity approach detailed in Ref. [12].We hope that our work could motivate theoretical and experimental inquiries into this kind of disorder.
The effect of the positions of the molecules within the cavity has been investigated experimentally [45].The experimental findings agree well with the current theoretical understanding presented here and in the main text.That is, local variations in the cavity field strength affect the couplings.We also know from Maxwell's equations that such vacuum electromagnetic fields may change sign as a function of position within the cavity.This is often neglected since observables such as the energy and the absorption/fluorescence spectrum are independent of the sign.Also, for a single molecule, the sign is always an irrelevant phase factor of the cavity field.The phase difference between different positions is important for the chemistry because one would expect that the molecules are distributed somewhat uniformly within the cavity.The average of the coupling constants should be zero for antisymmetric cavity modes and decreases as the order of the cavity mode is increased for symmetric modes.Thus, it would seem beneficial to limit the variation in the molecule positions by cavity design and/or using the lowest order cavity mode.The reader should be aware that this message is in contradiction with the published experiments in polaritonic chemistry: For instance, the second-order cavity mode was put in resonance with a vibrational mode in Ref. [17] and the tenth-order mode was used in Ref. [15].It is an intriguing possibility that some effects would be, in fact, caused by a coupling to the highly detuned first cavity mode and, consequently, the interpretation of polaritonic chemistry as a resonant effect would be challenged.
In the presence of great position and orientation disorder, we expect no change in the rates.Nevertheless, if Rabi splitting is large enough and α > 1, a phase separation based on the value of coupling constants is possible.If we imagine a case where we have somehow aligned the transition dipole moments and the relevant cavity mode is the second lowest order mode, the phase separation would manifest as a physical separation.The molecules on one side of the cavity would then go to another state than those on the other side.This could be achieved, alternatively, with two films of oriented dipole moments using the same measurement setup as in Ref. [45] so that the films are at the opposite sides of the cavity (that is, mirror-film-spacerfilm-mirror).

G.2 Coupled anharmonic LC oscillators
Although our work was initially motivated by polaritonic chemistry, our generic stochastic model is applicable to a large variety of systems.Here we illustrate how similar physics could be studied in an ensemble of superconducting flux qubits all coupled to the same cavity.We also show how in this setup one can in principle realize the alternating sign of the couplings.
Consider the electronic circuit drawn n Fig. 8.It describes an LC circuit coupled via mutual inductance with N (ideally identical) flux qubits [62].The figure depicts a simplified version of the flux qubit, but the idea can be generalized to the usually used configurations where the inductance L q is replaced by several Josephson junctions.Such systems have two macroscopic quantum states corresponding to clockwise and counterclockwise persistent currents, used as the logical qubit states.When the flux Φ ≈ h/2e[1 + f /(2π)], the flux qubit Hamiltonian is Here Qq is the displacement charge operator of the capacitor with capacitance C q .Qq is canonically conjugate with the phase difference φ across the junction.They thus take the roles of canonical momentum and position in the electronic circuit.E J is the Josephson energy of the junction and b = 2π 2 h 2 /(2e L E J ) is a dimensionless parameter characterizing the inductor with inductance L. For f ≪ π and b < 0.5, the effective potential V (φ) has two minima corresponding to the two directions of the persistent current.The relative bias between the minima can be controlled with f such that the two minima are degenerate for f = 0.For small b and f , the minima are around φ ± ≡ ±π 1+2b + f , and both have eigenfrequencies ω L = ω R ≈ ω p (1 + 2b), where ω p = 2e E J /C q /ħ h is the Josephson plasma frequency.
Then, the effective barrier height is Flux qubits operate in the regime where the potential barrier between the states is low, E b ∼ ħ hω L/R , and thereby the tunnel coupling is large, providing coherent oscillations of the quantum state between the two minima.This limit is reached when the charging energy E C = e 2 /(2C q ) is not much smaller than the Josephson energy E J .On the contrary, we assume E J ≫ E C and thereby a high barrier, such that the transitions between the minima are rare and mostly driven by thermal noise in the environment of the system.
In such systems, the mutual inductance coupling of the individual qubits to a common LC circuit has been used as a means to realize controllable coupling between the qubits, as they are tuned on and off resonance with each other [62].However, there the resonance condition derives from the resonance of the qubit energy splitting, whereas we assume that the resonance is between the bistable energy minima.Let us then consider a mutual inductance coupling of the flux qubits to a common LC circuit.We model this coupling in the limit where the qubits are in one of their minima, in which case we describe them as LC oscillators.The Lagrangian of the system is thus where we have introduced capacitance C p and inductance L p such that ω p = 1/ L p C p for the flux qubits.Note that this form of expression is possible because ω L = ω R = ω p ; otherwise the state of the flux qubits affects the capacitance and inductance.The mutual inductance M j can be characterized by a dimensionless parameter k j ∈ [−1, 1] so that M j = k j L L p .Here, one can understand the possibility of a negative k j as a phase shift of the induced current; the sign depends on the details of the inductive coupling.We define the conjugate momentum using the notation of the flux quantum Φ 0 = h/2e as which can be also formally written as ⃗ ϕ ∝ L ⃗ Q.With this definition, the left hand side represents magnetic flux.The Hamiltonian of the full systems is obtained by the Legendre transform.We find The inversion of the inductance matrix L presents a difficult analytical problem.If the mutual inductances M j are large, that is k j ≈ 1, every element of the matrix L −1 is generally nonzero.However, we assume that they are small.To first order, the solution is obtained by noting that we can obtain the conjugate relation for ϕ only in the zeroth order of M 's and substitute the result in the original Lagrangian L in Eq. (G.2) as it is first order in M 's.We find for an LC circuit of resonant frequency ω c = 1/ LC C p ϕ j .Thus, in the lowest order approximation, the Hamiltonian corresponds to the potential we use in Appendix A with λ 2 j = k j ω c ω p .Note that both of these models are within the harmonic approximation, as we expanded the flux qubit Hamiltonian around a minimum point.
In other words, the cavity induced bifurcation could be observed in flux qubits by studying their collective evolution towards a state where even in the absence of a bias f , a majority of the qubits would be in the state with, say, clockwise circulating persistent current.

G.3 Cold atoms
Recently, there has notable progress in the field of cold polar molecules [63].Following the long tradition of optical trapping of cold molecules, this field provides an interesting opportunity to realize a controllable polaritonic system.This is because the many internal degrees of freedom in molecules allow for metastability in the energy, as described in this supplement and in the main text.Thus, in the future, it would seem possible to manufacture a system that is topologically the same as in polaritonics but with controllable coupling constants, and then investigate the noise activated processes or even quantum tunneling between the possible molecular states.

j
µν = Γ 0 µν r j µν (N L ) where r j µν (N L ) represents the rate modification caused by the vacuum field which may be expected to be a function of N L and not of the individual p j 's.The subsequent solution of the rate modifications r j µν is the main result of this work from which the cavity-induced collective effects follow.

Figure 2 :
Figure 2: (a) Steady state solutions of Eq. (2) exhibiting bifurcation for equal couplings.The dotted line represents an unstable solution whereas the solid lines are stable; the average tends to one of the stable solutions after a long time.(b) Phase separation of a polaritonic system that is divided into two equally large partitions (N 1 = N 2 = N /2) with couplings related by g 1 = −2g 2 .The solid lines represent the dynamics of the partitions alone whereas the dashed line describes the dynamics of the full system.Here, the potential V (q) is symmetric with Γ 0 LR = Γ 0 RL and σ = 0.The system has undergone the phase transition as α = 1.3.

Figure 3 :
Figure 3: Kinematics for different Gaussian distributions of coupling constants g j .The other parameters are as in Fig.2(b).The solid lines are obtained by the Markov chain simulation of 1000 molecules averaged over 100 realizations, the dashed line by evaluating Eq. (2), while the dotted black lines represent both the stationary state [Eqs.(11)-(12)] and the initial behaviour [Eq.(14)].Inset: The rate modification in the initial state N L = N [Eq.(14)] as a function of α.From bottom to top, the different lines are obtained by setting 〈g〉 2 = c g 2 with c = 1, 0.75, 0.5, 0.25, 0.

2 j 〈λ 4 〉
One can derive a self-consistency equation for ∆ in the stationary state.Let us assume that, even though there are N systems, there are M < N different values for the coupling constants λ 2 j .Let us denote c j = λ and α = 2Pβ E b here, for brevity.We divide N into M subsystems or partitions so that N = M k=1 N k .Each subsystem k ∈ {1, 2, . . .M } has N k systems which all share the same coupling constant c j = c k .Then, ∆ = 1

Figure 5 :
Figure 5: A single realization of the Markov chain simulation for N = 1000 systems and 1000 time steps (Γ 0 ∆t = 0.01).Other parameters are α = 1.3, σ = 0, andΓ 0 LR = Γ 0 RL = Γ 0 .Here, the coupling constants are sampled from a Gaussian distribution with zero mean.a) The time evolution of the quantities δN and ∆.As discussed in the main text, ∆ has a solution other than zero when α > 1 while δN fluctuates around zero. b) As a visual representation, the states "left" and "right" are here represented by blue and orange, respectively, for each time step.Here, the different systems are indexed in the ascending order in the coupling constant (i.e., g 1 ≤ g j ≤ g N ).The correlation becomes visible; the systems with largest values of coupling are predominantly found in the left well after a transient period.

Figure 6 :
Figure 6: Monte Carlo simulation of a polaritonic system with nonlinear dipole moments, schematized on the right side of the figure.The transition dipole moment vanishes on the right state in the blue line and on the left state in the brown line.The orange line from Fig. 3 is the dashed line here for reference.The simulation uses N = 1000 molecules averaged over 100 realizations, and Γ 0 LR = Γ 0 RL = Γ 0 .

Figure 7 :
Figure 7: Bifurcation diagram for a symmetric potential without any variation in the left and right light-matter couplings.Here, we fix λ4 = λ 4 L within α and the different curves represent different values of λ 2 R .

Figure 8 :
Figure 8: Electronic circuit model of an LC oscillator coupled to a set of flux qubits, where the Josephson junctions are marked with the triangle.There is a flux Φ = h/(2e) + δΦ applied across each loop.

0 2π
X j + ω c ω p j k j X c X j , (G.5)where, in the last step, we denoted P j = Q j / C p and P c = Q c / C and correspondinglyX c = Φ Cϕ c and X j = Φ 0 2π